% This used the same script as for butterworth.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.


% Create a white box on a black background image
M= 256; N = 256;
image = zeros(M,N);
box = ones(64,64);
%box at centre
image(97:160,97:160) = box;

% 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)));
title('Fourier Spectra of Image');

%compute Butterworth Low Pass Filter

u0 = 20; % 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('FFT 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('Low Pass 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(6)
imshow(fhigh)
title('Original image High Pass Filtered');

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


% now addnoise to the image

image_noisy = imnoise(image,'gaussian');

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

%%%% 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(8)
imshow(g_highpass);
title('High Pass Filter Noisy Image')
figure(8)


% 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(9)
imshow(g_highpass);
title('Deconvolved Noisy Image')
figure(9)


