% Program to demonstrate inverse filtering
% inversefilter_new.m

% need to be in directory d:\jl\4F8_imageprocessing\imageprocessing5\imageproc5mat

clear,clc,close all,home

% Read image from disc and gaussian
load lenna256.mat;


xx=x(65:192,65:192);
[a1,a2]=size(xx)
% take the central 128x128 bit of Lenna
xnew=zeros(192,192);
% embed it in a 192x192 matrix of zeros -- i.e. zero pad the edges
xnew(33:160,33:160)=xx;


% uncomment the lines below to show what happens when we do not put a guard
% band of zeros around the edges of the image
% xx = x(33:224,33:224);
% xnew = xx;


[n1,n2]=size(xnew)
figure(1)
grayimage(xnew);
title('zero padded original image');
pause

% now convolve this with a gaussian -- no noise --
% and deconvolve using inverse filtering

%disp('Use 192 data points and variance=10');
N=input('Number of samples in horizontal and vertical direction (even), N = ');
var=input('Variance of gaussian = ');
N1=N/2;
vec1=(-N1:N1-1)';
for i=1:N
  for j=1:N
    g(i,j) = 256*exp(-((i-(N1+1))^2+(j-(N1+1))^2)/var);
  end
end


%  now convolve the gaussian and Lenna

figure(2)
grayimage(g);
title('gaussian blurring function');
pause;

Y = conv2(xnew,g,'same');
[m1,m2]=size(Y)
% now add noise: first 0, then 0.01, then 0.1, then 0.5

noise=input('Noise = ');

for i=1:m1
  for j=1:m2
    if Y(i,j)>0
      Ytmp(i,j)= Y(i,j) + noise*randn;
    else
      Ytmp(i,j)=0;
    end
  end
end

Y=Ytmp;
figure(3);
grayimage(Y);
title('blurred image (plus noise if any added)');
pause

% shows the blurring produced by this particular convolution -- in
% this case the gaussian is an example of a point spread function

% form the FFT of convolved Lenna and the convolving function

% let threshold be gamma=0.001
gamma=0.001;
YFT = fft2(fftshift(Y));

XFT = fft2(fftshift(xnew));

GFT = fft2(fftshift(g));

for i=1:N
  for j=1:N
    if (abs(GFT(i,j))) > gamma
      GFTinv(i,j)=1/GFT(i,j);
    else
      GFTinv(i,j)=0;
    end
  end
end
size(YFT)
size(GFTinv)       
figure(4)
grayimage(abs(GFT));
title('FT of convolving function');
pause;
figure(5)
grayimage(abs(GFTinv));
title('Inverse of FT of convolving function');
pause;


XFTnew = YFT.*GFTinv;

Xnew = fftshift(ifft2((XFTnew)));

figure(6);
grayimage(abs(Xnew));
title('image restored via inverse filter');
pause;


