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

% original girl2 image (aka tiffany)
X = ReadBin('girl2.bin',256);

% with broadband noise
XN = ReadBin('girl2Noise32.bin',256);

% with hi pass noise
XNhi = ReadBin('girl2Noise32Hi.bin',256);


%----------------------------------------------------------------------
%
% P2(a): Display images and compute the MSE of each noisy image
%

xx = (XNhi - X).^2;
MSE_Nhi = mean(xx(:));
disp(['MSE girl2Noise32Hi.bin:  ' num2str(MSE_Nhi)])
xx = (XN - X).^2;
MSE_N = mean(xx(:));
disp(['MSE girl2Noise32.bin:    ' num2str(MSE_N)])
disp(' ');

figure(1);image(X);colormap(gray(256));axis('image','off');
title('Original Tiffany Image','FontSize',18);
print -deps Mgirl2.eps;
figure(2);image(XNhi);colormap(gray(256));axis('image','off');
title('girl2Noise32Hi','FontSize',18);
print -deps Mgirl2Noise32Hi.eps;
figure(3);image(XN);colormap(gray(256));axis('image','off');
title('girl2Noise32','FontSize',18);
print -deps Mgirl2Noise32.eps;

%----------------------------------------------------------------------
%
% P2(b): Apply isotropic ideal LPF with U_cutoff = 64.
%   Use circular convolution for this.
%

U_cutoff = 64;
[U,V] = meshgrid(-128:127,-128:127);
HLtildeCenter = double(sqrt(U.^2 + V.^2) <= U_cutoff);
HLtilde = fftshift(HLtildeCenter);

% apply to original girl2 image and compute MSE
Y1 = ifft2(fft2(X) .* HLtilde);
yy = (Y1 - X).^2;
MSE_Y1 = mean(yy(:));
disp(['MSE:  ideal LPF on girl2:      ' num2str(MSE_Y1)])

% apply to image with hi pass noise; compute MSE and ISNR
Y1Nhi = ifft2(fft2(XNhi) .* HLtilde);
yy = (Y1Nhi - X).^2;
MSE_Y1Nhi = mean(yy(:));
disp(['MSE:  ideal LPF on Noise32Hi:  ' num2str(MSE_Y1Nhi)])
ISNR_Y1Nhi = 10*log10(MSE_Nhi/MSE_Y1Nhi);
disp(['ISNR: ideal LPF on Noise32Hi:  ' num2str(ISNR_Y1Nhi) ' dB'])

% apply to image with broadband noise; compute MSE and ISNR
Y1N = ifft2(fft2(XN) .* HLtilde);
yy = (Y1N - X).^2;
MSE_Y1N = mean(yy(:));
disp(['MSE:  ideal LPF on Noise32:    ' num2str(MSE_Y1N)])
ISNR_Y1N = 10*log10(MSE_N/MSE_Y1N);
disp(['ISNR: ideal LPF on Noise32:    ' num2str(ISNR_Y1N) ' dB'])
disp(' ');

figure(4);image(stretch(Y1));colormap(gray(256));axis('image','off');
title('LPF on girl2','FontSize',18);
print -deps MY1.eps;
figure(5);image(stretch(Y1Nhi));colormap(gray(256));axis('image','off');
title('LPF on Noise32Hi','FontSize',18);
print -deps MY1Nhi.eps;
figure(6);image(stretch(Y1N));colormap(gray(256));axis('image','off');
title('LPF on Noise32','FontSize',18);
print -deps MY1N.eps;

%----------------------------------------------------------------------
%
% P2(c): Use the method of Example 4 p. 5.65 to apply Gaussian
% LPF with U_cutoff = 64.
%

U_cutoff_G = 64;
SigmaG = 0.19 * 256 / U_cutoff_G;
[U,V] = meshgrid(-128:127,-128:127);
GtildeCenter = exp((-2*pi^2*SigmaG^2)/(256.^2)*(U.^2 + V.^2));
Gtilde = fftshift(GtildeCenter);
G  = ifft2(Gtilde);
G2 = fftshift(G);
ZPG2 = zeros(512,512);
ZPG2(1:256,1:256) = G2;

% apply to original girl2 image and compute MSE
ZPX = zeros(512,512);
ZPX(1:256,1:256) = X;
yy = ifft2(fft2(ZPX).*fft2(ZPG2));
Y2 = yy(129:384,129:384);
yy = (Y2 - X).^2;
MSE_Y2 = mean(yy(:));
disp(['MSE:  Gaussian LPF on girl2:     ' num2str(MSE_Y2)])

% apply to image with hi pass noise; compute MSE and ISNR
ZPX = zeros(512,512);
ZPX(1:256,1:256) = XNhi;
yy = ifft2(fft2(ZPX).*fft2(ZPG2));
Y2Nhi = yy(129:384,129:384);
yy = (Y2Nhi - X).^2;
MSE_Y2Nhi = mean(yy(:));
disp(['MSE:  Gaussian LPF on Noise32Hi: ' num2str(MSE_Y2Nhi)])
ISNR_Y2Nhi = 10*log10(MSE_Nhi/MSE_Y2Nhi);
disp(['ISNR: Gaussian LPF on Noise32Hi: ' num2str(ISNR_Y2Nhi) ' dB'])

% apply to image with broadband noise; compute MSE and ISNR
ZPX = zeros(512,512);
ZPX(1:256,1:256) = XN;
yy = ifft2(fft2(ZPX).*fft2(ZPG2));
Y2N = yy(129:384,129:384);
yy = (Y2N - X).^2;
MSE_Y2N = mean(yy(:));
disp(['MSE:  Gaussian LPF on Noise32:   ' num2str(MSE_Y2N)])
ISNR_Y2N = 10*log10(MSE_N/MSE_Y2N);
disp(['ISNR: Gaussian LPF on Noise32:   ' num2str(ISNR_Y2N) ' dB'])
disp(' ');

figure(7);image(stretch(Y2));colormap(gray(256));axis('image','off');
title('Gauss1 on girl2','FontSize',18);
print -deps MY2.eps;
figure(8);image(stretch(Y2Nhi));colormap(gray(256));axis('image','off');
title('Gauss1 on Noise32Hi','FontSize',18);
print -deps MY2Nhi.eps;
figure(9);image(stretch(Y2N));colormap(gray(256));axis('image','off');
title('Gauss1 on Noise32','FontSize',18);
print -deps MY2N.eps;

%----------------------------------------------------------------------
%
% P2(d): Use the method of Example 4 p. 5.65 to apply Gaussian
% LPF with U_cutoff = 34.
%

U_cutoff_G = 77.5;
SigmaG = 0.19 * 256 / U_cutoff_G;
GtildeCenter = exp((-2*pi^2*SigmaG^2)/(256.^2)*(U.^2 + V.^2));
Gtilde = fftshift(GtildeCenter);
G  = ifft2(Gtilde);
G2 = fftshift(G);
ZPG2 = zeros(512,512);
ZPG2(1:256,1:256) = G2;

% apply to original girl2 image and compute MSE
ZPX = zeros(512,512);
ZPX(1:256,1:256) = X;
yy = ifft2(fft2(ZPX).*fft2(ZPG2));
Y3 = yy(129:384,129:384);
yy = (Y3 - X).^2;
MSE_Y3 = mean(yy(:));
disp(['MSE:  Gaussian2 LPF on girl2:     ' num2str(MSE_Y3)])

% apply to image with hi pass noise; compute MSE and ISNR
ZPX = zeros(512,512);
ZPX(1:256,1:256) = XNhi;
yy = ifft2(fft2(ZPX).*fft2(ZPG2));
Y3Nhi = yy(129:384,129:384);
yy = (Y3Nhi - X).^2;
MSE_Y3Nhi = mean(yy(:));
disp(['MSE:  Gaussian2 LPF on Noise32Hi: ' num2str(MSE_Y3Nhi)])
ISNR_Y3Nhi = 10*log10(MSE_Nhi/MSE_Y3Nhi);
disp(['ISNR: Gaussian2 LPF on Noise32Hi: ' num2str(ISNR_Y3Nhi) ' dB'])

% apply to image with broadband noise; compute MSE and ISNR
ZPX = zeros(512,512);
ZPX(1:256,1:256) = XN;
yy = ifft2(fft2(ZPX).*fft2(ZPG2));
Y3N = yy(129:384,129:384);
yy = (Y3N - X).^2;
MSE_Y3N = mean(yy(:));
disp(['MSE:  Gaussian2 LPF on Noise32:   ' num2str(MSE_Y3N)])
ISNR_Y3N = 10*log10(MSE_N/MSE_Y3N);
disp(['ISNR: Gaussian2 LPF on Noise32:   ' num2str(ISNR_Y3N) ' dB'])

figure(10);image(stretch(Y3));colormap(gray(256));axis('image','off');
title('Gauss2 on girl2','FontSize',18);
print -deps MY3.eps;
figure(11);image(stretch(Y3Nhi));colormap(gray(256));axis('image','off');
title('Gauss2 on Noise32Hi','FontSize',18);
print -deps MY3Nhi.eps;
figure(12);image(stretch(Y3N));colormap(gray(256));axis('image','off');
title('Gauss2 on Noise32','FontSize',18);
print -deps MY3N.eps;

















