% A program that performs spatial filtering.
% Written by Amita B. Deb. Oxford 2008.

clc;
clear all;
close all;


M = 2^nextpow2(2000);
N = 2^nextpow2(2000);
deltax = 1;
deltay = 1;

x = -floor((M-1)/2):deltax:floor(M/2);
y = -floor((N-1)/2):deltay:floor(N/2);

[X,Y] = meshgrid(x,y);
d = 20; % Width of the Gaussian.

Z = 0;

sites = 10; % Number of sites in one direction.

for count1=1:1:sites,
    for count2 = 1:1:sites,
        
        Q = exp(-(( X - M/2 + count1*M/(sites+1)).^2 + (Y - N/2 + count2*M/(sites+1)).^2)/d^2);
        %% Defects.
        %if count1==4 & count2==4, Q = 0*exp(-(( X - M/2 + count1*M/(sites+1)).^2 + (Y - N/2 + count2*M/(sites+1)).^2)/d^2); end
        %% Defects added.
        Z = Z + Q;
        
    end
end

Z = Z;



%colormap gray;
subplot(2, 2, 1);
imagesc(x,y,Z);
colorbar;
axis square;
set(gca, 'TickDir', 'out');
title('Input function');
xlabel('x');
ylabel('y');





kx1 = mod( 1/2 + (0:(M-1))/M , 1 ) - 1/2;
kx = kx1*(2*pi*deltax);
kx = sort(kx);


ky1 = mod( 1/2 + (0:(N-1))/N , 1 ) - 1/2;
ky = ky1*(2*pi*deltay);
ky = sort(ky);


%% Fourier Transform.

A = fftshift(fft2(Z,M,N));
subplot(2, 2, 2);
imagesc(kx,ky,abs(A));
colorbar;
axis square;
set(gca, 'TickDir', 'out');
title('Fourier (focal) plane');
imwrite(A,'3.filtered.jpg');


%% Define the filter on the Fourier plane.

[KX,KY] = meshgrid(kx,ky);

% Gaussain filter
%Filter = exp((-(KX - 0).^2 + -(KY - 0).^2)/0.015^2);
%%%%%%%%%%%%%%%%%%%

% toroidal filter (between the zeroth and first order)
wew = length(kx);
Qw = zeros(wew,wew);

for count1 = 1:1:length(kx),
   for count2 = 1:1:length(ky),
       if (kx(count1)^2 + ky(count2)^2)>(.003)^2&&(kx(count1)^2 + ky(count2)^2)<(.03)^2,
           Qw(count1,count2) = 1;
       end
   end
end
Filter = Qw;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% No filter;
Filter = 1;
%%%%%%%%%%%%%%

subplot(2, 2, 3);
imagesc(kx,ky,Filter);
colorbar;
axis square;
set(gca, 'TickDir', 'out');
title('The filter');
imwrite(A,'4.filtered.jpg');

A_filtered = Filter.*A;
L_filtered = real(ifft2(ifftshift(A_filtered)));
subplot(2, 2,4);
imagesc(x,y,abs(L_filtered));
colorbar;
axis square;
set(gca, 'TickDir', 'out');
title('Image plane');