% 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;