2D Deconvolution using FFT in Matlab Problems

Question

I have convoluted an image I created in matlab with a 2D Gaussian function which I have also defined in matlab and now I am trying to deconvolve the resultant matrix to see if I get the 2D Gaussian function back using the fft2 and ifft2 commands. However the matrix I get as a result is incorrect (to my knowledge). Here is the code for what I have done thus far:

% Code for input image (img) [300x300 array]

N = 100;
t = linspace(0,2*pi,50);
r = (N-10)/2;
circle = poly2mask(r*cos(t)+N/2+0.5, r*sin(t)+N/2+0.5,N,N);
img = repmat(circle,3,3);

% Code for 2D Gaussian Function with c = 0 sig = 1/64 (Z) [300x300 array]

x = linspace(-3,3,300);
y = x';
[X Y] = meshgrid(x,y);
Z = exp(-((X.^2)+(Y.^2))/(2*1/64));

% Code for 2D Convolution of img with Z (C) [599x599 array]

C = conv2(img,Z);

% I have tested that this convolution is correct using cross section profile vectors for img and C and the resulting x-y plots are what i expect from the convolution.

% From my knowledge of convolution, the algorithm works as a multiplier in Fourier space, therefore by dividing the Fourier transform of my output (convoluted image) by my input (img) I should get back the point spread function (Z - 2D Gaussian function) after the inverse Fourier transform is applied to this result by division.

% Code for attempted 2D deconvolution

Fimg = fft2(img,599,599);

% zero padding added to increase result to 599x599 array

FC = fft2(C);
R = FC/Fimg;

% I now get this error prompt: Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.551432e-22

iFR = ifft2(R);

I'm expecting iFR to be close to Z but I'm getting something completely different. It may be an approximation of Z with complex values but I can't seem to check it since I don't know how to plot a 3D complex matrix in matlab. So if anyone can tell me whether my answer is correct or incorrect and how to get this deconvolution to work? I'd be much appreciated.

Menado

Slightly off topic nitpick: you convolve a signal, not convolute – Porta

Since the division results in a 599x599 matrix, how do you get the 2D Gaussian Function which is 300x300? – Rb

@Rb the resultant image is 599x599 with the 2D Gaussian contained within the first 300x300 rows. I used matrix indexing to get back a 300x300 array – Menado

Answer 1

R = FC/Fimg needs to be R = FC./Fimg; You need to do division element-wise.

Alignment

Of course, this makes a lot more sense; what a rookie mistake. Is there a way to plot the resultant complex matrix as a surface plot so I can compare this with my gaussian? – Menado

@SeanJamesJamieson I think it there shouldn't be any complex components... Not sure though. But just incase you can wrap the array with an abs before plotting it. – Alignment

Yeah I did this and I indexed the matrix for the first 300 rows and columns and I recovered my function back with an maximum pointwise defined uncertainty of 4 x 10^-13. Thanks for the help guys – Menado

Nice, but how do you know that the first 300 rows and columns represent the Gaussian function? For example, what happens if you take the last 300 rows and columns? I am doing something similar and this part is confusing me :/ – Rb

@Rb Can you elaborate what you mean? – Alignment

I was referring to @SeanJamesJamieson last step in getting the Gaussian function. In his words, "I indexed the matrix for the first 300 rows and columns and I recovered my function." The division results in a 599x599 matrix, but how do we get the 300x300 matrix representing the Gaussian function? – Rb

@Rb the recovered Gaussian is contained within the first quadrant of the resultant array after deconvolution. I think its to do with the default way matlab pads the initial array since now when i use a vector input in the padarray function to add an equal specific amount of zeros before and after the rows and columns of the array, The Gaussian is contained within the centre of the array – Menado

Answer 2

Here are some Octave (version 3.6.2) plots of that deconvolved Gaussian.

% deconvolve in frequency domain
Fimg = fft2(img,599,599);
FC = fft2(C);
R = FC ./ Fimg;
r = ifft2(R);

% show deconvolved Gaussian
figure(1);
subplot(2,3,1), imshow(img), title('image');
subplot(2,3,2), imshow(Z), title('Gaussian');
subplot(2,3,3), imshow(C), title('image blurred by Gaussian');
subplot(2,3,4), mesh(X,Y,Z), title('initial Gaussian');
subplot(2,3,5), imagesc(real(r(1:300,1:300))), colormap gray, title('deconvolved Gaussian');
subplot(2,3,6), mesh(X,Y,real(r(1:300,1:300))), title('deconvolved Gaussian');

% show difference between Gaussian and deconvolved Gaussian
figure(2);
gdiff = Z - real(r(1:300,1:300));
imagesc(gdiff), colorbar, colormap gray, title('difference between initial Gaussian and deconvolved Guassian');

Finned