Direct Fourier Reconstruction of a Tomographic Slice
Experiments of reconstruction using Fourier Slice Theorem (rather than filtered back projection, FBP).
The implementation reconstructs a tomographic image (i.e., a slice) using the Fourier Slice Theorem (also referred to as Fourier Projection Theorem or Central Slice Theorem) rather than the most common Filtered Back Projection (FBP) method. The Fourier Slice Theorem (FST) holds for parallel X-ray beams and does not hold for divergent sources. The Direct Fourier Reconstruction (DFR) code uses a phantom image, computes its radon transform (i.e., builds the image sinogram), symmetrically zero-pads the sinogram columns, computes the fft1 of the projections (on sinogram columns), filters in frequency the fft1 (with a selectable low-pass window) and allocates the filtered fft1s radially in order to build the (pseudo) fft2 of the original phantom. Finally, the phantom is recovered (reconstructed) via ifft2. The code explores different interpolation methods for the radial/polar lay-out of the fft1. The polar raster resulting from the allocation of the DFT coefficients is sparse. The crucial section of the Direct Fourier Reconstruction (DFR) method is the binning (gridding) of the fft1 Fourier coefficients of the sinogram columns (projection angle theta, sensor pixel position) in the Cartesian grid (omega, zeta) of the fft2. This is accomplished via interpolation of scattered fft1 data into a Cartesian grid that should be fully populated. This requires also extrapolation in addition to interpolation.
We test different strategies and routines for interpolating Fourier coefficients in the frequency domain of the pseudo fft2. While in the real space a local interpolation error affects only the grey level of a few neighbour pixels, in the Fourier domain a coefficient impacts to all the pixels of the image recovered with the ifft2 transform giving rise to undesirable artifacts.
If the projection is zero padded before fft1, the resulting reconstructed slice is cropped down (i.e. unpadded) to the original size of the phantom image. Throughout the code the DC coefficient is kept to the centre of the Fourier lines. High frequency artifact in the reconstructed slice ( when present !) are mitigated by zeroing the ’corners’ of the fft2 array before calling the ifft2 inverse transform ( see function Ifft2_2_Img). In other terms the 2D Fourier samples are set to zero with a low pass filter (that is referred to as 'mask' in the section of the code that recovers the image slice via inverse 2D Fourier transformation). Alternatively to the binary mask (i.e. ideal cutter of high frequencies) a 2D Butterworth filter is available by editing the code.
The FBP reconstruction (see the matlab function iradon) does a "smearing back" of the (low-pass filtered) projection across the pixel grid of the image at the angle the projection was acquired by the detector. This is repeated for each projection and each projection contribution summed up.; the FBP computational complexity (CC) is O(N^3).
The computational complexity of the fft of a -N values- projection-vector is O(N log N). Thus the CC of the Direct Fourier Reconstruction (DFR) is O(N^2 log N). Consequently, in principle, the DFR is N/ log(N) more computationally efficient than the FBP, where N is length of the projection in pixels. Notice that -since DFR requires zero padding the projection- N is the size of the padded projection!
The practice is another thing; on one side, the advantage can be largely eroded by the (mentioned) computing intensive interpolation of the Fourier coefficients in 2D that is not accounted for in the comparison of complexities. The CC of the interpolation of the complex Fourier coefficients in the Cartesian 2d space depends on the selected interpolation scheme (i.e., NN, linear, cubic, spline, etc.). Nearest-neighbours or linear interpolation schemes are generally inadequate and higher order interpolation is preferable especially if the number of projection-views is below that dictated by the Nyquist principle. On the other side, fortunately, most of the computational work of the interpolation (e.g., triangulation of the set of N^2 coordinate points) is done once for all the slices in the stack-of-slices that one should reconstruct from projections acquired with 2d detectors to form a tomographic volume.
The reconstruction with FBP is also shown in the figures generated by the code.
The code requires the Image Processing Toolbox and the Statistics toolbox
Tested only with the 2016a version
At the prompt “>” type
or simply :
to run a test example.
Keywords: Tomography, Fourier slice theorem, Tomographic reconstruction
Gianni Schena (2021). Direct Fourier Reconstruction of a Tomographic Slice (https://www.mathworks.com/matlabcentral/fileexchange/60257-direct-fourier-reconstruction-of-a-tomographic-slice), MATLAB Central File Exchange. Retrieved .
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!