Matlab Program Listing

% This is matlab code which does the ct reconstruction for
% parallel beam x-ray tomography. The code has been written to verify
% the correctness of the fortran algorithm. 
% The program will give an image pertaining to the projection 
% data by making gray scaled image of the final data.
% It takes about 15 minutes on a Sun Sparc 10 for an image of 64x64 pixels.

clear; tau=1.0;  rays=16;  angles=128;  flops(0);

% In this program the variable, tau is set to unity, but is defined
% in more general terms on pgs. 71-72 (Fig 3.14) "Kak and Slaney".

mp=rays/2.0;
x_min=1 
x_max=rays;
x_cen=((x_max)/2.0)
y_min=1 
y_max=rays;
y_cen=((y_max)/2.0)
x_max2 = 64;
y_max2= 64;
x_cen2 = x_max2/2;
y_cen2 = y_max2/2;
pi = acos(-1.0)
disp('Calculations in progress...........');

%initialize the object, the original pattern, to zeros
disp('initlalizing the object');
 for i = 1:x_max,
   for j = 1:y_max,
     object(i,j) = 0;
   end;
 end;

% initalize the projection matrix to zero: 
% simulated parallel xray projections [eqn(3) pg. 50]
disp('initlalizing the projection matrix');
 for i = 1:angles,
   for j = 1:rays,
     A(i,j) = 0;
   end;
 end;

nf1=flops;

% This routine generates data of the original object: (x_max, y_max),
% The pattern is supposed to look like a checker board. 
% Make a checker board pattern, see Fig. 4.6 pg. 36 of this report.
  blank =0;
    for i =1:4:x_max,
          if blank == 0 blank = 255;
          else
            blank = 0;
          end;

      for j=1:4:y_max,
              if blank ==0 blank = 255;
              else
                blank = 0;
              end;

                  for i1=i:i+4,
                      for j1=j:j+4,
                         object(i1,j1) = blank;
                      end;
                  end;

      end;
    end;

nf2=flops;

% Chris Henze's ramp function: (Biology)
disp('generating the ramp');
a1 = linspace(0,32,32);             % create linear ramp, pg 168 MATLAB
hf = [a1 fliplr(a1(1:32))];          % flip left-to-right, pg 77 MATLAB 
H = hf .* hamming(64)';           % Hamming funct, pg 168 MATLAB

nf3=flops;

disp('Using Hamming Window N =64');
% projecting the data points. The projections are 
% 64 angles X 16 rays per angle. 

disp('Projecting the object');
% The factor=1.15 in the denominator requires some explanation.
% Ideally this should be just 1/sqrt(2).  This factor,
% 1/sqrt(2), is the worst case at a rotation of 45 degrees where
% the object always fits inside a square whose edge dimension is
% the largest dimension of the original object.  Pragmatically,
% this factor 1.15 has been put in the denominator so that rays
% always fit within the object.  When the rays fall outside of the
% object, although the result should be near zero, a numerical
% error results.  If 1.15 is changed to 1.0 this program will not
% run.  An alternate technique may work if the uneven object fits
% into a circle instead of a square, where the largest object 
% dimension becomes the diameter of the circle.  

tw=1/(1.15*sqrt(2));

% Calculate the forward projections: Kak & Slaney pgs. 49-56
% This is where the CM5 does not parallelize well, but we can
% use the Paragon to parallelize this section by assigning 
% rows(angles) to each cpu(node).

  for x=x_min:x_max,
     for y=y_min:y_max,

               sum=0.0;

               for i=1:angles,

                   theta = (i-1)*pi/angles;
r = cos(theta)*((x-x_cen)/(x_cen))+sin(theta)*((y-y_cen)/(y_cen));
                       mb = (r*mp*tw)+mp;
                       lb = floor(mb);hb=ceil(mb);frac=mb-lb;
                   A(i,lb) = A(i,lb) + ((1-frac)*object(x,y));
                   A(i,hb) = A(i,hb) + (frac*object(x,y));
               end;
     end;
  end;

disp('forward projections completed');
nf4=flops;

% Calculate the 1D FFT, Step 1, pg 15 Report, also eqn(2.4)
for i=1:angles,
    for j=1:rays,
      p1(j)=A(i,j);
    end;

% pad with zeros to the right of the matrix p1
      p = [p1 zeros(1,48)];  % Extend p1(16) to 64 by padding with zeros

      a2 = fft(p);

% The next three lines could be eliminated but they
% were included here to view intermediate values
%   for j=1:length(a2),
%     A2(i,j)=a2(j);
%   end;

% Recall, H, Hamming operation, see the inside integral of eqn(33) pg 64
% of Kak & Slaney.  Also see this same operation on the bottom line of
% eqn(69) on pg 75, H = [FFT h (n t) with ZP]xsmoothing-window}.
  dtime= fft(p) .* H;

% Again the next three lines could be eliminated but they
% were included here to view intermediate values
%  a3 = dtime;
%  for j=1:length(a3),
%    A3(i,j)=a3(j);
%  end;

% The left operation of eqn(69)
  d = ifft(dtime);
  for j=1:length(d),
    c(i,j)=d(j);
  end;

% This is the end of the loop on angles
end;

disp('Filtered Projections Completed');
nf5 = flops;

disp('backprojecting...');
% This is the final step of backprojecting the results into array f(x,y).
% This is the summation of step 4 eqn(2.21) on pg 15 in this report, 
% also eqn(45) on pg 67 of Kak & Slaney.  
% This operation is similar to the forward projection
% except now we have to sum (this sum can be divided over cpu's) instead of 
% the arbitrary "random" accumulation of values of the forward projection.
% Comment:  This part is important for both the CM5 & the Paragon: each
% of these machines parallelize this sum differently but this is the most 
% important part of the total time on either machine and represent the 
% section that benefits the most from the parallelization on the CM5 or 
% the Paragon. 

  for x=x_min:x_max2,
   for y=y_min:y_max2,

         sum=0.0;
         for i=1:angles,
             theta = (i-1)*pi/angles;
r = cos(theta)*((x-x_cen2)/(x_cen2))+sin(theta)*((y-y_cen2)/(y_cen2));
                 mb = (r*mp*tw)+mp;
                 lb = floor(mb);hb=ceil(mb);frac=mb-lb;
             sum = sum + real(c(i,lb))*(1-frac) + frac*real(c(i,hb));
          end;

     f(x,y)=sum;

   end;
  end;
disp('backward projection completed');

disp('store reconstructed image to f.dat');
save f.dat f /ascii

disp('DONE')

%fid = fopen('mag.gray','wb');
%fwrite(fid,fgray,'real');
%nf6=flops;
%save ct26_0117.mat;
%quit;