% Program to demonstrate Wiener filtering

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


clear,clc,close all,home

% Read image 
load lenna256.mat;
figure(1)
grayimage(x);
pause


%%%%%% Design blurring filter%%%%%%%
% use FL=0, FU=0.25, N = 7 (fu=0.1 for lots of blur)
FL=input('Lower normalised cut-off frequency = ');
FU=input('Upper normalised cut-off frequency = ');
N=input('Number of samples in single-sided impulse response = ');

% Set up rotated Hamming window
U=N;
string1=['Rotated Hamming window, N= '...
,num2str(2*U+1),'*',num2str(2*U+1)];
string2=['Rotated Hamming window, N= '...
,num2str(2*U+1),'*',num2str(2*U+1)];
w2=zeros(2*U+1,2*U+1);
for n1=1:2*U+1
	for n2=1:2*U+1
		u=sqrt((n1-U)^2+(n2-U)^2);
		if u<=U
			w2(n1,n2)=0.54+0.46*cos(pi*u/U);
		end
	end
end

% Design circularly symmetric filter
FL=2*pi*FL;
FU=2*pi*FU;

n=(-N:N)+1e-8;
h=((n.^2)'*ones(1,2*N+1)+ones(2*N+1,1)*(n.^2)).^(1/2);
h=(FU/(2*pi)*besselj(1,FU.*h)-FL/(2*pi)*besselj(1,FL.*h))./h;
h=h.*w2;


% Apply filter to image
y=filter2(h,x);
figure(2)
grayimage(y);
title('Blurred image');
pause

% Add noise to image
energy=sum(dot(x,x));
[m,n]=size(x);
energy=energy/(m*n);
%% use SN2 = 75 for the very low noise and SN2=25 for the high noise case
%S2N=input('Signal to noise ratio (dB) = ');
var=input('input noise variance = ');
% use 0.1 for small noise and 10.0 for large noise
%S2N=10^(S2N/10);
%var=energy/S2N;
noise=sqrt(var)*randn([m,n]);

y=y+noise;
figure(3)
grayimage(y);
title('Blurred image + noise');
pause

% Calculate Wiener filter % Periodogram
m=64;
n=64;
%Pxx=(abs(fft2(x)).^2)./(m*n);
% NB here we are going to use the true values of Pxx and Pdd -- can't
% of course do this in real life :(
Pxx=(spectrum2(x,m));
Pdd=spectrum2(noise,m);
H=abs(fft2(h,m,n));
% this pads h with zeros up to size nxm before transforming
%Pdd=var;
%G=(conj(H))./((abs(H).^2)+.001*max(max(abs(H))));
G=(conj(H).*Pxx)./((abs(H).^2).*Pxx+Pdd);
g=real(ifft2(G));
g=fftshift(g);
gwin=g(m/2+1-15:m/2+1+15,n/2+1-15:n/2+1+15);

% Set up rotated Hamming window
U=15;
w2=zeros(2*U+1,2*U+1);
for n1=1:2*U+1
	for n2=1:2*U+1
		u=sqrt((n1-U)^2+(n2-U)^2);
		if u<=U
			w2(n1,n2)=0.54+0.46*cos(pi*u/U);
		end
	end
end
gwin=gwin.*w2;

xest=filter2(gwin,y,'valid');
figure(4)
grayimage(xest);
title('Restored image');

