Reconstruction from Projections
Reconstruction of a two-dimensional (2-D) function (image, in
our case), from its one-dimensional (1-D) projections at all
possible angles, is possible because of the Projection-Slice
theorem.
The 1-D Fourier transform of the projection of a 2-D function, g(x,y), at an angle
'w' is equal to the slice of the 2-D Fourier transform of the function
at an angle 'w'.
Projection-Slice Theorem
Figure 3.1 Projection-Slice Theorem
This means that if we know the 1-D projections of a 2-D function at all possible angles, we know the 2-D Fourier transform of the function which is equivalent to knowing the function. Therefore, we can reconstruct the function exactly if we know the projections at all possible angles.
This reconstruction of the original image from its projections can be done in two steps:
- Filtering
- Back-projection
The image can be reconstructed from the projections using the Radon inversion formula given by
 |
(3.1) |
where P(f) is the Fourier transform of the 1-D projection p(t).
Further details on the reconstruction algorithm can be found in references [3-5, 7].
Filtering the projection
It is clear from the inversion formula that the first step in the reconstruction (inner integral) is to filter each of the projections using a filter with frequency response |f|.
Figure 3.2 Filtering
Back-projection
After computing the filtered projections, they must be sampled at (x*cos(w) + y*sin(w)) for each 'w' and integrated over 'w'. In the MATLAB implementation, since we can have only a finite number of sampled projections, an additional interpolation step is required to obtain the sample for back-projection.
Figure 3.3 Back-projection
Implementation in MATLAB
The reconstruction is performed in MATLAB using the function backproj.m. This function uses two other functions, convfilter.m (or fftfilter.m) and bproj.m, to perform the filtering and back-projection respectively.
Convfilter.m implements the filtering operation by direct convolution while fftfilter.m implements the filtering operation using the Fast Fourier Transform (FFT). The impulse response of the filter needed to implement direct convolution in convfilter.m is obtained by analytically calculating the impulse response and sampling it at the detector spacing.
Bproj.m determines the reconstructed image pixel by pixel by back-projecting the filtered projections. To reconstruct each pixel, the sample from each filtered projection corresponding to the pixel is back-projected, and these samples are summed together. Since the filtered projections are available only in sampled form and the required sampling instants may lie in between the available sampling instants, the required samples for back-projection have to be obtained by interpolation. A linear interpolation scheme is used in our program.