Deblurring an Image using inverse filtering

I am trying to deblur an image using inverse filtering that was blurred using a 25x25 gaussian blur function with sigma = 15. I am extracting the blurred image from a .mat file, displaying it which works correctly.

Next I define my gaussian filter and then compute frequency reponse of the filter. To deblur the image, I divide blurred image by frequency response of the filter and take ifft.

The blurred image displays correctly in figure 1, but figure 2 which should display deblurred image displays all purple. I am trying to keep my code as simple and minimal as possible.

What I am doing wrong here? I will appreciate any hints or inputs

images = load('project_images.mat');   % Load the mat file containing images

m_blur = images.mandrill_blurred;  % Extract the first image 
imagesc(m_blur);                   % display the blurred image

h = fspecial('gaussian',[25 25],15);   % 25x25 Gaussian blur function with sigma = 15
hf = fft2(h,size(m_blur,1),size(m_blur,2));      


m_deblur = real(ifft2(m_blur)./hf);      %inverse filter 

figure(2)
imagesc(m_deblur)                        % Display deblurred image

[Rik]

Please attach your data.

[sid101]

Not uploading project_images.mat since it is a big file.

I uploaded the mat file for the blurred image

[RAGNAR SA]

Are still have the final code for your question? after making the changes with help of Rik and the image analyst.

Answer 1

[Rik]

A bit of data exploration shows that you have quite an outlier in your image:

figure(3),clf(3)
histogram(m_deblur)
set(gca,'YScale','log')
axis([-10 140 0.1 max(ylim)])

Once you replace that with a 0, the automatic scaling should work as expected again. In the code below I went a bit further and set the caxis value manually to something that felt about right.

figure(2)
imagesc(m_deblur)                        % Display deblurred image
caxis([-0.15 0.15])

So in conlusion: this image is not ready yet.

The reason for this is that you didn't put the blurred image in the Fourier domain yet, so the division doesn't make a lot of sense.

images = load('mandrill_blurred.mat');   % Load the mat file containing images

m_blur = images.mandrill_blurred;  % Extract the first image 
imagesc(m_blur);                   % display the blurred image

h = fspecial('gaussian',[25 25],15);   % 25x25 Gaussian blur function with sigma = 15
hf = fft2(h,size(m_blur,1),size(m_blur,2));      

m_deblur = real(ifft2(fft2(m_blur)./hf));      %inverse filter 
%                     ^^^^^      ^    ^
%                      you forgot this

figure(2)
imagesc(m_deblur)                        % Display deblurred image

Answer 2

[Image Analyst]

You have to do the division in Fourier space. You're doing it in the spatial domain, AFTER you've converted your image back to the spatial domain with ifft2().

m_deblur = real(ifft2(m_blur)./hf);

So if m_blur is the spectrum, then ifft2(m_blur) is the spatial domain image. Then hf is the spectrum of h, so in the line above you're dividing a spatial domain array by a frequency domain array. That doesn't make sense. You need to do it like this (adapt names as needed):

freqImage = fft2(spaceDomainImage);
freqFilter = fft2(spaceDomainFilter);
inverseFilteredImage = freqImage ./ freqFilter; % This is still in the frequency domain.
spatialFilteredImage = ifft2(inverseFilteredImage); % Convert it back to spatial domain.

[sid101]

Thanks both answers say the same thing and work. I wish I could accept two answers.

[Ash4Matlab]

Hello sid101 and Image Analyst,

I used a Point Spread function (PSF) to blur my image. Now I am trying to deblur image by dividing my blurred image with the PSF. But I am not getting anything meaningful as you can see below

I have added the code here too. data.mat contains the the image and compex binary pupil used for blurring.

It will be helpful if you can teel me the issue here. Thanks in advance.

load data.mat                                % Contains Image and Pupil (h)
ImageFT    = fftshift(fft2(Image));          % Fourier Transform of the Image
PSF        = abs(h).^2;                      % Point spread function (PSF)
PSF        = PSF/sum(PSF(:));                % Normalizing PSF
PSFFT      = fftshift(fft2(ifftshift(PSF))); % Fouirer Transform of PSF
Blurred    = ImageFT.*PSFFT;                 % Blurred Image in frequecny domain
%% Inverse Filtering
Deblurred = (Blurred./PSFFT);                % Deblurring in frequency domain
%% Show Images
figure(1)
subplot(1,3,1)
imshow(Image)
title('Original Image')
subplot(1,3,2)
imshow(mat2gray(real(ifft2(ifftshift(Blurred))))) % Blurred image in spatial domain
title('Blurred Image')
subplot(1,3,3)
imshow(mat2gray(real(ifft2(ifftshift(Deblurred))))) % Deblurred image in spatial domain
title('Deblurred Image')

[Image Analyst]

There is no way that you can do an inverse filter to return that extremely blurred image to anything close to the original.