%% Fourier Transform General Information

%see:   http://fftw.org/
%most fft codes use fftw libraries in the implmentation


%convolution theorm:
%   http://mathworld.wolfram.com/ConvolutionTheorem.html

%Fourier transforms are separable in dimensions



%Beware the placement of the 2pi! 
%if the FT is using exp(-i*w*t), then the 2pi is in the w value and must
%come into the normalization as well
%as opposed to exp(-i 2*pi* f*t)
%see:   http://en.wikipedia.org/wiki/Fourier_transform

%Note, we are really doing discrete fourier transforms
%   http://en.wikipedia.org/wiki/Discrete_Fourier_transform


%we can program a discrete fourier transform ourselves:

%we'll use the definition that 'k' isn't really 2pi/lambda, but a spatial
%'frequency', so FT is using exp(-i*2pi*k*x)

%the maximum 'frequency', is then 'N/L', where 'N' is the number of
%elements and 'L' is the physical size (length units) 

%In the discrete transform the frequency -k/n is the same as the frequency
%(n-k)/n.


%% For illustration purposes, make direct calulations of DFTs


% some different patterns to try 

a=cos(pi/2*xx(16,16))
%a=sin(pi*xx(16,16))
%a=cos(pi*xx(16,16))
%a=noise(10*(rr(a)<5))
%a=newim([16 16])+1

% ZERO FREQUENCY AT CENTER (Dipimage)
x=xx(a);
y=yy(a);
% FT in x direction
ft_a_x=newim([16 16],'complex');
k_im=x/size(a,1); %this will define 'k'

for ii=0:size(a,1)-1
    k=k_im(ii,0);
    ft_a_x(ii,:)=sum(exp(-i*2*pi*k*x).*a,[],[1]);
end
ft_a_x;

% FT in y direction
ft_a=newim([16 16],'complex');
k_im=y/size(a,2); %this will define 'k'

for jj=0:size(a,2)-1
    k=k_im(0,jj);
    ft_a(:,jj)=sum(exp(-i*2*pi*k*y).*ft_a_x,[],[2]);
end
ft_a;

%normalize to give us values in ft_a that correspond to the amplitude of
%the waves
ft_a=ft_a/size(a,1)/size(a,2)


% ZERO FREQUENCY AT TOP LEFT (MATLAB)

x=xx(a)-min(xx(a));
y=yy(a)-min(yy(a));
% FT in x direction
ft_a_x=newim([16 16],'complex');
k_im=x/size(a,1); %this will define 'k'

for ii=0:size(a,1)-1
    k=k_im(ii,0);
    ft_a_x(ii,:)=sum(exp(-i*2*pi*k*x).*a,[],[1]);
end
ft_a_x;

% FT in y direction
ft_a=newim([16 16],'complex');
k_im=y/size(a,2); %this will define 'k'

for jj=0:size(a,2)-1
    k=k_im(0,jj);
    ft_a(:,jj)=sum(exp(-i*2*pi*k*y).*ft_a_x,[],[2]);
end
ft_a;

%normalize to give us values in ft_a that coorespond to the amplitude of
%the waves
ft_a=ft_a/size(a,1)/size(a,2)



%% Dipimage Fourier transforms

k=pi/2;
sin_kx=sin(k*xx)

% 'ft' is a quick way to do forward transforms in dipimage in all
% dimensions
ft_sin_kx=ft(sin_kx)

%Beware normalization conventions !!!

%dipimage has the folowing normalization convention:

%F( f(x) )=1/sqrt(N) Sum f(x)*exp(-i*k*x)
%where 'N' is the number of elements in the transform direction.

ft_sin_kx/sqrt(size(ft_sin_kx,1))/sqrt(size(ft_sin_kx,2))

ft_sin_kx=dip_fouriertransform(sin_kx,'forward',[1 1])

sin_kx_2=dip_fouriertransform(ft_sin_kx,'inverse',[1 1])



%% MATLAB FTs

a=cos(pi*xx(16,16))
fft2_a=fft2(single(a));

dipshow(fft2_a)
dipshow(ifft2(fft2_a))

b=noise(10*(rr(a)<5))
fft2_b=fft2(single(b))

dipshow(fft2_b)
dipshow(ifft2(fft2_b))



%% Some typical applications

%convolution
kern=(rr<5)/sum(rr<5)
im=dip_image(abs(xx)<10&abs(yy)<20,'sfloat')

blur=real(ift(ft(im).*ft(kern)))

%de-convolution?
deblur=real(ift(ft(im)./ft(kern)))
%doesn't work- look at 1./ft(kern)

%image shift
shiftvector=[5.2 10.2]
sz=size(im)
sx=xx(im)*shiftvector(1)/sz(1)*2*pi;
sy=yy(im)*shiftvector(2)/sz(2)*2*pi;

phase=exp(-i*(sx+sy));
shiftim=real(ift(ft(im).*phase))
joinchannels('RGB',im,shiftim)

%auto-correlation functions
blur=real(ift(conj(ft(im)).*ft(im)))

% a quick way to calculate an Airy pattern
c=rr<25
psf=real(ft(c).*conj(ft(c)))


%% Finding Fit Regions for SPT (No FFT relevance)




