% 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 Image Spectra');
%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 in Fourier Space');

% 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 (Inverted) Low Pass Filtered Image = Original');