%
% Make images for HW5 P2
%

X = ReadBin('girl2.bin',256);
figure(1);image(X);colormap(gray(256));axis('image','off');
title('Original Image','FontSize',18);

Uc = 64;
SigmaN = 32;

%
% Ideal low pass and hi pass filters
%
[U,V] = meshgrid(-128:127,-128:127);
HLtildeCenter = double(sqrt(U.^2 + V.^2) <= Uc);
HLtilde = fftshift(HLtildeCenter);

HHtildeCenter = double(sqrt(U.^2 + V.^2) > Uc);
HHtilde = fftshift(HHtildeCenter);

%
% Noise
N1 = random('norm',0,SigmaN,256,256);
Nhi = ifft2(HHtilde .* fft2(N1));
Nhi = SigmaN * Nhi / std(Nhi(:));

%
% Make noisy images
%
girl2Noise32   = X + N1;
girl2Noise32Hi = X + Nhi;

%
% Write the noisy images
fid = fopen('girl2Noise32.bin','w+');
fwrite(fid,stretch(girl2Noise32'),'uchar');
fclose(fid);
fid = fopen('girl2Noise32Hi.bin','w+');
fwrite(fid,stretch(girl2Noise32Hi'),'uchar');
fclose(fid);

%
% Apply idea LPF
%
Ypure = ifft2(fft2(X) .* HLtilde);
Y32   = ifft2(fft2(girl2Noise32) .* HLtilde);
Y32Hi = ifft2(fft2(girl2Noise32Hi) .* HLtilde);

%
% Gaussian filter
%
Ucg = 64;
%Ucg = 34;
SigmaG = 0.19 * 256 / Ucg;
GtildeCenter = exp((-2*pi*SigmaG^2)/(256.^2)*(U.^2 + V.^2));
Gtilde = fftshift(GtildeCenter);

% this is circular
YGpure = ifft2(fft2(X) .* Gtilde);
YG32   = ifft2(fft2(girl2Noise32) .* Gtilde);
YG32Hi = ifft2(fft2(girl2Noise32Hi) .* Gtilde);

%
% Use Example 4 p. 5.65
G2 = ifft2(Gtilde);
G2 = fftshift(G2);
ZPG2 = zeros(512,512);
ZPG2(1:256,1:256) = G2;
ZPX = zeros(512,512);
ZPX(1:256,1:256) = X;
xx = ifft2(fft2(ZPX).*fft2(ZPG2));
YG2pure = xx(129:384,129:384);

ZPX32 = zeros(512,512);
ZPX32(1:256,1:256) = girl2Noise32;
xx = ifft2(fft2(ZPX32).*fft2(ZPG2));
YG2_32 = xx(129:384,129:384);

ZPX32Hi = zeros(512,512);
ZPX32Hi(1:256,1:256) = girl2Noise32Hi;
xx = ifft2(fft2(ZPX32Hi).*fft2(ZPG2));
YG2_32Hi = xx(129:384,129:384);







xx = (girl2Noise32 - X).^2;
MSEg32 = mean(xx(:));
disp(['MSE g32     ' num2str(MSEg32)])
xx = (girl2Noise32Hi - X).^2;
MSEg32Hi = mean(xx(:));
disp(['MSE g32Hi   ' num2str(MSEg32Hi)])

xx = (Ypure - X).^2;
MSEYpure = mean(xx(:));
disp(['MSE Ypure   ' num2str(MSEYpure)])
xx = (YGpure - X).^2;
MSEYGpure = mean(xx(:));
disp(['MSE YGpure  ' num2str(MSEYGpure)])
xx = (YG2pure - X).^2;
MSEYG2pure = mean(xx(:));
disp(['MSE YG2pure ' num2str(MSEYG2pure)])

xx = (Y32 - X).^2;
MSEY32 = mean(xx(:));
disp(['MSE Y32     ' num2str(MSEY32)])
xx = (YG32 - X).^2;
MSEYG32 = mean(xx(:));
disp(['MSE YG32    ' num2str(MSEYG32)])
xx = (YG2_32 - X).^2;
MSEYG2_32 = mean(xx(:));
disp(['MSE YG2 32  ' num2str(MSEYG2_32)])

xx = (Y32Hi - X).^2;
MSEY32Hi = mean(xx(:));
disp(['MSE Y32Hi   ' num2str(MSEY32Hi)])
xx = (YG32Hi - X).^2;
MSEYG32Hi = mean(xx(:));
disp(['MSE YG32Hi  ' num2str(MSEYG32Hi)])
xx = (YG2_32Hi - X).^2;
MSEYG2_32Hi = mean(xx(:));
disp(['MSE YG2 32Hi' num2str(MSEYG2_32Hi)])


figure(2);image(stretch(HLtildeCenter));colormap(gray(256));axis('image','off');
title('HLtildeCenter','FontSize',18);
figure(3);image(stretch(HHtildeCenter));colormap(gray(256));axis('image','off');
title('HHtildeCenter','FontSize',18);
figure(4);image(stretch(N1));colormap(gray(256));axis('image','off');
title('N1','FontSize',18);
figure(5);image(stretch(Nhi));colormap(gray(256));axis('image','off');
title('Nhi','FontSize',18);
figure(6);image(stretch(girl2Noise32));colormap(gray(256));axis('image','off');
title('girl2Noise32','FontSize',18);
figure(7);image(stretch(girl2Noise32Hi));colormap(gray(256));axis('image','off');
title('girl2Noise32Hi','FontSize',18);

figure(8);image(stretch(Y32));colormap(gray(256));axis('image','off');
title('Y32','FontSize',18);
figure(9);image(stretch(Y32Hi));colormap(gray(256));axis('image','off');
title('Y32Hi','FontSize',18);
figure(10);image(stretch(Ypure));colormap(gray(256));axis('image','off');
title('Ypure','FontSize',18);

figure(11);image(stretch(GtildeCenter));colormap(gray(256));axis('image','off');
title('GtildeCenter','FontSize',18);
figure(12);image(stretch(YGpure));colormap(gray(256));axis('image','off');
title('YGpure','FontSize',18);
figure(13);image(stretch(YG32));colormap(gray(256));axis('image','off');
title('YG32','FontSize',18);
figure(14);image(stretch(YG32Hi));colormap(gray(256));axis('image','off');
title('YG32Hi','FontSize',18);

figure(15);image(stretch(YG2pure));colormap(gray(256));axis('image','off');
title('YG2pure','FontSize',18);
figure(16);image(stretch(YG2_32));colormap(gray(256));axis('image','off');
title('YG2\_32','FontSize',18);
figure(17);image(stretch(YG2_32Hi));colormap(gray(256));axis('image','off');
title('YG2\_32Hi','FontSize',18);




%figure(5);
%image(stretch(log(1 + abs(fftshift(fft2(X2))))));
%colormap(gray(256));axis('image','off');
%title('X2 FFT','FontSize',18);

