% Testing the fast Fourier transform
% in dimension 2 (for images)
%
% Samuli Siltanen November 2014

% Load image 
im = imread('Tokyo2048.tif','tif');
im = double(im);

% Downsample image
im = MyDS2(im); % 1024x1024
im = MyDS2(im); % 512x512
im = MyDS2(im); % 256x256

% Normalize image: pixel values between zero and one
im = im-min(im(:));
im = im/max(im(:));

% Compute the FFT
Fim = fft2(im);

% Reorganize the FFT elements
Fim = fftshift(Fim);

% Record the size of the image (We assume row=col)
[row,col] = size(Fim);

% Create coordinates for the Fourier transformed image
D = 1/col;
t = D/2 + [0:(col-1)]*D;
[X,Y] = meshgrid(t);

% Put some of the frequency components to zero
index        = sqrt((X-0.5).^2+(Y-0.5).^2)>1/20; 
Fim_filtered = Fim;
Fim_filtered(index) = 0;
Fim_filtered2 = Fim;
Fim_filtered2(~index) = 0;

% Calculate filtered image
%im_filtered  = real(ifft2(fftshift(Fim)));
im_filtered  = real(ifft2(fftshift(Fim_filtered)));
im_filtered = im_filtered-min(im_filtered(:));
im_filtered = im_filtered/max(im_filtered(:));
im_filtered2 = real(ifft2(fftshift(Fim_filtered2)));
im_filtered2 = im_filtered2-min(im_filtered2(:));
im_filtered2 = im_filtered2/max(im_filtered2(:));

% Prepare the FFT side for plotting
Fim_filtered_show = log(abs(Fim_filtered)+.1);
Fim_filtered_show = Fim_filtered_show-min(Fim_filtered_show(:));
Fim_filtered_show = Fim_filtered_show/max(Fim_filtered_show(:));
Fim_filtered2_show = log(abs(Fim_filtered2)+.1);
Fim_filtered2_show = Fim_filtered2_show-min(Fim_filtered2_show(:));
Fim_filtered2_show = Fim_filtered2_show/max(Fim_filtered2_show(:));

% Plot the original image and the low-pass filtered image
figure(1)
clf
imagesc([im,Fim_filtered_show,im_filtered])
colormap gray
axis equal
axis off

% Plot the original image and the high-pass filtered image
figure(2)
clf
imagesc([im,Fim_filtered2_show,im_filtered2])
colormap gray
axis equal
axis off


