% This used the same script as for butterworth2.m example
% Low Pass Filter to create a blurred image
%
% Simply do deconvoution of butterworth low pass (divide rather multiply)
% to invert the result
%
% Note in practice the ideal forward blurring function will not be known 
% SO results will not be so perfect.


% read in MATLAB demo text image
image = imread('text.png');
[M N] = size(image)

% Show Image

figure(1)
imshow(image)
title('Original Image')

% compute fft and display its spectra

F=fft2(double(image));
figure(2)
imshow(abs(fftshift(F))/256);
title('FFT Spectra of Image');

%compute Ideal Butterworth Pass Filter

u0 = 50; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);

for i = 1: M
    for j = 1:N
      %Apply a 2nd order Butterworth  
      UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
      
      H(i,j) = 1/(1 + UVw*UVw);
    end
end    
% display
figure(3)
imshow(fftshift(H))
title('Butterworth Low Pass Filter');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure(4)
imshow(g)
title('Butterworth Low Pass Filtered Image');


%%%%%% NOW ILLUSTRATE DECONVOLUTION   %%%%
% Transfer (blurring function) is H --- the low pass filter
% Simply do deconvoution of butterworth low pass (divide rather multiply)
% to invert the result --- effectively now a high pass filter.


Ghigh = G./H;
ghigh=real(ifft2(double(Ghigh)));
figure(5)
imshow(ghigh)
title('Deconvolved Image');

%%% Demo now to show noisy enhancing properties of the High pass filter


% Filter using High pass filter H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure(6)
imshow(g_highpass);
title('Origial Image Deconvolved');



% High Pass Filter Original Image  Strictly a High Pass Filter is 1 - H 

Fhigh = F.*(1-H);

fhigh=real(ifft2(double(Fhigh)));
figure(7)
imshow(fhigh)
title('Original image High Pass Filtered');

% now addnoise to the image

image_noisy = imnoise(double(image),'gaussian');

% show noisy im
figure(8)
imshow(image_noisy)
title('Image with Noise Added');
figure(8)

%%%% compute FFT of noisy image

F = fft2(double(image_noisy));

% Filter using high pass filter H as above

G_highpass=F.*(1 - H);

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure(9)
imshow(g_highpass);
title('High Pass Filter Noisy Image');
figure(9)

% Filter using high pass filter H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure(10)
imshow(g_highpass);
title('Deconvolved Noisy Image');
figure(10)



% now addnoise to the blurred low pass image g and do the same

image_noisy = imnoise(double(g),'gaussian');

% show noisy im
figure(11)
imshow(image_noisy)
title('Low Pass Filtered (Blurred) Image with Noise')
figure(11)

%%%% compute FFT of noisy image

F = fft2(double(image_noisy));



% Filter using deconvolution H as above

G_highpass=F./H;

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure(12)
imshow(g_highpass);
title('Deconvolved Blurred Noisy Image');
figure(12)


% Filter using High Pass Filter (1 - H) as above

G_highpass=F.*(1 - H);

% convert back to image spce with iFFT

g_highpass = real(ifft2(double(G_highpass)));

figure(12)
imshow(g_highpass);
title('High Pass Filtered Blurred Noisy Image');
figure(12)