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;