% dip08_5.m           gr_nov16
%    Motion blur and turbulence blur restoration with inverse filter and Wiener filter
clear all;

% f = double(imread('.\img\lenar.tif'));
 f = double(imread('.\img\lena.tif'));  f = f(101:450,51:400);
[m,n] = size(f);

[u,v] = meshgrid(1:n,1:m);
m2 = m/2+1;  n2 = n/2+1;

        %%%%%%%%%%%%%%%%% Atmospheric turbulence, k = .0002 .. .002
 PSF = exp(-.0002*((u-m2).^2+(v-n2).^2).^(5/6));
 
        %%%%%%%%%%%%%%% Motion degradation model
a = 1/8;     % average of a*m shots, each with a one-pixel shift
              % for m=256, try with a = 1/128 ... 1/32 ... 1/8
%  PSF = sin(pi*(u-m2)*a) .* exp(-j*pi*(u-m2)*a) ./ (pi*(u-m2)*a);   
%  PSF(:,m2)=1;   %  when u=m2 the PSF is sin(0)/0, a Matlab NaN

figure(1);  mesh(abs(PSF));

F = fftshift(fft2(f));
Fb = F .* PSF;
fb = real(ifft2(ifftshift(Fb)));      %  blurred image
sdn = 0.0000000000001;   fb = fb + randn(m,n)*sdn;  %  ...plus noise
% fb = double(uint8(fb));  % quantization noise

G1 = Fb ./ PSF;      %  not realistic, only for reference
G2 = fftshift(fft2(fb)) ./ PSF;      %  inverse filter
H = PSF;  H(abs(PSF)<.01)=.01;  % 'raise' zeros or very small values in the PSF
G3 = fftshift(fft2(fb)) ./ H;       %  clipped inverse filter

g1 = real(ifft2(ifftshift(G1)));
g2 = real(ifft2(ifftshift(G2)));
g3 = real(ifft2(ifftshift(G3)));
figure(2);  imshow([fb, g1; g2, g3],[0 255]);
return;

%%  Wiener filter
Hw = (conj(PSF) .* F.^2) ./ (PSF.^2 .* F.^2 + sdn^2); 
% Hw(abs(Hw)>5)=5;  % clip peaks in the Wiener filter
 figure(3);  mesh(abs(Hw)); 
Gw = fftshift(fft2(fb-mean(fb(:)))) .* Hw;
gw = real(ifft2(ifftshift(Gw))) + mean(fb(:));
figure(4);  imshow([fb, gw],[0 255]); 

%%  Matlab's Wiener filter
psf = real(fftshift(ifft2(ifftshift(PSF))));
nsr = sdn.^2 ./ sum( (f(:)/255).^2);  % noise-to-signal power ratio
gwm = deconvwnr(fb/255,psf,nsr);  % Matlab's Wiener filter
% figure(5);  mesh(abs(psf)); 
 figure(6);  imshow(gwm*255,[0 255]); 

return;
