% EE4830 Digital Image Processing
% Class Demo: MATLAB basics
%				  and image toolbox
% Feb 21, 2001
% <xlx@ctr.columbia.edu>

%read image Lena
lena=imread('lena.jpg');
%display
clf;
imshow(lena);
disp('press any key to continue ... ');
pause

% ---------- Section 1: basic matrix operation ---------------------
% use image as our sample matrix
% show simple matrix manipulation: transpose, concatenation, sort, flip, ... 
% and common operators: ', :, [], ...

lena=double(imread('lena.png'));
% Lena  versions
% transposed, flipped, and rotated
colormap(gray(256));
%clf RESET
image( [ lena, fliplr( lena ) ; rot90(lena'), flipud(fliplr(lena))] );
axis off;
pause

% Still Lena ?                                                                          
% negative, crop and interpolation, quant, etc.                                                                   
croplena=imresize(lena(65:192,65:192),2);
step=32;
image( [ lena, ( 255 - lena ); croplena, step*floor(lena/step)+step/2] );
axis off;
pause

% ----------- Section 2: FFT --------------------------------------
% intuition on image transform 
% example: natrual image, i.e. lena
% 			  and artificial image: stripes

flena = fft2(lena);
%show the amplitude and phase
subplot(221); image(lena);
subplot(222); imagesc(fftshift(log10(abs(flena))));
subplot(223); imagesc(angle(flena));
subplot(224); mesh(fftshift(log10(abs(flena)))); axis tight;
% Here the magnitude is log-ed so that high-freqquency components will properly show up 
pause
%zoom in to see a color version
figure(2); mesh(fftshift(log10(abs(flena))));
pause
%build artificial stripes
close(2); clf;
%horizontal stripe
hstripe = repmat([ones(32,256);zeros(32,256)], 4, 1);
subplot(221); imagesc(hstripe);
%diagonal stripe
dstripe=zeros(256);
onerow = [ones(1,32),zeros(1,32)];
onerow = repmat(onerow, 1, 4);
for i=1:256
   dstripe(i,:)= onerow;
   onerow = [onerow(2:256),onerow(1)];
end
subplot(223); imagesc(dstripe);
pause

%Do FFT
%log and add some number for proper display
subplot(222); imagesc(fftshift(log10(log10(abs(fft2(hstripe))+11))));
subplot(224); imagesc(fftshift(log10(log10(abs(fft2(dstripe))+11))));
pause

%Explanation
%miltiple 1-D step function
clf;
step=[zeros(1,64), ones(1,128),zeros(1,64)];
subplot(221); plot(step); axis tight;
title('1-D step function');
subplot(222); plot((fftshift(abs(fft(step))))); axis tight;
title('sinc');
pause
%zoom in
axis([60,100,0,5]);
title('sinc (zoom)');
pause
%2-D case: seperable transform
subplot(223); imagesc(hstripe);
subplot(224); imagesc(fftshift(log10(log10(abs(fft2(hstripe))+11))));
pause
%another explaination
subplot(223); imagesc(dstripe);
subplot(224); imagesc(fftshift(log10(log10(abs(fft2(dstripe))+11))));
pause

%Take a non-trivial example
figure(2);
lowlib=imread('lowlib.jpg');
image(lowlib); axis off;
pause
%Architects often talk about the ballance between horizontal and vertical lines
%Now see how this will show up in the frequency domain
lowlib=imread('lowlib_crop.jpg');
clf reset;
subplot(121); image(lowlib); axis off;
pause
%take the luminance channel
lowlib=double(lowlib(:,:,1));
%Do FFT and display
subplot(122); imagesc(fftshift(log10(abs(fft2(lowlib))+11))); axis off;

%End of Demo
pause
close all;