The MATLAB code below is used for filtered backprojection (for $x$-rays). In this example program, it reconstructs an image of a simulated checker board, generated by the program. It uses Fourier transforms and Fourier slice theorem to reconstruct a 2D slice of a sample object.

Can anyone explain what the $z$-axis of the reconstructed image represents?

I would assume for $x$-ray tomography it represents the $x$-ray attenuation co-efficient as a cross-section.

I have applied this code to straight-line-flux 'magnetic induction tomography', using Helmholtz coils and an array of sensors; with rotation of the metallic sample, for each projection. The result I obtain for the $z$-axis does not seem to correlate with the measured values of the initial images (taken at each rotation). However, I do appear to obtain the cross-sectional shape of the sample. Any help on this would be appreciated.

% from: Parallel Implementation of the Filtered Back Projection Algorithm for Tomographic Imaging
% http://www.sv.vt.edu/xray_ct/parallel/Parallel_CT.html

% Or see: A. C. Kak and M. Slaney, “ 3 Algorithms for Reconstruction with Nondiffracting Sources, ” in
% Principles of Computerized Tomographic Imaging, Society of Industrial and Applied
% Mathematics, 2001, pp. 203–273.
% http://www.slaney.org/pct/pct-toc.html

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

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;

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


% 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;
            blank = 0;

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

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



% 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


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.  


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


               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);
                 A(i,lb) = A(i,lb) + ((1-frac)*object(x,y));
                 A(i,hb) = A(i,hb) + (frac*object(x,y));


disp('forward projections completed');

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

% 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),

% This is the end of the loop on angles

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

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

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


disp('backward projection completed');

disp('store reconstructed image to f.txt');

%save f.txt


% making positional data
c = 0;
for i = 1:64
    for j = 1:64
        c = c + 1;
        position(c, 1) = j - 1;   % x's
        position(c, 2) = i - 1;  % y's 
dlmwrite('C:\Users\MIT1\Documents\Brendan_Darrer\UCL_AMOPP\MATLAB\CAT_and_backprojection\parallel_beam_x-ray\position.txt', position, 'delimiter', '\t', ...
         'precision', 6)

%make 2D arrays into one row array
p0 = f';
f2 = p0(:)'   % make 2D array A into single row

f3 = f2'   %make into single column

result = [position f3]

dlmwrite('C:\Users\MIT1\Documents\Brendan_Darrer\UCL_AMOPP\MATLAB\CAT_and_backprojection\parallel_beam_x-ray\result.txt', result, 'delimiter', '\t', ...
         'precision', 6)

x=result(:,1); y=result(:,2); z=result(:,3);

Fig2Handle = figure('Position', [100, 100, 1049, 895]);

xlin=linspace(min(x),max(x),100);  % was 50
ylin=linspace(min(y),max(y),100); % was 50

surf(X,Y,Z) % interpolated

axis tight; hold on


% Remove gridlines.
shading flat
shading interp

colormap hsv


title({'Filtered back projection of a';'simulated check borad - from 16 rays and 12 angles'})
xlabel('x axis / mm')
ylabel('y axis / mm')
ylabel(c,'z axis / mm')

filtered backprojection of simulated check board


2 Answers 2


The Z axis in your picture is a color map. It's a predefined map which maps gray-scale values to colors.

It has nothing to do with Computational tomography and used only for visualization.

Also in your experiment if your sensors measure summation of some parameter along different lines (like the attenuation coefficient in X-Ray tomography) from different views you could use these kind of reconstruction algorithms.


The cross-sectional image in x-ray tomography is a 2D model of the attenuation co-efficient $$\mu (x,y)$$

Therefore, whatever medium is being measured &/or calculated, from the projections, will be the medium in the reconstructed cross-sections.

E.g. "The decay of radioactive nucleoids in the body as in emission tomography, or the refractive index variations as in ultrasonic tomography." (1)

(1) A. C. Kak and M. Slaney, “4 Measurement of Projection Data - The Nondiffracting Case,” in Principles of Computerized Tomographic Imaging, Society of Industrial and Applied Mathematics, 2001, pp. 113.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.