Naive Deconvolution with Identity Matrix and Image (noiseless)

Question

Background

I'm trying to compare noisy and noiseless naive deconvolution and I kind of ran into a wall on the noiselesss deconvolution step that I feel is very particular and I have no idea why it is occurring. Essentially I have a 256x256 input image, x and I want to convolve it with a filter h, that is a identity matrix of size 21x21.

Problem

In the deconvolution step:

X = Y/H

I'll obviously run into problems because H has zeros in it and my x will have nan values. So I tried to add a very small term like 1e-20 to H. For some reason only 1*10^-9.5 works. I'm not sure if its my code or if its the way matlab works.

Code 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)

Code 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)

Results

For version V1.0, below is h buffered with values of 1*10^-9.5 on the right and h buffered with values of 1*10^-9.5 on the left.

Between Code V1.0 and V2.0 all I did was change a few variables (for clarity when writing this post) and now in V2.0 I'm strictly limited by 1*10^-12. So I believe I might have a coding error? When I don't use 1*10^-12 in V2.0 I get a black image and the following warning message:

"Warning: Displaying real part of complex input."

Question

I'm very curious to why this is occurring? I can't really tell if its because of how fft behaves or if its a bug with my coding, but I obviously need to do some research. I'm just posting here in case someone knows or if there is a related question that I have not found yet. Here's a the closest related post I could find, post, but I haven't even gotten to the step where I will eventually add noise so the problems are probably completely separate.

Updates/Next steps taken so far

According to this post, fft and zero padding, I could try to use padarray to increase resolution?

Aside I'm now realizing this isn't a crucial problem since nobody does deconvolution this way and I should probably just move on, but I'm just so confused why this happening.

score 1 · Answer 1 · answered Feb 21 '20 at 10:26

Note that there is numerical imprecision in the computations you perform. Numbers are stored as double-precision floating point values, which have a limited precision. The result of each operation is rounded to a representable value. Especially in the FFT, which does a lot of computations, these rounding errors add up, leading to a very low level of noise in the output.

For example, examine the result of

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

err has values in the order of 1e-16.

The naive deconvolution greatly enhances this "noise", leading to a much stronger pattern. The closer to zero the values in H are (not h!), the stronger this pattern will be.

If you play around with the size of h and the padding to obtain H, you'll see that the same h can have values much closer to 0 (or even identical to 0) for some combinations than for others:

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

With your choices for n and p, you'll find min(abs(H(:))) to be identical to 0, but change p=256, and you'll find it is 0.0180, much more distance. That is, if your FFT size is 256, you won't see any artifact patterns, even without adding an offset to h. The reason is that the same diagonal sinc wave in H is sampled at different locations depending on the transform size. For some samplings we happen to stay far away from the zero crossings, for some samplings we exactly hit those zero crossings.

score -1 · Answer 2 · answered Feb 21 '20 at 08:58

Warning: Displaying real part of complex input.

When dealing with images in Fourier domain, you should take the absolute value of the inverse transform. For example, the following example does nothing but fftshift, which corresponds phase shift (multiply by e^jt) in spatial domain. Note that if you don't take absolute value, patterns are formed like in your result.

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')

Naive Deconvolution with Identity Matrix and Image (noiseless)

2 Answers2