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 Xray 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 zeropads the sinogram columns, computes the fft1 of the projections (on sinogram columns), filters in frequency the fft1 (with a selectable lowpass 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 layout 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 (lowpass 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 projectionvector 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.). Nearestneighbours or linear interpolation schemes are generally inadequate and higher order interpolation is preferable especially if the number of projectionviews 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 stackofslices 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
Direct_Fourier_Reconst(phantom(256),(0:1:179),3,1);
or simply :
Direct_Fourier_Reconst;
to run a test example.
Keywords: Tomography, Fourier slice theorem, Tomographic reconstruction
1.0.0.0  new image shows the value of the similarity index 

1.0.0.0  computes Structural Similarity Index for measuring image quality : ssim(reconstructed image, phantom image) 

1.0.0.0  further checks on input variables 

1.0.0.0  updated the "Description" 

1.0.0.0  High frequency artifact in the reconstructed slice are mitigated by zeroing the ’corners’ of the pseudo fft2 

1.0.0.0  Minor bug corrected. 

1.0.0.0  Removed the RadialBasis method of interpolation and replaced with a method that uses 2D filtering (rather than interpolation). Also histcounts is used to allocate the Fourier coefficients in a square grid. 

1.0.0.0  corrected typos 

1.0.0.0  added possibility to filter the projections by editing the code 

1.0.0.0  added "bowtie double wedge" figure : fft2(sinogram) 

1.0.0.0  updated description 

1.0.0.0  added simple alternative that uses griddata 

1.0.0.0  new output figure 

1.0.0.0  improved description 

1.0.0.0  none 

1.0.0.0  added possibility to use on half Fourier space to reconstruct the image 

1.0.0.0  crops down ( to the original phantom size) the image reconstructed by padding the sinogram 

1.0.0.0  minor changes in the output 
Create scripts with code, output, and formatted text in a single executable document.
Jonathan Alis Lima (view profile)
Thank you so much.
First implemetation of central slice theorem I found.