%% Fourier filtering of an image file and example on Fraunhofer diffraction
% Toni Laurila 7.9.2012

clear all
close all

%% load file image file
load('C:\Users\Toni\Optiikan kurssi\Optiikka 2014\2014\circles.mat')

img=circles;
imgsize=size(img);
imgsize=imgsize(1);
figure(1)
imshow(img)
title('{\bf Original Image}')

%% Section A; Explain what happens here?

img_shift = fftshift(img);

img_transform=fft2(img_shift);
img_fft_intensity = img_transform.*conj(img_transform);
img_fft = fftshift(img_fft_intensity)/max(img_fft_intensity(:));

% next line enhances the image contrast with logarithmic amplification in order 
% to bring out details
img_fft = 1+log10(img_fft)/4;

figure(2)
imshow(img_fft)

% End of section A

%% Create aperture function for filtering
% size of the filter radius R (expressed as pixels)
R=100;

imgsize_power2=log2(size(img_fft))
n = 2^imgsize_power2(1);
F = zeros(n);
I = 1:n;
x = I-n/2;
y = n/2-I;
[X,Y] = meshgrid(x,y);
A = (X.^2 + Y.^2 <= R^2);
F(A) = 1;

%% Section B: Explain what is F & Figure 3? What about Finv?

F_ones=ones(imgsize,imgsize)
Finv=F_ones-F;

figure(3)
imshow(F)
%imshow(Finv)

% End of section B

%% Fraunhofer diffraction pattern
% Demonstrate how to calculate diffraction pattern through the aperture F

D1 = fft2(F);
diffracted_intensity = fftshift(D1);

figure(4)
imagesc(abs(diffracted_intensity))
colormap(hot)
title('{\bf Diffraction Pattern}')


%% Section C: Explain what is done here?

gF = fft2(img_shift);

gFfilt=fftshift(F).*real(gF)+i*imag(gF);

gBfilt = ifft2(gFfilt);

gB= fftshift(gBfilt);

gB=gB./max(gB(:));
gB=gB.*conj(gB);

figure(5)
imshow(gB)

% End of section C