%
% hw04.m
%
% 04/05/04 jph
%  03/27/2013 the LTI filtering problem is now postponed to a later homework
%  02/15/2015 added problems 6 and 7. jph.
%

diary hw04.diary

%---------------------------------------------------------------------------
% Problem 1
%
% Here the images are of size 8x8.  For creating these images, it is
% important that the columns and rows use zero-based indexing.
%
% Make an 8x8 matrix COLS where each entry gives the zero-based index of the
% column and an 8x8 matrix ROWS where each entry gives the zero-based index
% of the row.
%
[COLS,ROWS] = meshgrid(0:7,0:7);  

%
% Make the image I1.  Display the real and imaginary parts as 8 bpp images
% with full-scale contrast.
%
u0 = 2.0;
v0 = 2.0;
I1 = 0.5 * exp(j*2*pi/8*(u0*COLS + v0*ROWS));

%
% Perform full-scale contrast stretch on real and imaginary parts of I1 and
%   display as images.
%
I1R = stretch(real(I1));
I1I = stretch(imag(I1));
figure(1);
image(I1R);colormap(gray(256));axis('image','off');
title('Re[{\fontname{times}I}_1]','FontSize',18);
print -deps MI1R.eps;
print -dpdf MI1R.pdf;

figure(2);
image(I1I);colormap(gray(256));axis('image','off');
title('Im[{\fontname{times}I}_1]','FontSize',18);
print -deps MI1I.eps;
print -dpdf MI1I.pdf;

%
% Compute the centered 2D DFT of I1.  Print out real and imaginary parts
% as ascii arrays.
%
Itilde1 = fftshift(fft2(I1));
fprintf(1,'%s\n','Re[DFT(I1)]:');
disp(round(real(Itilde1) * 10^4)*10^(-4));
fprintf(1,'%s\n','Im[DFT(I1)]:');
disp(round(imag(Itilde1) * 10^4)*10^(-4));

%---------------------------------------------------------------------------
% Problem 2
%

%
% Make the image I2.
%
I2 = 0.5 * exp(-j*2*pi/8*(u0*COLS + v0*ROWS));

%
% Display real and imaginary parts as 8 bpp images with full scale contrast.
%
I2R = stretch(real(I2));
I2I = stretch(imag(I2));

figure(3);
image(I2R);colormap(gray(256));axis('image','off');
title('Re[{\fontname{times}I}_2]','FontSize',18);
print -deps MI2R.eps;
print -dpdf MI2R.pdf;

figure(4);
image(I2I);colormap(gray(256));axis('image','off');
title('Im[{\fontname{times}I}_2]','FontSize',18);
print -deps MI2I.eps;
print -dpdf MI2I.pdf;

%
% Compute the centered 2D DFT of I2.  Print out real and imaginary parts
% as ascii arrays.
%
Itilde2 = fftshift(fft2(I2));
fprintf(1,'%s\n','Re[DFT(I2)]:');
disp(round(real(Itilde2) * 10^4)*10^(-4));
fprintf(1,'%s\n','Im[DFT(I2)]:');
disp(round(imag(Itilde2) * 10^4)*10^(-4));

%---------------------------------------------------------------------------
% Problem 3
%

%
% Make the image I3.
%
I3 = cos(2*pi/8*(u0*COLS + v0*ROWS));

%
% Display as 8 bpp image with full scale contrast.
%
figure(5);
image(stretch(I3));colormap(gray(256));axis('image','off');
title('{\fontname{times}I}_3','FontSize',18);
print -deps MI3R.eps;
print -dpdf MI3R.pdf;

%
% Compute the centered 2D DFT of I3.  Print out real and imaginary parts
% as ascii arrays.
%
Itilde3 = fftshift(fft2(I3));
fprintf(1,'%s\n','Re[DFT(I3)]:');
disp(round(real(Itilde3) * 10^4)*10^(-4));
fprintf(1,'%s\n','Im[DFT(I3)]:');
disp(round(imag(Itilde3) * 10^4)*10^(-4));

%---------------------------------------------------------------------------
% Problem 4
%

%
% Make the image I4.
%
I4 = sin(2*pi/8*(u0*COLS + v0*ROWS));

%
% Display as 8 bpp image with full scale contrast.
%
figure(6);
image(stretch(I4));colormap(gray(256));axis('image','off');
title('{\fontname{times}I}_4','FontSize',18);
print -deps MI4R.eps;
print -dpdf MI4R.pdf;

%
% Compute the centered 2D DFT of I4.  Print out real and imaginary parts
% as ascii arrays.
%
Itilde4 = fftshift(fft2(I4));
fprintf(1,'%s\n','Re[DFT(I4)]:');
disp(round(real(Itilde4) * 10^4)*10^(-4));
fprintf(1,'%s\n','Im[DFT(I4)]:');
disp(round(imag(Itilde4) * 10^4)*10^(-4));

%---------------------------------------------------------------------------
% Problem 5
%

%
% Make the image I5.
%
u1 = 1.5;
v1 = 1.5;
I5 = cos(2*pi/8*(u1*COLS + v1*ROWS));

%
% Display as 8 bpp image with full scale contrast.
%
figure(7);
image(stretch(I5));colormap(gray(256));axis('image','off');
title('{\fontname{times}I}_5','FontSize',18);
print -deps MI5R.eps;
print -dpdf MI5R.pdf;

%
% Compute the centered 2D DFT of I5.  Print out real and imaginary parts
% as ascii arrays.
%
Itilde5 = fftshift(fft2(I5));
fprintf(1,'%s\n','Re[DFT(I5)]:');
disp(round(real(Itilde5) * 10^4)*10^(-4));
fprintf(1,'%s\n','Im[DFT(I5)]:');
disp(round(imag(Itilde5) * 10^4)*10^(-4));

diary off

%---------------------------------------------------------------------------
% Problem 6
%

%
% Cameraman
%
X = ReadBin('camera.bin',256);
Xtilde = fftshift(fft2(X));
XtildeR = stretch(real(Xtilde));
XtildeI = stretch(imag(Xtilde));
XtildeMag = stretch(log(1 + abs(Xtilde)));
XtildePhase = stretch(angle(Xtilde));

figure(8);
image(X);colormap(gray(256));axis('image','off');
title('{\fontname{times}camera}','FontSize',18);
print -deps Mcamera.eps;
print -dpdf Mcamera.pdf;

figure(9);
image(XtildeR);colormap(gray(256));axis('image','off');
title('{\fontname{times}Re[DFT(camera)]}','FontSize',18);
print -deps McameraFFTR.eps;
print -dpdf McameraFFTR.pdf;

figure(10);
image(XtildeI);colormap(gray(256));axis('image','off');
title('{\fontname{times}Im[DFT(camera)]}','FontSize',18);
print -deps McameraFFTI.eps;
print -dpdf McameraFFTI.pdf;

figure(11);
image(XtildeMag);colormap(gray(256));axis('image','off');
title('{\fontname{times}camera log magnitude spectrum}','FontSize',18);
print -deps McameraFFTMag.eps;
print -dpdf McameraFFTMag.pdf;

figure(12);
image(XtildePhase);colormap(gray(256));axis('image','off');
title('{\fontname{times}arg[DFT(camera)]}','FontSize',18);
print -deps McameraFFTPhase.eps;
print -dpdf McameraFFTPhase.pdf;

%
% Salesman
%
X = ReadBin('salesman.bin',256);
Xtilde = fftshift(fft2(X));
XtildeR = stretch(real(Xtilde));
XtildeI = stretch(imag(Xtilde));
XtildeMag = stretch(log(1 + abs(Xtilde)));
XtildePhase = stretch(angle(Xtilde));

figure(13);
image(X);colormap(gray(256));axis('image','off');
title('{\fontname{times}salesman}','FontSize',18);
print -deps Msalesman.eps;
print -dpdf Msalesman.pdf;

figure(14);
image(XtildeR);colormap(gray(256));axis('image','off');
title('{\fontname{times}Re[DFT(salesman)]}','FontSize',18);
print -deps MsalesmanFFTR.eps;
print -dpdf MsalesmanFFTR.pdf;

figure(15);
image(XtildeI);colormap(gray(256));axis('image','off');
title('{\fontname{times}Im[DFT(salesman)]}','FontSize',18);
print -deps MsalesmanFFTI.eps;
print -dpdf MsalesmanFFTI.pdf;

figure(16);
image(XtildeMag);colormap(gray(256));axis('image','off');
title('{\fontname{times}salesman log magnitude spectrum}','FontSize',18);
print -deps MsalesmanFFTMag.eps;
print -dpdf MsalesmanFFTMag.pdf;

figure(17);
image(XtildePhase);colormap(gray(256));axis('image','off');
title('{\fontname{times}arg[DFT(salesman)]}','FontSize',18);
print -deps MsalesmanFFTPhase.eps;
print -dpdf MsalesmanFFTPhase.pdf;

%
% Head
%
X = ReadBin('head.bin',256);
Xtilde = fftshift(fft2(X));
XtildeR = stretch(real(Xtilde));
XtildeI = stretch(imag(Xtilde));
XtildeMag = stretch(log(1 + abs(Xtilde)));
XtildePhase = stretch(angle(Xtilde));

figure(18);
image(X);colormap(gray(256));axis('image','off');
title('{\fontname{times}head}','FontSize',18);
print -deps Mhead.eps;
print -dpdf Mhead.pdf;

figure(19);
image(XtildeR);colormap(gray(256));axis('image','off');
title('{\fontname{times}Re[DFT(head)]}','FontSize',18);
print -deps MheadFFTR.eps;
print -dpdf MheadFFTR.pdf;

figure(20);
image(XtildeI);colormap(gray(256));axis('image','off');
title('{\fontname{times}Im[DFT(head)]}','FontSize',18);
print -deps MheadFFTI.eps;
print -dpdf MheadFFTI.pdf;

figure(21);
image(XtildeMag);colormap(gray(256));axis('image','off');
title('{\fontname{times}head log magnitude spectrum}','FontSize',18);
print -deps MheadFFTMag.eps;
print -dpdf MheadFFTMag.pdf;

figure(22);
image(XtildePhase);colormap(gray(256));axis('image','off');
title('{\fontname{times}arg[DFT(head)]}','FontSize',18);
print -deps MheadFFTPhase.eps;
print -dpdf MheadFFTPhase.pdf;

%
% eyeR
%
X = ReadBin('eyeR.bin',256);
Xtilde = fftshift(fft2(X));
XtildeR = stretch(real(Xtilde));
XtildeI = stretch(imag(Xtilde));
XtildeMag = stretch(log(1 + abs(Xtilde)));
XtildePhase = stretch(angle(Xtilde));

figure(23);
image(X);colormap(gray(256));axis('image','off');
title('{\fontname{times}eyeR}','FontSize',18);
print -deps MeyeR.eps;
print -dpdf MeyeR.pdf;

figure(24);
image(XtildeR);colormap(gray(256));axis('image','off');
title('{\fontname{times}Re[DFT(eyeR)]}','FontSize',18);
print -deps MeyeRFFTR.eps;
print -dpdf MeyeRFFTR.pdf;

figure(25);
image(XtildeI);colormap(gray(256));axis('image','off');
title('{\fontname{times}Im[DFT(eyeR)]}','FontSize',18);
print -deps MeyeRFFTI.eps;
print -dpdf MeyeRFFTI.pdf;

figure(26);
image(XtildeMag);colormap(gray(256));axis('image','off');
title('{\fontname{times}eyeR log magnitude spectrum}','FontSize',18);
print -deps MeyeRFFTMag.eps;
print -dpdf MeyeRFFTMag.pdf;

figure(27);
image(XtildePhase);colormap(gray(256));axis('image','off');
title('{\fontname{times}arg[DFT(eyeR)]}','FontSize',18);
print -deps MeyeRFFTPhase.eps;
print -dpdf MeyeRFFTPhase.pdf;

%---------------------------------------------------------------------------
% Problem 7
%

X = ReadBin('camera.bin',256);
Xtilde = fft2(X);
%
% since cos(0)=1 and sin(0)=0, Re[J1tilde] = |Xtilde|
%    and Im[J1tilde] = 0.
%
J1tilde = abs(Xtilde);
J1 = real(ifft2(J1tilde));
J1prime = stretch(log(1 + J1));
figure(28);
image(J1prime);colormap(gray(256));axis('image','off');
title('{\fontname{times}J1''}','FontSize',18);
print -deps MJ1prime.eps;
print -dpdf MJ1prime.pdf;

J2tilde = exp(j*angle(Xtilde));
J2 = stretch(real(ifft2(J2tilde)));
figure(29);
image(J2);colormap(gray(256));axis('image','off');
title('{\fontname{times}J2}','FontSize',18);
print -deps MJ2.eps;
print -dpdf MJ2.pdf;


