%%% ImagePhaseDemo.m
%%% 
%%% Programmer:    Prof. Brian L. Evans
%%% Organization:  The University of Texas at Austin
%%% Created:       February 21, 2008
%%%
%%% Requires the image processing toolbox in Matlab

aerialview = imread('westconcordorthophoto.png');

grayramp = 0 : 1/255 : 1;
graymap = [grayramp; grayramp; grayramp]';

figure(1);
imshow(aerialview);
colormap(graymap);
title('Original image');

aerialviewfft = fft2(single(aerialview));

aerialviewfftmag = abs(aerialviewfft);
aerialview2fp = ifft2(aerialviewfftmag);

%%% Show the contribution in the spatial domain of the magnitude of the FFT 
%%% of the original image
figure(2);
imagesc(aerialview2fp);
colormap(graymap);
title('Image corresponding to magnitude of FFT of the image');

%%% Show the contribution in the spatial domain of the phase of the FFT 
%%% of the original image
aerialviewfftphase = unwrap(angle(aerialviewfft));
aerialview2phase = ifft2(exp(j*aerialviewfftphase));

figure(3);
imagesc(real(aerialview2phase));
colormap(graymap);
title('First image corresponding to phase of FFT of the image: real component');

figure(4);
imagesc(imag(aerialview2phase));
colormap(graymap);
title('Second image corresponding to phase of FFT of the image: imaginary component');
