%%%
%2D Fourier transform Circle
%%%
Nx =   64;					      %   number   of	grid points in x
Ny =   64;					      %   number   of	grid points in y
Nr =   100;					      %   number   of	detector elements
Nphi   = 500;					  %   number   of	projection angles

xLen =   2.0;					      %   side of reconstructed rectangle
yLen =   2.0;					      %   side of reconstructed rectangle
phiLen   = pi;  				      %   angular interval
rLen =   xLen * 1.5;				  %   length of detector

oversampling = 3;                 % anti-aliasing factor
sphereRadius = 0.01

indexVec_x = (1:Nx)-33;indexVec_y = (1:Ny)-33
dx   = 2/Nx;dy=2/Ny;
xVec = indexVec_x * dx; yVec = indexVec_y * dy;
dfx  = 1/(dx*Nx);dfy  = 1/(dy*Ny);
XVec = indexVec_x * dfx; YVec = indexVec_y * dfy;

   
%%%
% PHANTOM DEFINITION
%%%
E = [ 0.0   0.0     sphereRadius    1.0];

%comment next line if you want to reuse plots
figure(1);clf
phm = pixelize(E, xVec, yVec, oversampling);
subplot(2,2,1); imagesc(yVec, xVec, phm); title('Phantom');
axis image; colorbar; colormap gray; pause(0);

% Digital Fourier transform
f_phm =(abs(fftshift(fft2(phm))));
subplot(2,2,2); imagesc(YVec, XVec, f_phm); title('Digital FT ');
axis image; colorbar; colormap gray; pause(0);

% Continuous Fourier transform
a=2.*pi
[X, Y] = meshgrid(linspace(-dfx*Nx/2,dfx*Nx/2,Nx),linspace(-dfy*Ny/2,dfy*Ny/2,Ny));
R = sqrt(X.^2+Y.^2);
f = abs(a*sphereRadius*sphereRadius*besselj(1,sphereRadius*R*a)./(sphereRadius*R));
subplot(2,2,3); imagesc(YVec, XVec,f); title('Continuous FT');
%mesh(X,Y,reshape(f,size(X)));
axis image;

