clear
%Given signal-to-noise ratio (SNR)
rho = 17;

%Read image data and convert it to double
g = imread('\siva\course\cs445\spring_02\assign3\peppers_degraded.tif');
g = im2double(g);
[height,width] = size(g);

figure(1)
imshow(g);

%DFT of image data
G = fft2(g);

%Define the Butterworth filter H(u,v)


r0 = 36;

u = -(height/2):1:((height/2)-1);
u = (u.*u)';
v = -(width/2):1:((width/2)-1);
v = (v.*v)';
eu = ones(height,1);
ev = ones(width,1);
H = (u*eu' + ev*v')/(r0^2);
H = 1 + H.^5;
H1 = ones(size(H))./H;

%%%%% The following can also be used for creating H matrix %%%%%
%%%%% but it uses for loops and takes considerably longer %%%%%%%
%for u=1:height
%   for v=1:width
%      H(u,v) = 1/(1+ ((((u-1-height/2)^2+(v-1-width/2)^2)/(r0^2))^5));
%   end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Defined filter has origin of frequency axes at the matrix center
%Use fftshift to get the origin to the top-right corner of matrix
H = fftshift(H);

%Wiener restoration filter 
R = conj(H)./(H.*conj(H) + 1/rho);

%Restoration
Fhat = R.*G;

%Spatial domain output
fhat = ifft2(Fhat);

%This must be real but there can be small imaginary 
%component due to round-off
fhat = real(fhat);

figure(2)
imshow(fhat,[])

%Compare with ideal image
f = imread('\siva\course\cs445\images\peppers_256.tif');
f = im2double(f);

error = f - fhat;
error = error.*error;
sum(error(:))/prod(size(f))
