% Finding well-focused areas in a photo using 2D FFT.
%
% The image was prepared in the routine PrepareImage.m
% which needs to be run before this file.
%
% Samuli Siltanen Nov 2015

% Read in the image 
load cladonia im
[row,col] = size(im);

% Scale the image to the interval [0,1]
im = MyScaleIm(im);

% Patch size
N = 64;

% Examine specific patches of the image
patch1 = im(1:N,300+[1:N]);
patch2 = im(1:N,col-N+[1:N]);

% Fourier transform
Fpatch1 = fft2(patch1-mean(patch1(:)));
Fpatch2 = fft2(patch2-mean(patch2(:)));

% Build a frequency space representation good for 
% comparison
F1 = fftshift(Fpatch1);
%F1 = log(abs(F1));
F1 = abs(F1);
%F1 = MyScaleIm(F1);
F2 = fftshift(Fpatch2);
%F2 = log(abs(F2));
F2 = abs(F2);
%F2 = MyScaleIm(F2);

% Create a band-pass mask
t = [-N/2:(N/2-1)];
[X,Y] = meshgrid(t);
mask  = ones(size(patch1));
mask(sqrt(X.^2+Y.^2)<(N/2)/5)  = 0;
mask(sqrt(X.^2+Y.^2)>(N/2)*.8) = 0;

% Define quantitative measures for
% "focusedness"
mA1 = sum(F1(:));
mA2 = sum(F2(:));
mB1 = var(F1(:));
mB2 = var(F2(:));
mC1 = sum(mask(:).*F1(:));
mC2 = sum(mask(:).*F2(:));
disp([mA1,mB1,mC1;mA2,mB2,mC2])

figure(3)
clf
subplot(2,1,1)
imagesc([patch1,patch2])
axis equal
axis off
colormap gray
subplot(2,1,2)
F1(N/2+1,N/2+1) = NaN;
F2(N/2+1,N/2+1) = NaN;
imagesc([F1,F2])
axis equal
axis off
colormap gray


% Take a look
figure(1)
clf
colormap gray
subplot(2,2,1)
imagesc(patch1,[0,1])
axis equal
axis off
colormap gray
subplot(2,2,3)
imagesc(F1)
axis equal
axis off
colormap gray
subplot(2,2,2)
imagesc(patch2,[0,1])
axis equal
axis off
colormap gray
subplot(2,2,4)
imagesc(F2)
axis equal
axis off
colormap gray





% Take a look
figure(2)
clf
subplot(1,3,1)
imagesc(im,[0,1])
axis equal
axis off
colormap gray
subplot(1,3,2)
imagesc(patch1,[0,1])
axis equal
axis off
colormap gray
subplot(1,3,3)
imagesc(patch2,[0,1])
axis equal
axis off
colormap gray




