Deconvolution of an image with a gaussian point spread function.

32 views (last 30 days)
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')

Answers (1)

Image Analyst
Image Analyst on 2 Oct 2022
Maybe try some of the numerous built-in deconvolution methods, such as deconvlucy

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by