% Fourier pairs demo
% Matlab examples for cs484-f13
% Richard Mann, Jan 30, 2013

%Im=double(imread('microserf.tif','tiff'));
Im=double(imread('einstein.pgm','pgm'));

close all;
figure; colormap gray; imagesc(Im); axis equal; title('original');

% All filters will have size 256x256
n=size(Im); n=n(1); halfn=n/2;
xi=-halfn:(halfn-1); % full range of window
xx= repmat(xi,length(xi),1);
yy= repmat(xi',1,length(xi));
rr= xx.^2+yy.^2;

fbox0=zeros(n); fbox0(halfnp1,halfnp1)=1;
fbox1=zeros(n); fbox1(halfnp1-1:halfnp1+1, halfnp1-1:halfnp1+1)=1;
fbox5=zeros(n); fbox5(halfnp1-2:halfnp1+2, halfnp1-2:halfnp1+2)=1;
fbox9=zeros(n);fbox9(halfnp1-4:halfnp1+4, halfnp1-4:halfnp1+4)=1;
fcirc5=zeros(n); fcirc5(rr<=25)=1;
fcirc10=zeros(n); fcirc10(rr<=10^2)=1;
Sigma=8; fgauss8=1/(2*pi)/(Sigma^2)*exp(-0.5*rr/(Sigma*Sigma));
Sigma=2; fgauss2=1/(2*pi)/(Sigma^2)*exp(-0.5*rr/(Sigma*Sigma));
Sigma=1; fgauss1=1/(2*pi)/(Sigma^2)*exp(-0.5*rr/(Sigma*Sigma));

N=256; k=[10, 3]; k1=k(1);k2=k(2);
Sigma=16; Sigma2=Sigma^2;
fgabor1=exp(j*2*pi*(k1*xx+k2*yy)/N).*exp(-rr/2/Sigma2);

N=256; k=[10,3]; k1=k(1);k2=k(2);
Sigma=8; Sigma2=Sigma^2;
fgabor2=exp(j*2*pi*(k1*xx+k2*yy)/N).*exp(-rr/2/Sigma2);

N=256;
Sigma=16; Sigma2=Sigma^2;
fdGauss1=exp(-rr/2/Sigma2);
fdGauss1=diff(fdGauss1);

% Display pairs: ff=input filter
ff=fbox1;
ff=fcirc5;
ff=fcirc10;
ff=fgauss1;
ff=fgabor1;
ff=fgabor1; ff=real(ff);
ff=fgabor1; ff=imag(ff);

ff=fgabor2;
ff=fgabor2; ff=real(ff);
ff=fgabor2; ff=imag(ff);

ff=fdgauss1;

% FFT and display
% Note: ff is filter, with center pixel at n/2+1,n/2+1.
% FFTSHIFT moves it to (zero,zero) (with wrap around)
FF=fftshift(fft2(fftshift(ff))); FF=FF./max(abs(FF(:)));
ff=ff./max(ff(:));
figure;
imagesc([ real(ff), imag(ff), abs(ff); real(FF), imag(FF), abs(FF)]); colormap gray;
axis equal; axis off;
title('[ Re(f), Im(f), Abs(f); Re(F), Im(F), Abs(F) ]','FontSize',18);
set(gcf,'PaperOrientation','landscape');
set(gcf,'PaperUnits','normalized');
set(gcf,'PaperPosition', [0 0 1 1]);
set(gcf,'DefaultAxesFontSize',16)
print(gcf, '-dpdf', 'Output.pdf');

figure; mesh(abs(FF));
set(gcf,'PaperUnits','normalized');
set(gcf,'PaperPosition', [0 0 1 1]);
set(gcf,'DefaultAxesFontSize',16)
print(gcf, '-dpdf', 'Output_mesh.pdf');

% Filtering example, Box vs Gaussian filter
ff=fgauss1;
ff=fbox1;

figure;
subplot(1,4,1); imagesc(Im); colormap gray; axis equal; axis off;
Imf=conv2(Im,ff,'same'); Imfs=Im(1:2:n,1:2:n);
subplot(1,4,2); imagesc(Imfs); colormap gray; axis equal; axis off;
Imfsf=conv2(Imfs,ff,'same'); Imfsfs=Imfsf(1:2:halfn, 1:2:halfn);
subplot(1,4,3); imagesc(Imfsfs); colormap gray; axis equal; axis off;
Imfsfsf=conv2(Imfsfs,ff,'same'); Imfsfsfs=Imfsfsf(1:2:halfn/2, 1:2:halfn/2);
subplot(1,4,4); imagesc(Imfsfsfs); colormap gray; axis equal; axis off;
set(gcf,'PaperOrientation','landscape');
set(gcf,'PaperUnits','normalized');
set(gcf,'PaperPosition', [0 0 1 1]);
set(gcf,'DefaultAxesFontSize',16)
print(gcf, '-dpdf', 'Output.pdf');

% "Ideal" filter, zero out components of FFT above some frequency
width=n;
hanning2d=0.5*(1-cos(2*pi*[0:width-1]'/(width-1)))*0.5*(1-cos(2*pi*[0:width-1]/(width-1)));
ImF=fftshift(fft2(Im));
ImF(rr>20^2)=0;
Imfilt=real(ifft2(fftshift(ImF)));
figure;
imagesc([Im,Imfilt]); colormap gray; axis equal; axis off;
set(gcf,'PaperOrientation','landscape');
set(gcf,'PaperUnits','normalized');
set(gcf,'PaperPosition', [0 0 1 1]);
set(gcf,'DefaultAxesFontSize',16)
print(gcf, '-dpdf', 'Output.pdf');