O(N3 log N) backprojection algorithm for the 3-D radon transform

We present a novel backprojection algorithm for three-dimensional (3-D) radon transform data that requires O(N3 log2 N) operations for reconstruction of an N x N x N volume from O(N2) plane-integral projections. Our algorithm uses a hierarchical decomposition of the 3-D radon transform to recursively decompose the backprojection operation. Simulations are presented demonstrating reconstruction quality comparable to the standard filtered backprojection, which requires O(N5) computations under the same circumstances.     

VOIR: a volumetric image reconstruction algorithm based on Fouriertechniques for inversion of the 3-D Radon transform

A novel volumetric image reconstruction algorithm known as VOIR is presented for inversion of the 3-D Radon transform or its radial derivative. The algorithm is a direct implementation of the projection slice theorem for plane integrals. It generalizes one of the most successful methods in 2-D Fourier image reconstruction involving concentric-square rasters to 3-D; in VOIR, the spectral data, which is calculated by fast Fourier techniques, lie on concentric cubes and are interpolated by a bilinear method on the sides of these concentric cubes. The algorithm has great computational advantages over filtered-backprojection algorithms; for images of side dimension N, the numerical complexity of VOIR is O(N(3) log N) instead of O(N (4)) for backprojection techniques. An evaluation of the image processing performance is reported by comparison of reconstructed images from simulated cone-beam scans of a contrast and resolution test object. The image processing performance is also characterized by an analysis of the edge response from the reconstructed images.     

Rotate-and-slant projector for fast LOR-based fully-3-D iterative PET reconstruction

One of the greatest challenges facing iterative fully-3-D positron emission tomography (PET) reconstruction is the issue of long reconstruction times due to the large number of measurements for 3-D mode as compared to 2-D mode. A rotate-and-slant projector has been developed that takes advantage of symmetries in the geometry to compute volumetric projections to multiple oblique sinograms in a computationally efficient manner. It is based upon the 2-D rotation-based projector using the three-pass method of shears, and it conserves the 2-D rotator computations for multiple projections to each oblique sinogram set. The projector is equally applicable to both conventional evenly-spaced projections and unevenly-spaced line-of-response (LOR) data. The LOR-based version models the location and orientation of the individual LORs (i.e., the arc-correction), providing an ordinary Poisson reconstruction framework. The projector was implemented in C with several optimizations for speed, exploiting the vertical symmetry of the oblique projection process, depth compression, and array indexing schemes which maximize serial memory access. The new projector was evaluated and compared to ray-driven and distance-driven projectors using both analytical and experimental phantoms, and fully-3-D iterative reconstructions with each projector were also compared to Fourier rebinning with 2-D iterative reconstruction. In terms of spatial resolution, contrast, and background noise measures, 3-D LOR-based iterative reconstruction with the rotate-and-slant projector performed as well as or better than the other methods. Total processing times, measured on a single cpu Linux workstation, were approximately 10x faster for the rotate-and-slant projector than for the other 3-D projectors studied. The new projector provided four iterations fully-3-D ordered-subsets reconstruction in as little as 15 s--approximately the same time as FORE + 2-D reconstruction. We conclude that the rotate-and-slant projector is a viable option for fully-3-D PET, offering quality statistical reconstruction in times only marginally slower than 2-D or rebinning methods.     

Partial Radon transforms

This article formally defines partial Radon transforms for functions of more than two dimensions. It shows that a generalized projection-slice theorem exists which connects planar and hyperplanar projections of a function to its Fourier transform. In addition, a general theoretical framework is provided for carrying out n-dimensional backprojection reconstruction in a multistage fashion through the use of the partial Radon transform.     

Fast predictions of variance images for fan-beam transmission tomography with quadratic regularization

Accurate predictions of image variances can be useful for reconstruction algorithm analysis and for the design of regularization methods. Computing the predicted variance at every pixel using matrix-based approximations [1] is impractical. Even most recently adopted methods that are based on local discrete Fourier approximations are impractical since they would require a forward and backprojection and two fast Fourier transform (FFT) calculations for every pixel, particularly for shift-variant systems like fan-beam tomography. This paper describes new "analytical" approaches to predicting the approximate variance maps of 2-D images that are reconstructed by penalized-likelihood estimation with quadratic regularization in fan-beam geometries. The simplest of the proposed analytical approaches requires computation equivalent to one backprojection and some summations, so it is computationally practical even for the data sizes in X-ray computed tomography (CT). Simulation results show that it gives accurate predictions of the variance maps. The parallel-beam geometry is a simple special case of the fan-beam analysis. The analysis is also applicable to 2-D positron emission tomography (PET).     

Reconstruction algorithms: Transform methods

Transform methods for image reconstruction from projections are based on analytic inversion formulas. In this tutorial paper, the inversion formula for the case of two-dimensional (2-D) reconstruction from line integrals is manipulated into a number of different forms, each of which may be discretized to obtain different algorithms for reconstruction from sampled data. For the convolution-backprojection algorithm and the direct Fourier algorithm the emphasis is placed on understanding the relationship between the discrete operations specified by the algorithm and the functional operations expressed by the inversion formula. The performance of the Fourier algorithm may be improved, with negligible extra computation, by interleaving two polar sampling grids in Fourier space. The convolution-backprojection formulas are adapted for the fan-beam geometry, and other reconstruction methods are summarized, including the rho-filtered layergram method, and methods involving expansions in angular harmonics. A standard mathematical process leads to a known formula for iterative reconstruction from projections at a finite number of angles. A new iterative reconstruction algorithm is obtained from this formula by introducing one-dimensional (1-D) and 2-D interpolating functions, applied to sampled projections and images, respectively. These interpolating functions are derived by the same Fourier approach which aids in the development and understanding of the more conventional transform methods.     

Iterative image reconstruction using inverse Fourier rebinning for fully 3-D PET

We describe a fast forward and back projector pair based on inverse Fourier rebinning for use in iterative image reconstruction for fully 3-D positron emission tomography (PET). The projector pair is used as part of a factored system matrix that takes into account detector-pair response by using shift-variant sinogram blur kernels, thereby combining the computational advantages of Fourier rebinning with iterative reconstruction using accurate system models. The forward projector consists of a 2-D projector, which maps 3-D images into 2-D direct sinograms, followed by exact inverse rebinning which maps the 2-D into fully 3-D sinograms. The back projector is implemented as the transpose of the forward projector and differs from the true exact rebinning operator in the sense that it does not require reprojection to compute missing lines of response (LORs). We compensate for two types of inaccuracies that arise in a cylindrical PET scanner when using inverse Fourier rebinning: 1) nonuniform radial sampling and 2) nonconstant oblique angles in the radial direction in a single oblique sinogram. We examine the effects of these corrections on sinogram accuracy and reconstructed image quality. We evaluate performance of the new projector pair for maximum a posteriori (MAP) reconstruction of simulated and in vivo data. The new projector results in only a small loss in resolution towards the edge of the field-of-view when compared to the fully 3-D geometric projector and requires an order of magnitude less computation.     

Continuous and discrete data rebinning in time-of-flight PET

This paper investigates data compression methods for time-of-flight (TOF) positron emission tomography (PET), which rebin the 3-D TOF measurements into a set of 2-D TOF data for a stack of transaxial slices. The goal of this work is to develop rebinning algorithms that are more accurate than the TOF single-slice-rebinning (TOF-SSRB) method proposed by Mullani in 1982. Two approaches are explored. The first one is based on a partial differential equation, which expresses a consistency condition for TOF-PET data with a Gaussian TOF profile. From this equation we derive an analytical rebinning algorithm, which is unbiased in the limit of continuous sampling. The second approach is discrete: each 2-D rebinned data sample is calculated as a linear combination of the 3-D TOF samples in the same axial plane parallel to the axis of the scanner. The coefficients of the linear combination are precomputed by optimizing a cost function which enforces both accuracy and good variance reduction, models the TOF profile, the axial PSF of the LORs, and the specific sampling scheme of the scanner. Measurements of a thorax phantom on a prototype TOF-PET scanner with a resolution of 550 ps show that the proposed discrete method improves the bias-variance trade-off and is a promising alternative to TOF-SSRB when data compression is required to achieve clinically acceptable reconstruction time.     

Rectification for cone-beam projection and backprojection

The purpose of this paper is to derive a technique for accelerating the computation of cone-beam forward and backward projections that are the basic steps of tomographic reconstruction. The cone-beam geometry of C-arm systems is commonly described with projection matrices. Such matrices provide a continuous framework for analyzing the flow of operations needed to compute backprojection for analytical reconstruction, as well as the combination of forward and backward projections for iterative reconstruction. The proposed rectification technique resampies the original data to planes that are aligned with two of the reconstructed volume main axes, so that the original cone-beam geometry can be replaced by a simpler geometry, where succession of plane magnifications are involved only. Rectification generalizes previous independent results to the cone-beam backprojection of preprocessed data as well as to cone-beam iterative reconstruction. The memory access pattern of simple magnifications provides superior predictability and is, therefore, easier to optimize, independently of the choice of the interpolation technique. Rectification is also shown to provide control over interpolation errors through oversampling, allowing tradeoffs between computation speed and precision to be set. Experimental results are provided for linear and nearest neighbor interpolations, based on simulations, as well as phantom and patient data acquired on a digital C-arm system.     

Iterative image reconstruction using inverse fourier rebinning for fully 3-D PET

We describe a fast forward and back projector pair based on inverse Fourier rebinning for use in iterative image reconstruction for fully three-dimensional (3-D) positron emission tomography (PET). The projector pair is used as part of a factored system matrix that takes into account detector-pair response by using shift-variant sinogram blur kernels, thereby combining the computational advantages of Fourier rebinning with iterative reconstruction using accurate system models. The forward projector consists of a two-dimensional (2-D) projector, which maps 3-D images into 2-D direct sinograms, followed by exact inverse rebinning which maps the 2-D into fully 3-D sinograms. The back projector is implemented as the transpose of the forward projector and differs from the true exact rebinning operator in the sense that it does not require reprojection to compute missing line of responses (LORs). We compensate for two types of inaccuracies that arise in a cylindrical PET scanner when using inverse Fourier rebinning: 1) nonuniform radial sampling and 2) nonconstant oblique angles in the radial direction in a single oblique sinogram. We examine the effects of these corrections on sinogram accuracy and reconstructed image quality. We evaluate performance of the new projector pair for maximum a posteriori (MAP) reconstruction of simulated and in vivo data. The new projector results in only a small loss in resolution towards the edge of the field-of-view when compared to the fully 3-D geometric projector and requires an order of magnitude less computation.     

Fourier rebinning applied to multiplanar circular-orbit cone-beam SPECT

We study the application of Fourier rebinning methods to dual-planar cone-beam SPECT. Dual-planar cone-beam SPECT involves the use of a pair of dissimilar cone-beam collimators on a dual-camera SPECT system. Each collimator has its focus in a different axial plane. While dual-planar data is best reconstructed with fully three-dimensional (3-D) iterative methods, these methods are slow and have prompted a search for faster reconstruction techniques. Fourier rebinning was developed to estimate equivalent parallel projections from 3-D PET data, but it simply expresses a relationship between oblique projections taken in planes not perpendicular to the axis of rotation and direct projections taken in those that are. We find that it is possible to put cone-beam data in this context as well. The rebinned data can then be reconstructed using either filtered backprojection (FBP) or parallel iterative algorithms such as OS-EM. We compare the Feldkamp algorithm and fully 3-D OSEM reconstruction with Fourier-rebinned reconstructions on realistically-simulated Tc-99m HMPAO brain SPECT data. We find that the Fourier-rebinned reconstructions exhibit much less image noise and lower variance in region-of-interest (ROI) estimates than Feldkamp. Also, Fourier-rebinning followed by OSEM with nonuniform attenuation correction exhibits less bias in ROI estimates than Feldkamp with Chang attenuation correction. The Fourier-rebinned ROI estimates exhibit bias and variance comparable to those from fully 3-D OSEM and require considerably less processing time. However, in areas off the axis of rotation, the axial-direction resolution of FORE-reconstructed images is poorer than that of images reconstructed with 3-D OSEM. We conclude that Fourier rebinning is a practical and potentially useful approach to reconstructing data from dual-planar circular-orbit cone-beam systems.     

Fourier-based reconstruction for fully 3-D PET: optimization of interpolation parameters

Fourier-based approaches for three-dimensional (3-D) reconstruction are based on the relationship between the 3-D Fourier transform (FT) of the volume and the two-dimensional (2-D) FT of a parallel-ray projection of the volume. The critical step in the Fourier-based methods is the estimation of the samples of the 3-D transform of the image from the samples of the 2-D transforms of the projections on the planes through the origin of Fourier space, and vice versa for forward-projection (reprojection). The Fourier-based approaches have the potential for very fast reconstruction, but their straightforward implementation might lead to unsatisfactory results if careful attention is not paid to interpolation and weighting functions. In our previous work, we have investigated optimal interpolation parameters for the Fourier-based forward and back-projectors for iterative image reconstruction. The optimized interpolation kernels were shown to provide excellent quality comparable to the ideal sinc interpolator. This work presents an optimization of interpolation parameters of the 3-D direct Fourier method with Fourier reprojection (3D-FRP) for fully 3-D positron emission tomography (PET) data with incomplete oblique projections. The reprojection step is needed for the estimation (from an initial image) of the missing portions of the oblique data. In the 3D-FRP implementation, we use the gridding interpolation strategy, combined with proper weighting approaches in the transform and image domains. We have found that while the 3-D reprojection step requires similar optimal interpolation parameters as found in our previous studies on Fourier-based iterative approaches, the optimal interpolation parameters for the main 3D-FRP reconstruction stage are quite different. Our experimental results confirm that for the optimal interpolation parameters a very good image accuracy can be achieved even without any extra spectral oversampling, which is a common practice to decrease errors caused by interpolation in Fourier reconstruction.     

Evaluating performance of reconstruction algorithms for 3-D [15O] water PET using subtraction analysis

Positron emission tomography (PET) [15O] activation studies have benefited significantly from three-dimensional (3-D) data acquisition. However, they have been slow to take advantage of new 3-D reconstruction techniques. Compared with the widely used 3-D reprojection reconstruction (3DRP), the advantage of signal and noise for iterative algorithms has been outweighed by concern about long and complicated reconstruction procedures and inconsistent performance. Most pseudo-3-D algorithms, such as rebinning methods, aim at increasing the speed of reconstruction but lack further resolution improvement or noise control. Although many evaluations have been conducted through simulations and phantom experiments, the spatially varying nature of signal and noise and the complexity of biological effects have complicated the interpretation of real data based on simulation or phantom results. We have taken a different approach and used the analysis of real data directly as a measure with which to compare three reconstruction algorithms: 3DRP, iterative filtered backprojection with median root prior (IFBP-MRP), and Fourier rebinning followed by two-dimensional (2-D) filtered backprojection (FORE-FBP) for [15O] PET. Two subjects, each with 32 scans acquired in four sessions during a finger opposition motor task, are analyzed using subtraction. A fixed volume-of-interest (VOI) measurement in regions related to the task demonstrates that at high resolution, IFBP-MRP has the best signal-to-noise performance followed by 3DRP and FORE-FBP; however, this advantage gradually diminishes as the resolution decreases. For a voxel measurement derived from the image of each reconstruction, all three algorithms are capable of detecting highly activated regions. Although there are some differences in the size, shape, and center location of the activated foci, our preliminary results suggest that IFBP-MRP does offer enhanced signal with some noise control compared with 3DRP for the analysis of high-resolution images. If images are to be analyzed at an intermediate to lower resolution, FORE-FBP provides a significant reduction of reconstruction time compared with 3DRP.     

Planogram rebinning with the frequency-distance relationship

We present an efficient rebinning algorithm for positron emission tomography (PET) systems with panel detectors. The rebinning algorithm is derived in the planogram coordinate system which is the native data format for PET systems with panel detectors and is the 3-D extension of the 2-D linogram transform developed by Edholm. Theoretical error bounds and numerical results are included.      

A fast and accurate Fourier algorithm for iterative parallel-beamtomography

We use a series-expansion approach and an operator framework to derive a new, fast, and accurate Fourier algorithm for iterative tomographic reconstruction. This algorithm is applicable for parallel-ray projections collected at a finite number of arbitrary view angles and radially sampled at a rate high enough that aliasing errors are small. The conjugate gradient (CG) algorithm is used to minimize a regularized, spectrally weighted least-squares criterion, and we prove that the main step in each iteration is equivalent to a 2-D discrete convolution, which can be cheaply and exactly implemented via the fast Fourier transform (FFT). The proposed algorithm requires O(N²logN) floating-point operations per iteration to reconstruct an N×N image from P view angles, as compared to O(N ²P) floating-point operations per iteration for iterative convolution-backprojection algorithms or general algebraic algorithms that are based on a matrix formulation of the tomography problem. Numerical examples using simulated data demonstrate the effectiveness of the algorithm for sparse- and limited-angle tomography under realistic sampling scenarios. Although the proposed algorithm cannot explicitly account for noise with nonstationary statistics, additional simulations demonstrate that for low to moderate levels of nonstationary noise, the quality of reconstruction is almost unaffected by assuming that the noise is stationary     

Backprojection by upsampled Fourier series expansion andinterpolated FFT

A fast backprojection method through the use of interpolated fast Fourier transform (FFT) is presented. The computerized tomography (CT) reconstruction by the convolution backprojection (CBP) method has produced precise images. However, the backprojection part of the conventional CBP method is not very efficient. The authors propose an alternative approach to interpolating and backprojecting the convolved projections onto the image frame. First, the upsampled Fourier series expansion of the convolved projection is calculated. Then, using a Gaussian function, it is projected by the aliasing-free interpolation of FFT bins onto a rectangular grid in the frequency domain. The total amount of computation in this procedure for a 512x512 image is 1/5 of the conventional backprojection method with linear interpolation. This technique also allows the arbitrary control of the frequency characteristics.     

Four-dimensional spectral-spatial imaging using projection reconstruction

An n-dimensional (n-D) filtered backprojection image reconstruction algorithm has been developed and used in the reconstruction of 4-D spectral-spatial magnetic resonance imaging (MRI) data. The algorithm uses n-1 successive stages of 2-D filtered backprojection to reconstruct an n-D image. This approach results in a reduction in computational time on the order of N(n-2) relative to the single-stage technique, where N(n) is the number of elements in an n-D image. The authors describe implementation of the algorithm, including digital filtering and sampling requirements. Images obtained from simulated data are presented to illustrate the accuracy and potential utility of the technique.     

The integer transforms analogous to discrete trigonometric transforms

The integer transform (such as the Walsh transform) is the discrete transform that all the entries of the transform matrix are integer. It is much easier to implement because the real number multiplication operations can be avoided, but the performance is usually worse. On the other hand, the noninteger transform, such as the DFT and DCT, has a good performance, but real number multiplication is required. W derive the integer transforms analogous to some popular noninteger transforms. These integer transforms retain most of the performance quality of the original transform, but the implementation is much simpler. Especially, for the two-dimensional (2-D) block transform in image/video, the saving can be huge using integer operations. In 1989, Cham had derived the integer cosine transform. Here, we will derive the integer sine, Hartley, and Fourier transforms. We also introduce the general method to derive the integer transform from some noninteger transform. Besides, the integer transform derived by Cham still requires real number multiplication for the inverse transform. We modify the integer transform introduced by Cham and introduce the complete integer transform. It requires no real number multiplication operation, no matter what the forward or inverse transform. The integer transform we derive would be more efficient than the original transform. For example, for the 8-point DFT and IDFT, there are in total four real numbers and eight fixed-point multiplication operations required, but for the forward and inverse 8-point complete integer Fourier transforms, there are totally 20 fixed-point multiplication operations required. However, for the integer transform, the implementation is simpler, and many of the properties of the original transform are kept.     

An efficient Fourier method for 3-D radon inversion in exactcone-beam CT reconstruction

The radial derivative of the three-dimensional (3-D) radon transform of an object is an important intermediate result in many analytically exact cone-beam reconstruction algorithms. The authors briefly review Grangeat's (1991) approach for calculating radon derivative data from cone-beam projections and then present a new, efficient method for 3-D radon inversion, i.e., reconstruction of the image from the radial derivative of the 3-D radon transform, called direct Fourier inversion (DFI). The method is based directly on the 3-D Fourier slice theorem. From the 3-D radon derivative data, which is assumed to be sampled on a spherical grid, the 3-D Fourier transform of the object is calculated by performing fast Fourier transforms (FFTs) along radial lines in the radon space. Then, an interpolation is performed from the spherical to a Cartesian grid using a 3-D gridding step in the frequency domain. Finally, this 3-D Fourier transform is transformed back to the spatial domain via 3-D inverse FFT. The algorithm is computationally efficient with complexity in the order of N ³ log N. The authors have done reconstructions of simulated 3-D radon derivative data assuming sampling conditions and image quality requirements similar to those in medical computed tomography (CT)     

A high-speed reconstruction from projections using direct Fourier method with optimized parameters-an experimental analysis

Some issues of the direct Fourier method (DFM) implementation are discussed. A hybrid spline-linear interpolation for the DFM is proposed. The results of comprehensive simulation research are presented. The following reconstruction problems and parameters are emphasized: interpolation, increasing the radial density of the polar raster, filtering, the 2-D inverse Fourier transformation dimension, and considering the cases of noiseless and noisy input data. For the a priori prescribed resolution of the reconstructed image, values of reconstruction parameters have been determined which are optimal with regard to reconstruction quality and computation cost. The computational requirements of the DFM algorithm which correspond to distinct interpolation schemes are compared to one another for CT and MR tomography, respectively. The estimations obtained are compared to computational characteristics of the convolution backprojection method.