Direct Fourier Reconstruction of a Tomographic Slice

버전 1.0.0.0 (50.8 KB) 작성자: Gianni Schena
Experiments of reconstruction using Fourier Slice Theorem (rather than filtered back projection).
다운로드 수: 1.2K
업데이트 날짜: 2018/6/15

라이선스 보기

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

인용 양식

Gianni Schena (2024). 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. 검색됨 .

MATLAB 릴리스 호환 정보
개발 환경: R2016a
모든 릴리스와 호환
플랫폼 호환성
Windows macOS Linux
카테고리
Help CenterMATLAB Answers에서 Morphological Operations에 대해 자세히 알아보기
커뮤니티
 Power Electronics Control 커뮤니티에 더 많은 파일이 있습니다

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
버전 게시됨 릴리스 정보
1.0.0.0

updated the "Description"
crops down ( to the original phantom size) the image reconstructed by padding the sinogram
improved description
none
corrected typos
Minor bug corrected.
added simple alternative that uses griddata
added possibility to filter the projections by editing the code

new image shows the value of the similarity index
added possibility to use on half Fourier space to reconstruct the image
computes Structural Similarity Index for measuring image quality : ssim(reconstructed image, phantom image)
new output figure
minor changes in the output
Removed the Radial-Basis 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.
further checks on input variables
added "bow-tie double wedge" figure : fft2(sinogram)
updated description
High frequency artifact in the reconstructed slice are mitigated by zeroing the ’corners’ of the pseudo fft2