Appendix B: Matlab code for the filtered backprojection
algorithm 

Below we give a few MATLAB M-files which implement the parallel-beam 
 filtered backprojection algorithm for the standard lattice. 
 This source code is provided for illustrative purposes only and 
comes without warranties of any kind.
The main 
file fbp.m is supplemented by three function M-files.  Rad.m
computes line integrals for a mathematical phantom consisting of 
ellipses, slkernel.m computes the discrete Shepp-Logan 
convolution kernel (cf. (\ref{eq:kSL})), 
and window3.m allows to view the reconstructed image. 
It is automatically called at the end of the reconstruction.
However, with the example phantom given in fbp.m the picture shown 
does not display the most interesting details. 
It is better to 
 call window3 again with the parameters window3(-0.07,0.07,roi,P).
 
% --------------- cut here for M-file fbp.m -----------------------

%fbp.m
%Parallel-beam filtered backprojection algorithm
% for the standard lattice
% Last revision: August 29, 2001

%specify input parameters here 

p=200; %number of view angles between 0 and pi
q=64; %q=1/d, d = detector spacing
MX=128;  MY = 128; %matrix dimensions 
roi=[-1  1 -1  1]; %roi=[xmin xmax ymin ymax]
                   %region of interest where 
                   %reconstruction is computed 
circle = 1; % If circle = 1 image computed only inside 
            % circle inscribed in roi.


%Specify parameters of ellipses for mathematical phantom. 
% xe = vector of x-coordinates of centers 
% ye = vector of y-coordinates of centers 
% ae = vector of first half axes
% be = vector of second half axes
% alpha = vector of rotation angles (degrees)
% rho = vector of densities 
xe=[0 0 0.22 -0.22 0 0 0 -0.08 0 0.06 0.5538];
ye=[0 -0.0184 0 0 0.35 0.1 -0.1 -0.605 -0.605 -0.605 -0.3858];
ae=[0.69 0.6624 0.11 0.16 0.21 0.046 0.046 0.046 0.023 0.023 0.0333];
be=[0.92 0.874 0.31 0.41 0.25 0.046 0.046 0.023 0.023 0.046 0.206];
alpha=[0 0 -18 18 0 0 0 0 0 0 -18];
rho=  [1 -0.98 -0.02 -0.02 0.01 0.01 0.01 0.01 0.01 0.01 0.03];

%end of input section 

b=pi*q; rps=1/b;
alpha = alpha*pi/180;

if MX > 1 
 hx = (roi(2)-roi(1))/(MX-1);
 xrange = roi(1) + hx*[0:MX-1];
else 
 hx = 0; xrange = roi(1);
end

if MY > 1 
 hy = (roi(4)-roi(3))/(MY-1);
 yrange = flipud((roi(3) + hy*[0:MY-1])');
else 
 hx = 0; yrange = roi(3);
end

center = [(roi(1)+roi(2)), (roi(3)+roi(4))]/2;
x1 = ones(MY,1)*xrange; %x-coordinate matrix
x2 = yrange*ones(1,MX); %y-coordinate matrix
if circle == 1
  re = min([roi(2)-roi(1),roi(4)-roi(3)])/2;
  chi = ((x1-center(1)).^2 + (x2-center(2)).^2 <= re^2); %char. fct of roi;
else
  chi = isfinite(x1);
end
x1 = x1(chi); x2 = x2(chi);
P = zeros(MY,MX);Pchi = P(chi);


h = 1/q;		
s = h*[-q:q-1];
bs = [-2*q:2*q-1]/(q*rps);
wb = slkernel(bs)/(rps^2); %compute discrete convolution kernel.
for j = 1:p
  j
  phi = (pi*(j-1)/p);
  theta = [cos(phi);sin(phi)];
  RF = Rad(theta,phi,s,xe,ye,ae,be,alpha,rho); %compute line integrals

%    convolution 
  C = conv(RF,wb);
  Q = h*C(2*q+1:4*q); Q(2*q+1)=0; 
  

%    interpolation and backprojection 

     Q = [real(Q)'; 0];
     t = theta(1)*x1 + theta(2)*x2;
     k1 = floor(t/h);
     u = (t/h-k1);
     k = max(1,k1+q+1); k = min(k,2*q);
     Pupdate =  ((1-u).*Q(k)+u.*Q(k+1));
     Pchi=Pchi+ Pupdate;

end   % j-loop

P(chi) = Pchi*(2*pi/p);

pmin = min(min(P)); 
pmax = max(max(P));
window3(pmin,pmax,roi,P); % view the computed image


% --- Cut here for Rad.m 

function [RF] = Rad(theta,phi,s,x,y,u,v,alpha,rho)

%  This function computes the Radon transform of ellipses
% centered at (x,y) with major axis u, minor axis v,
% rotated through angle alpha, with weight rho.


RF = zeros(size(s));

for mu = 1:max(size(x))
   a = (u(mu)*cos(phi-alpha(mu)))^2+(v(mu)*sin(phi-alpha(mu)))^2;
   test = a-(s-[x(mu);y(mu)]'*theta).^2;
   ind = test>0;
   RF(ind) = RF(ind)+rho(mu)*(2*u(mu)*v(mu)*sqrt(test(ind)))/a;
end  % mu-loop


% --- Cut here for slkernel.m 

function u = slkernel(t)

u = zeros(size(t));
i1 = abs(abs(t)-pi/2)<=1.e-6;
u(i1) = ones(size(u(i1)))/pi;
t1 = t(abs(abs(t)-pi/2)>1.e-6);
v = (pi/2 - t1.*sin(t1))./((pi/2)^2 - t1.^2);
u(abs(abs(t)-pi/2)>1.e-6) = v;
u = u/(2*pi^3);

% --- Cut here for window3.m

function pic1 = window3(mi,ma,roi,pic);
%function pic1 = window3(mi,ma,roi,pic);
% displays image pic with coordinates given by roi
% roi = [xmin xmax ymin ymax]
x = [roi(1), roi(2)]; y = [roi(3), roi(4)];
colors = 128; co = colors-1;
pic1 = pic - mi*ones(size(pic));
pic1 = (co/(ma-mi))*pic1;
P = (pic1 >= 0);
pic1 = pic1.*P;
P = (pic1 <= co);
pic1 = pic1.*P + co*(ones(size(pic1)) - P);
colormap(gray(colors));
image(x,fliplr(y),flipud(pic1));
axis('square');