Deconvolution of an image with a gaussian point spread function.

Hello everyone

I am trying to deconvolute a signal, and i know that the convolution is was gaussian. I have tried to simply convert the resulting image signal to a 2D matrix and then apply a fourier transform. I then try to create an appropriatley sized gaussian point spread function and fourier transform that as well. Afterwards, i try to deconvolute, either just by reversing the equation in the Fourier domain, e.g. F(f*h) x (F(h))^-1 = F(f) where f is the original image, f*h the aquired signal and h the gaussian psf. I also try to add a regularization parameter, though im not sure im doing that correctly. However, even disregarding the regularization problem, it seems to me i have made a mistake beforehand, since the inverse fourier transorm i receive is simply a black image instead of at least noise. Im not asking you to fix my code for me, but it would be greatly appreciated if someone could tell where i went wrong (e.g. dimension reshaping etc.) and point me in the right direction. Below my code and the images i receive.

Additional information: The original image data is of size 100x362x3 uint8.

I also get and error code: Warning: Rank deficient, rank = 2, tol = 1.885847e-11.

> In GaussianDeconv (line 21)

This is my first time posting, please tell me if there is something i should do different or if this is appropriate. Thanks in advance!

____________________________________Plot image below___________________________________________

_______________________________________Code below______________________________________________

function GaussianDeconv(A)

name = num2str(A);                         %Get name of image file
imdataO = imread(name);
imdata = imdataO(:,:,2);                    %make the data 2d
[m,n] = size(imdata);                      %Get row and column size


y = linspace(-4,4,400);                       %Create Gaussian psf
x = linspace(-4,4,400);
[X,Y] = meshgrid(x,y);
Z1 = exp(-((X.^2)+(Y.^2))/(2*1/64));
Z = imresize(Z1,[m n]);                    %resize to be compatible with image


imdataf = real(fft2(imdata));             %fourier transform
Gaussianf = real(fftshift(fft2(Z)));      %transform and shift to middle


Deconv1f = imdataf/Gaussianf;              %Deconv without regularization 
Deconv1 = real(ifft2(Deconv1f));


B = eye(m,n);                                    %create unit matrix
lambda = 1;                                      %regularization parameter
GaussianReg = Gaussianf + (lambda*B);            %Regularize?


Deconv2f = (imdataf)/(GaussianReg);        %Deconv with regularization
Deconv2 = real(ifft2(Deconv2f));

GaussianPad = padarray(Gaussianf,[(500-size(Gaussianf,1))/2 (500-size(Gaussianf,2))/2],0,'both');  %Pad so that regularization is in middle

C = eye(m);
RegMatPad = padarray(C,[(500-size(C,1))/2 (500-size(C,2))/2],0,'both');
GaussianRegPad= GaussianPad*transpose(GaussianPad) + (lambda*RegMatPad)*transpose(lambda*RegMatPad);


GRPresize1 = GaussianRegPad(201:300,70:431);            %resize
Deconv3f = ((imdataf)/(GRPresize1));                    %try regularization
Deconv3 = real(ifft2(Deconv3f));


GRPresize2 = GaussianRegPad(201:300,201:300);                   %resize
Deconv4f = ((imdataf*transpose(Gaussianf))/(GRPresize2));       %try thikonov?
Deconv4 = real(ifft2(Deconv4f));


%________________________________plots below_________________________________________%
p = 6;
ax1a = subplot(p,3,1);               %Display data
imshow(imdataO);
title('Source image');
ax1b = subplot(p,3,3);
imshow(imdata);
title('Image data in 2d matrix form');
ax2 = subplot(p,3,4);
imshow(Z);
title('Gaussian PSF');
ax3 = subplot(p,3,5);
imshow(Gaussianf);
title('PSF fourier transformed and shifted')
ax4 = subplot(p,3,6);
imshow(imdataf);
title('Source image fourier transformed');
ax5 = subplot(p,3,7);
imshow(Deconv1f);
title('Deconv1 fourier transform')
ax6 = subplot(p,3,8);
imshow(Deconv1);
title('Deconv1 fourier retransformed')
ax7 = subplot(p,3,9);
imshow(Gaussianf);
title('PSF fourier transformed and shifted')
ax8 = subplot(p,3,10);
imshow(Deconv2f);
title('Deconv2 fourier transform')
ax9 = subplot(p,3,11);
imshow(Deconv2);
title('Deconv2 fourier retransformed')
ax10 = subplot(p,3,12);
imshow(GaussianReg);
title('Gaussian Ft regularized')
ax11 = subplot(p,3,13);
imshow(Deconv3f);
title('Deconv3 fourier transform')
ax12 = subplot(p,3,14);
imshow(Deconv3);
title('Deconv3 fourier retransformed')
ax13 = subplot(p,3,15);
imshow(GRPresize1);
title('Gaussian Ft regularized,padded and resized')
ax14 = subplot(p,3,16);
imshow(Deconv4f);
title('Deconv4 fourier transform (thikonov)')
ax15 = subplot(p,3,17);
imshow(Deconv4);
title('Deconv4 fourier retransformed (thikonov)')
ax16 = subplot(p,3,18);
imshow(GRPresize2);
title('Gaussian Ft regularized,padded and resized')

Answer 1

[Image Analyst]

Maybe try some of the numerous built-in deconvolution methods, such as deconvlucy

[Yannick Robert]

I know, the point is i want to learn how to do it this way, since it is an exercise. Should have mentioned that.

[Image Analyst]

I'll add the "homework" tag for you.

[Yannick Robert]

Thank you