%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  SIMULATION AND RECONSTRUCTION IN PARALLEL BEAM TOMOGRAPHY  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


colormap(gray(256))

cj=sqrt(-1); pi2=2*pi; cc=pi/180;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%        SHEPP-LOGAN HEAD PHANTOM PARAMETERS   %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

io=10;      % Number of ellipsoids
            % a(i): Horizontal (x) axis size of ellipsoid
            % b(i): Vertiacl (y) axis size of ellipsoid
            % A(i): Ellipsoid intensity
            % x(i): Shift of ellipsoid in horizontal (x) domain
            % y(i): Shift of ellipsoid in vertical (y) domain
            % to(i): rotation of ellipsoid in the (x,y) domain

i=1; a(i)=.92; b(i)=.69; A(i)=2; x(i)=.0; y(i)=.0; to(i)=0;

i=2; a(i)=.874; b(i)=.6624; A(i)=-.7; x(i)=-.0184; y(i)=0; to(i)=0;

i=3; a(i)=.31; b(i)=.11; A(i)=-.4; x(i)=.0; y(i)=-.22; to(i)=-18*cc;

i=4; a(i)=.41; b(i)=.16; A(i)=-.4; x(i)=.0; y(i)=.22; to(i)=18*cc;

i=5; a(i)=.25; b(i)=.21; A(i)=.1; x(i)=.35; y(i)=.0; to(i)=0;

i=6; a(i)=.046; b(i)=.046; A(i)=.12; x(i)=.1; y(i)=.0; to(i)=0;

i=7; a(i)=.046; b(i)=.046; A(i)=.12; x(i)=-.1; y(i)=.0; to(i)=0;

i=8; a(i)=.023; b(i)=.046; A(i)=.12; x(i)=-.605; y(i)=.08; to(i)=0;

i=9; a(i)=.023; b(i)=.023; A(i)=.12; x(i)=-.605; y(i)=.0; to(i)=0;

i=10; a(i)=.046; b(i)=.023; A(i)=.12; x(i)=-.605; y(i)=-.06; to(i)=0;

c=.0025*ones(1,io);   % FORGET IT; c(i) is ellipsoid's axis in z domain
% c=zeros(1,io);

X0=1;                    % target's relative size
nx=128;                  % number of detectors (receivers)
du=(2*X0)/nx;            % detector sample spacing
u=du*(-nx/2:nx/2-1);     % detector array
dku=1/(nx*du);           % freq. domain detector sample spacing
ku=dku*(-nx/2:nx/2-1);   % freq. domain detector array
%
N=2*ceil(.8*nx);         % number of rotation angles (theta)
N2=N/2;                  % half of number of rotation angles
dt=pi/N;                 % sample spacing in angle domain (d_theta)
t=dt*(0:N-1);            % angle domain (theta) array
%
z=zeros(1,N);            % FORGET IT

U=ones(N,1)*u;           % U(theta,u)
TE=t(:)*ones(1,nx);      % Theta(theta,u)

s=zeros(N,nx);           % s(theta,u)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%         SIMULATION            %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:io;
 s=s+A(i)*P(a(i),b(i),c(i),z,TE,U,x(i),y(i),to(i));
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       ZERO-PADDING         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

u2=du*(-nx:nx-1);         % detector array (zero-padded)
dku2=1/(2*nx*du);         % freq. domain detector sample spacing
ku2=dku2*(-nx:nx-1);      % freq. domain detector array
temp=s;
s=zeros(N,2*nx);
s(:,nx/2+1:nx/2+nx)=temp;
clear temp

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      RECONSTRUCTION           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s1=fty(s);               % FT[s(theta,u)]=S(theta,ku)
 
KX=cos(t(:))*ku2;         % kx(theta,ku)
KY=sin(t(:))*ku2;         % ky(theta,ku)


% INTERPOLATION FROM POLAR TO RECTANGULAR

mx=12; kmx=dku2*mx;
F=zeros(2*nx+4*mx,2*nx+4*mx);         % Target function F(kx,ky)
nx2=nx+2*mx+1;
KU=dku2*(-nx-2*mx:nx-1+2*mx);
XX=2*(.96*2*X0);


for i=1:N;
 i
 for j=1:2*nx;
  ix=nx2+round(KX(i,j)/dku2);
  jy=nx2+round(KY(i,j)/dku2);
  I=ix-mx:ix+mx; J=jy-mx:jy+mx;
  DD=KU(I)-KX(i,j);
  sx=sinc(XX*DD).*(.54+.46*cos(pi*(DD/kmx)));
  DD=KU(J)-KY(i,j);
  sy=sinc(XX*DD).*(.54+.46*cos(pi*(DD/kmx)));
  F(I,J)=F(I,J)+s1(i,j)*abs(ku2(j))*sx(:)*sy;
 end;
end;

F=F(2*mx+1:2*mx+2*nx,2*mx+1:2*mx+2*nx);



f=(dt*dku*(XX^2))*F;
f(nx+1,nx+1)=s1(1,nx+1);
f=f*((2*nx*dku2)^2)*du;
f=iftx(ifty(f));                  % Target function f(x,y)
G=real(f(nx/2+1:nx/2+nx,nx/2+1:nx/2+nx));imgc;
image(u,u(nx:-1:1),cg*(G-ng))
axis('square')