%
% P1.m
%
% 04/05/2015 jph
%


%----------------------------------------------------------------------
%
% P1(a): apply a 7x7 average filter to salesman image by doing linear 
% convolution in the image domain.
%
% Handle edge effects by zero padding.
%

X = ReadBin('salesman.bin',256);
figure(1);image(X);colormap(gray(256));axis('image','off');
title('Original Image','FontSize',18);
print -deps Msalesman.eps;

X2 = zeros(262,262);
X2(4:259,4:259) = X;
Y2 = zeros(262,262);

for row=4:259
  for col=4:259
    Y2(row,col) = sum(sum(X2(row-3:row+3,col-3:col+3)))/49;
  end
end

Y = stretch(Y2(4:259,4:259));
figure(2);image(Y);colormap(gray(256));axis('image','off');
title('Filtered Image','FontSize',18);
print -deps MY1a.eps;

% Save the result image for comparison with later results
Y1a = Y;


%----------------------------------------------------------------------
%
% P1(b): Use the method of Example 3 on page 5.61 of the Notes to do
% the same linear convolution by pointwise multiplication of DFT's.  
%
% The impulse response image H will be 128x128.
%
% The original image is already in the Matlab array X.
%

% Make the 128x128 impulse response image
H = zeros(128,128);
H(62:68,62:68) = 1/49;

% Zero pad the original image and the H image
Padsize = 256 + 128 - 1;
ZPX = zeros(Padsize,Padsize);
ZPX(1:256,1:256) = X;
figure(3);image(ZPX);colormap(gray(256));axis('image','off');
title('Zero Padded Image','FontSize',18);
print -deps MZPX1b.eps;

ZPH = zeros(Padsize,Padsize);
ZPH(1:128,1:128) = H;
figure(4);image(stretch(ZPH));colormap(gray(256));axis('image','off');
title('Zero Padded Impulse Resp','FontSize',18);
print -deps MZPH1b.eps;

% Compute DFT's of zero padded images
ZPXtilde = fft2(ZPX);
ZPHtilde = fft2(ZPH);

% Show centered log-magnitude spectra
ZPXtildeDisplay = stretch(log(1 + abs(fftshift(ZPXtilde))));
figure(5);image(ZPXtildeDisplay);colormap(gray(256));axis('image','off');
title('Log-magnitude spectrum of zero padded image','FontSize',12);
print -deps MZPXtilde1b.eps;
ZPHtildeDisplay = stretch(log(1 + abs(fftshift(ZPHtilde))));
figure(6);image(ZPHtildeDisplay);colormap(gray(256));axis('image','off');
title('Log-magnitude spectrum of zero padded H image','FontSize',12);
print -deps MZPHtilde1b.eps;

% Compute the convolution by pointwise multiplication of DFT's.
% Show the resulting zero padded image and it's centered log-magnitude
% spectrum.
ZPYtilde = ZPXtilde .* ZPHtilde;
ZPY = ifft2(ZPYtilde);

ZPYtildeDisplay = stretch(log(1 + abs(fftshift(ZPYtilde))));
figure(7);image(ZPYtildeDisplay);colormap(gray(256));axis('image','off');
title('Log-magnitude spectrum of zero padded result','FontSize',12);
print -deps MZPYtilde1b.eps;
figure(8);image(stretch(ZPY));colormap(gray(256));axis('image','off');
title('Zero Padded Result','FontSize',18);
print -deps MZPY1b.eps;

% Extract the final result image and display
Y = stretch(ZPY(65:320,65:320));
figure(9);image(Y);colormap(gray(256));axis('image','off');
title('Final Filtered Image','FontSize',18);
print -deps MY1b.eps;

% Compare this result image with the one from part (a)
disp(['(b): max difference from part (a): ' num2str(max(max(abs(Y-Y1a))))])


%----------------------------------------------------------------------
%
% P1(c): Use the method of Example 5 on page 5.76 of the Notes to do
% the same linear convolution again by pointwise multiplication of 
% DFT's, this time using a 256x256 true zero-phase impulse response.
%
% The impulse response image H will be 256x256.
%
% The original image is still in the Matlab array X.
%

% Make the 256x256 impulse response image H
% Put the "square" in the center
H = zeros(256,256);
H(126:132,126:132) = 1/49;

% Now use fftshift to get the true zero-phase impulse response image
H2 = fftshift(H);
figure(10);image(stretch(H2));colormap(gray(256));axis('image','off');
title('Zero Phase Impulse Resp','FontSize',18);
print -deps MH1c.eps;

% Zero pad the input image
ZPX = zeros(512,512);
ZPX(1:256,1:256) = X;

% Make the zero padded impulse response image as in Example 5
% on page 5.76 of the notes
ZPH2 = zeros(256,256);
ZPH2(1:128,1:128) = H2(1:128,1:128);
ZPH2(1:128,385:512) = H2(1:128,129:256);
ZPH2(385:512,1:128) = H2(129:256,1:128);
ZPH2(385:512,385:512) = H2(129:256,129:256);
figure(11);image(stretch(ZPH2));colormap(gray(256));axis('image','off');
title('Zero Padded zero-phase H','FontSize',18);
print -deps MZPH1c.eps;

% Compute the filtered result by pointwise multiplication of DFT's
Y = ifft2(fft2(ZPX) .* fft2(ZPH2));
Y = stretch(Y(1:256,1:256));
figure(12);image(Y);colormap(gray(256));axis('image','off');
title('Final Filtered Image','FontSize',18);
print -deps MY1c.eps;

% Compare this result image with the one from part (a)
disp(['(c): max difference from part (a): ' num2str(max(max(abs(Y-Y1a))))])





















