Image Processing Toolbox 

This example shows how to use radon, iradon, fanbeam, and ifanbeam to form projections from a sample image and then reconstruct the image from the projections. While radon and iradon use a parallelbeam geometry for the projections, fanbeam and ifanbeam use a fanbeam geometry. To compare parallelbeam and fanbeam geometries, the examples below create synthetic projections for each geometry and then use those synthetic projections to reconstruct the original image.
A realworld application that requires image reconstruction is Xray absorption tomography where projections are formed by measuring the attenuation of radiation that passes through a physical specimen at different angles. The original image can be thought of as a cross section through the specimen in which intensity values represent the density of the specimen. Projections are collected by special medical imaging devices and then an internal image of the specimen is reconstructed using iradon or ifanbeam.
The function iradon reconstructs an image from parallelbeam projections. In parallelbeam geometry, each projection is formed by combining a set of line integrals through an image at a specific angle. The function ifanbeam reconstructs an image from fanbeam projections, which have one emitter and multiple sensors.
See the Image Processing Toolbox™ User's Guide for diagrams that illustrate both geometries.
On this page… 

Parallel Beam  Calculate Synthetic Projections Parallel Beam  Reconstruct Head Phantom from Projection Data 
The test image is the SheppLogan head phantom which can be generated using the function phantom. The phantom image illustrates many qualities that are found in realworld tomographic imaging of human heads. The bright elliptical shell along the exterior is analogous to a skull and the many ellipses inside are analogous to brain features or tumors.
P = phantom(256); imshow(P)
Parallel Beam  Calculate Synthetic Projections
Calculate synthetic projections using parallelbeam geometry and vary the number of projection angles. For each of these calls to radon, the output is a matrix in which each column is the Radon transform for one of the angles in the corresponding theta.
theta1 = 0:10:170; [R1,~] = radon(P,theta1); num_angles_R1 = size(R1,2)
num_angles_R1 = 18
theta2 = 0:5:175; [R2,~] = radon(P,theta2); num_angles_R2 = size(R2,2)
num_angles_R2 = 36
theta3 = 0:2:178; [R3,xp] = radon(P,theta3); num_angles_R3 = size(R3,2)
num_angles_R3 = 90
Note that for each angle, the projection is computed at N points along the xpaxis, where N is a constant that depends on the diagonal distance of the image such that every pixel will be projected for all possible projection angles.
N_R1 = size(R1,1) N_R2 = size(R2,1) N_R3 = size(R3,1)
N_R1 = 367 N_R2 = 367 N_R3 = 367
So, if you use a smaller head phantom, the projection needs to be computed at fewer points along the xpaxis.
P_128 = phantom(128); [R_128,xp_128] = radon(P_128,theta1); N_128 = size(R_128,1)
N_128 = 185
Display the projection data R3. Some of the features of the original phantom image are visible in the image of R3. The first column of R3 corresponds to a projection at 0 degrees, which is integrating in the vertical direction. The centermost column corresponds to a projection at 90 degrees, which is integrating in the horizontal directions. The projection at 90 degrees has a wider profile than the projection at 0 degrees due to the large vertical semiaxis of the outermost ellipse of the phantom.
figure, imagesc(theta3,xp,R3) colormap(hot) colorbar xlabel('Parallel Rotation Angle  \theta (degrees)'); ylabel('Parallel Sensor Position  x\prime (pixels)');
Parallel Beam  Reconstruct Head Phantom from Projection Data
Match the parallel rotationincrement, dtheta, in each reconstruction with that used above to create the corresponding synthetic projections. In a realworld case, you would know the geometry of your transmitters and sensors, but not the source image, P.
The following three reconstructions (I1, I2, and I3) show the effect of varying the number of angles at which projections are made. For I1 and I2 some features that were visible in the original phantom are not clear. Specifically, look at the three ellipses at the bottom of each image. The result in I3 closely resembles the original image, P.
Notice the significant artifacts present in I1 and I2. To avoid these artifacts, use a larger number of angles.
% Constrain the output size of each reconstruction to be the same as the % size of the original image, P. output_size = max(size(P)); dtheta1 = theta1(2)  theta1(1); I1 = iradon(R1,dtheta1,output_size); figure, imshow(I1)
dtheta2 = theta2(2)  theta2(1); I2 = iradon(R2,dtheta2,output_size); figure, imshow(I2)
dtheta3 = theta3(2)  theta3(1); I3 = iradon(R3,dtheta3,output_size); figure, imshow(I3)
Fan Beam  Calculate Synthetic Projections
Calculate synthetic projections using fanbeam geometry and vary the 'FanSensorSpacing'.
D = 250; dsensor1 = 2; F1 = fanbeam(P,D,'FanSensorSpacing',dsensor1); dsensor2 = 1; F2 = fanbeam(P,D,'FanSensorSpacing',dsensor2); dsensor3 = 0.25; [F3, sensor_pos3, fan_rot_angles3] = fanbeam(P,D,... 'FanSensorSpacing',dsensor3);
Display the projection data F3. Notice that the fan rotation angles range from 0 to 360 degrees and the same patterns occur at an offset of 180 degrees because the same features are being sampled from both sides. You can correlate features in this image of fanbeam projections with the same features in the image of parallelbeam projections, above.
figure, imagesc(fan_rot_angles3, sensor_pos3, F3) colormap(hot) colorbar xlabel('Fan Rotation Angle (degrees)') ylabel('Fan Sensor Position (degrees)')
Fan Beam  Reconstruct Head Phantom from Projection Data
Match the fansensorspacing in each reconstruction with that used to create each of the synthetic projections. In a realworld case, you would know the geometry of your transmitters and sensors, but not the source image, P.
Changing the value of the 'FanSensorSpacing' effectively changes the number of sensors used at each rotation angle. For each of these fanbeam reconstructions, the same rotation angles are used. This is in contrast to the parallelbeam reconstructions which each used different rotation angles.
Note that 'FanSensorSpacing' is only one parameter of several that you can control when using fanbeam and ifanbeam. You can also convert back and forth between parallel and fanbeam projection data using the functions fan2para and para2fan.
Ifan1 = ifanbeam(F1,D,'FanSensorSpacing',dsensor1,'OutputSize',output_size); figure, imshow(Ifan1)
Ifan2 = ifanbeam(F2,D,'FanSensorSpacing',dsensor2,'OutputSize',output_size); figure, imshow(Ifan2)
Ifan3 = ifanbeam(F3,D,'FanSensorSpacing',dsensor3,'OutputSize',output_size); figure, imshow(Ifan3)