Наивная деконволюция с матрицей идентичности и изображением (без шума)

Question

Фон

Я пытаюсь сравнить шумную и бесшумную наивную деконволюцию, и я как-то наткнулся на стену на этапе деконволюции без шума, который, как мне кажется, очень специфичен, и у меня нет Идея, почему это происходит. По сути, у меня есть входное изображение 256x256, x, и я хочу свернуть его с фильтром h, то есть единичной матрицей размером 21x21.

Задача

В шаг деконволюции:

X = Y / H

Я, очевидно, столкнусь с проблемами, потому что H содержит нули, а мой x будет иметь значения nan. Поэтому я попытался добавить к H очень маленький термин, например 1e-20. По какой-то причине работает только 1 * 10 ^ -9,5. Я не уверен, что это мой код или как работает MATLAB.

Код V1.0

%% Naive Deconv
    x = imread('C:/Users/chang/OneDrive/Desktop/UCLA/Winter 2020/EE 211/BSDS300-images/BSDS300/images/train/94079.jpg');
    figure(10), imshow(x)
    disp(size(x))
    targetSize = [256 256];
    rect = centerCropWindow2d(size(x),targetSize);
    img_crop = imcrop(x,rect);
    figure(20), imshow(img_crop)


    x = rgb2gray((img_crop));
    %x = mat2gray(x);
    x = im2double(x);
    figure(30), imshow(x)

    n = 21;
    h = eye(n) + 1*10^(-9.5);

    % h = exp(h)/sum(exp(h))
    h = h./n;

    X = fft2(x, 276, 276);
    H = fft2(h, 276, 276);

    Y = X.*H;
    y = ifft2(Y);

    disp(size(y))
    targetSize = [276, 276];
    rect = centerCropWindow2d(size(y),targetSize);
    y = imcrop(y,rect);
    figure(40), imshow(y)

    Y = fft2(y);
    H = fft2(h, 276, 276);

    X = Y./(H);
    x = ifft2(X);

    figure(50), imshow(x)
    x = imcrop(x,[0 0 256 256]);
    figure(60), imshow(x)

Код V 2.0

%% Naive Deconv
x = imread('C:/Users/chang/OneDrive/Desktop/UCLA/Winter 2020/EE 211/BSDS300-images/BSDS300/images/train/94079.jpg');
figure(10), imshow(x)
disp(size(x))
targetSize = [256 256];
rect = centerCropWindow2d(size(x),targetSize);
img_crop = imcrop(x,rect);
figure(20), imshow(img_crop)


x = rgb2gray((img_crop));
%x = mat2gray(x);
x = im2double(x);
figure(30), imshow(x)

n = 21;
h = eye(n) + 1*10^(-12);

%h = exp(h)/sum(exp(h))
h = h./n;

X = fft2(x, 276, 276);
H = fft2(h, 276, 276);

Y = X.*H;
y = ifft2(Y);

disp(size(y))
targetSize = [276, 276];
rect = centerCropWindow2d(size(y),targetSize);
y = imcrop(y,rect);
figure(40), imshow(y)

Y_hat = fft2(y);
H = fft2(h, 276, 276);

X_hat = Y./(H);
x_hat = ifft2(X_hat);

figure(50), imshow(x_hat)
x_hat = imcrop(x_hat,[0 0 256 256]);
figure(60), imshow(x_hat)

Результаты

Для версии V1.0 ниже h буферизуется со значениями 1 * 10 ^ -9.5 справа и h буферизируется со значениями 1 * 10 ^ -9,5 слева.

Между кодами V1.0 и V2.0 все, что я сделал, это изменил несколько переменных (для ясности при написании этого поста), и теперь в V2.0 я строго ограничен 1 * 10 ^ -12. Поэтому я считаю, что у меня может быть ошибка кодирования? Когда я не использую 1 * 10 ^ -12 в V2.0, я получаю черное изображение и следующее предупреждающее сообщение:

"Предупреждение: отображение реальной части сложного ввода."

Вопрос

Мне очень интересно, почему это происходит? Я не могу точно сказать, из-за того, как ведет себя fft или из-за ошибки в моем кодировании, но мне, очевидно, нужно провести некоторое исследование. Я просто пишу здесь, если кто-то знает, или если есть связанный вопрос, который я еще не нашел. Вот ближайший связанный пост, который я смог найти, post , но я даже не дошел до шага, на котором я в конечном итоге добавлю шум, поэтому проблемы, вероятно, будут совершенно разными.

Обновления / Следующие шаги, предпринятые до сих пор

Согласно этому сообщению, fft и нулевое заполнение , я мог бы попытаться использовать padarray для увеличения разрешения ?

Кроме того, теперь я понимаю, что это не является ключевой проблемой, поскольку никто не деконволюции таким образом, и я, вероятно, должен просто идти дальше, но я просто так растерялся, почему это происходит.

Answer 1

Обратите внимание, что в выполняемых вами вычислениях есть неточности. Числа хранятся в виде значений с плавающей запятой двойной точности, которые имеют ограниченную точность. Результат каждой операции округляется до представимого значения. Особенно в FFT, который выполняет много вычислений, эти ошибки округления складываются, приводя к очень низкому уровню шума на выходе.

Например, проверьте результат

err = x - ifft(fft(x));

err имеет значения порядка 1e-16.

Наивная деконволюция значительно усиливает этот «шум», что приводит к гораздо более сильному паттерну. Чем ближе к нулю значения в H (не h!), Тем сильнее будет этот паттерн.

Если вы поиграете с размером h и отступом для получения H, вы увидите, что один и тот же h может иметь значения, намного более близкие к 0 (или даже идентичные 0) для некоторых комбинаций, чем для других:

n = 21;    % your kernel size
p = 276;   % your padding
h = eye(n);
H = fft2(h, p, p);

При вашем выборе n и p, вы найдете min(abs(H(:))) идентичным 0, но измените p=256, и вы обнаружите, что это 0,0180, намного большее расстояние. То есть, если ваш размер FFT равен 256, вы не увидите никаких образцов артефактов, даже без добавления смещения к h. Причина в том, что одна и та же диагональная волна sin c в H дискретизируется в разных местах в зависимости от размера преобразования. Для некоторых выборок мы оказались далеко от пересечения нуля, а для некоторых выборок мы точно достигли этих пересечений нуля.

Answer 2

Предупреждение: Отображение реальной части комплексного ввода.

При работе с изображениями в области Фурье следует принять абсолютное значение обратного преобразования. Например, следующий пример ничего не делает, кроме fftshift, что соответствует сдвигу фаз (умножьте на e ^jt ) в пространственной области. Обратите внимание, что если вы не берете абсолютное значение, шаблоны формируются как в вашем результате.

clear; close all; clc
img = imread('peppers.png');
img = rgb2gray(img);
imshow(img)
title('original')

figure
Fimg = fft2(img);
FFimg = fftshift(Fimg);
IFimg = ifft2(FFimg);
imshow(uint8(IFimg))
title('inverse Fourier')

figure
aIFimg = abs(IFimg);
imshow(uint8(aIFimg))
title('absolute value of inverse Fourier')