%% Using the FT to deblur an image
clear
close all

% Load in an image
img = im2double(imread('cameraman.tif'));
figure
imshow(img)

% We will take this as ground truth at the sample.  If we image this on the 
% microscope, it will get blurred by the PSF. 

% To 
IMG=fft2(img); 
IMG= fftshift(IMG); 

figure
imshow(log(abs(IMG)),[])


% Construct Gaussian PSF.  It has a PSF with radius of 2 pixels.
filt = fspecial('Gaussian',13,2);
figure
imshow(filt,[])

% We can use the command psf2otf (point spread function to optical transfer
% function) to turn our blur into the frequency domain.
FILT = psf2otf(filt,size(img)); 
FILT = fftshift(FILT);

figure
imshow(abs(FILT),[])

size(FILT)
size(IMG)

% To get our blurred image, we just have to take multiply the FFT by the OTF 
imgBlur = ifft2(IMG.* FILT); 
imgBlur = abs(imgBlur); 

figure
imshow([img imgBlur],[])

%% Now we have a blurred image
%Generate blurred image via Fourier domain
IMGBLUR = fft2(imgBlur); 
IMGBLUR = fftshift(IMGBLUR); 

% We can display the two power spectra to see what they look like.
figure
subplot(1,2,1)
imshow(log(abs(IMG)),[])
title('Initial')
subplot(1,2,2)
imshow(log(abs(IMGBLUR)),[])
title('Blurred')

% Now lets divide the blurred image FFT by the FILT FFT.
IMGdeBLUR = IMGBLUR./FILT;

% Now to reconstruct our image, we just have to take the ifft, and then we
% only want the real parts of it.
imgdeblur = abs(ifft2(IMGdeBLUR));

figure
imshow([img imgBlur imgdeblur],[])


%% What happens if we try with a wrong filter

% We can make a new Gaussian filter with SD of 1.9 instead of 2.
filt2 = fspecial('Gaussian',13,1.8);
FILT2 = psf2otf(filt2,size(img)); 
FILT2 = fftshift(FILT2);

% Now we deblur exaclty the same way
IMGdeBLUR2 = IMGBLUR./FILT2;
imgdeblur2 = abs(ifft2(IMGdeBLUR2));
figure
imshow([imgdeblur imgBlur imgdeblur2],[])
