clear
%Size of "window" around the interfering frequency that 
%need to be filtered out 
N=7;  

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

figure(1)
imshow(g,[]);

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

%Shift origin of freq. axes to matrix center for display purposes
G = fftshift(G);

%Display magnitude spectrum in log-scale
figure(2)
imshow(log(1+abs(G)),[]);

%Interactively locate the frequency (location in frequency domain) of interference
pixval

keyboard

%Initialize restoration filter to ones and set some values to 
%zero, depending on the location of interference
R = ones(height,width);

%From the spectrum, the interference is centered around (39,59) and (219,199)
%these values, denoted by variables (p1,q1) and (p2,q2) are to be entered 
%in the previous steps (keyboard command). 
R(p1-N:p1+N,q1-N:q1+N) = 0;
R(p2-N:p2+N,q2-N:q2+N) = 0;

%The above filtering was somewhat inadequate
%So decided to "zero-out" some some DFT values
R(p1,:) = 0;
R(p2,:) = 0;
R(:,q1) = 0;
R(:,q2) = 0;

%Check to see if the designed filter is what was intended
figure(3)
imshow(R,[])

%Filtering
Fhat = R .* G;

%Move origin of frequency axes back to top-right corner befor inverse DFT
Fhat = fftshift(Fhat);

%Spatial domain output
fhat = ifft2(Fhat);
fhat = real(fhat);

figure(4)
imshow(fhat,[]);

%Compare with ideal image
f = imread('\matlab\toolbox\images\imdemos\ic.tif');
f = im2double(f);

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