Backprojection Operator Research Articles

A distributed ASTRA toolbox.

While iterative reconstruction algorithms for tomography have several advantages compared to standard backprojection methods, the adoption of such algorithms in large-scale imaging facilities is still limited, one of the key obstacles being their high computational load. Although GPU-enabled computing clusters are, in principle, powerful enough to carry out iterative reconstructions on large datasets in reasonable time, creating efficient distributed algorithms has so far remained a complex task, requiring low-level programming to deal with memory management and network communication. The ASTRA toolbox is a software toolbox that enables rapid development of GPU accelerated tomography algorithms. It contains GPU implementations of forward and backprojection operations for many scanning geometries, as well as a set of algorithms for iterative reconstruction. These algorithms are currently limited to using GPUs in a single workstation. In this paper, we present an extension of the ASTRA toolbox and its Python interface with implementations of forward projection, backprojection and the SIRT algorithm that can be distributed over multiple GPUs and multiple workstations, as well as the tools to write distributed versions of custom reconstruction algorithms, to make processing larger datasets with ASTRA feasible. As a result, algorithms that are implemented in a high-level conceptual script can run seamlessly on GPU-enabled computing clusters, up to 32 GPUs or more. Our approach is not limited to slice-based reconstruction, facilitating a direct portability of algorithms coded for parallel-beam synchrotron tomography to cone-beam laboratory tomography setups without making changes to the reconstruction algorithm.

An iterative method to reconstruct the refractive index of a medium from time-of-flight measurements**The authors dedicate this article to Prof Todd Quinto on the occasion of his 65th birthday.

The article deals with a classical inverse problem: the computation of the refractive index of a medium from ultrasound time-of-flight measurements. This problem is very popular in seismics but also for tomographic problems in inhomogeneous media. For example ultrasound vector field tomography needs a priori knowledge of the sound speed. According to Fermat’s principle ultrasound signals travel along geodesic curves of a Riemannian metric which is associated with the refractive index. The inverse problem thus consists of determining the index of refraction from integrals along geodesics curves associated with the integrand leading to a nonlinear problem. In this article we describe a numerical solver for this problem scheme based on an iterative minimization method for an appropriate Tikhonov functional. The outcome of the method is a stable approximation of the sought index of refraction as well as a corresponding set of geodesic curves. We prove some analytical convergence results for this method and demonstrate its performance by means of several numerical experiments. Another novelty in this article is the explicit representation of the backprojection operator for the ray transform in Riemannian geometry and its numerical realization relying on a corresponding phase function that is determined by the metric. This gives a natural extension of the conventional backprojection from 2D computerized tomography to inhomogeneous geometries.

SU‐D‐206‐02: Evaluation of Partial Storage of the System Matrix for Cone Beam Computed Tomography Using a GPU Platform

Purpose:Iterative reconstruction algorithms in computed tomography (CT) require a fast method for computing the intersections between the photons’ trajectories and the object, also called ray‐tracing or system matrix computation. This work evaluates different ways to store the system matrix, aiming to reconstruct dense image grids in reasonable time.Methods:We propose an optimized implementation of the Siddon's algorithm using graphics processing units (GPUs) with a novel data storage scheme. The algorithm computes a part of the system matrix on demand, typically, for one projection angle. The proposed method was enhanced with accelerating options: storage of larger subsets of the system matrix, systematic reuse of data via geometric symmetries, an arithmetic‐rich parallel code and code configuration via machine learning. It was tested on geometries mimicking a cone beam CT acquisition of a human head. To realistically assess the execution time, the ray‐tracing routines were integrated into a regularized Poisson‐based reconstruction algorithm. The proposed scheme was also compared to a different approach, where the system matrix is fully pre‐computed and loaded at reconstruction time.Results:Fast ray‐tracing of realistic acquisition geometries, which often lack spatial symmetry properties, was enabled via the proposed method. Ray‐tracing interleaved with projection and backprojection operations required significant additional time. In most cases, ray‐tracing was shown to use about 66 % of the total reconstruction time. In absolute terms, tracing times varied from 3.6 s to 7.5 min, depending on the problem size. The presence of geometrical symmetries allowed for non‐negligible ray‐tracing and reconstruction time reduction. Arithmetic‐rich parallel code and machine learning permitted a modest reconstruction time reduction, in the order of 1 %.Conclusion:Partial system matrix storage permitted the reconstruction of higher 3D image grid sizes and larger projection datasets at the cost of additional time, when compared to the fully pre‐computed approach.This work was supported in part by the Fonds de recherche du Quebec ‐ Nature et technologies (FRQ‐NT). The authors acknowledge partial support by the CREATE Medical Physics Research Training Network grant of the Natural Sciences and Engineering Research Council of Canada (Grant No. 432290).

TH‐EF‐BRA‐03: Assessment of Data‐Driven Respiratory Motion‐Compensation Methods for 4D‐CBCT Image Registration and Reconstruction Using Clinical Datasets

Purpose:To evaluate the performance of a 4D‐CBCT registration and reconstruction method that corrects for respiratory motion and enhances image quality under clinically relevant conditions.Methods:Building on previous work, which tested feasibility of a motion‐compensation workflow using image datasets superior to clinical acquisitions, this study assesses workflow performance under clinical conditions in terms of image quality improvement. Evaluated workflows utilized a combination of groupwise deformable image registration (DIR) and image reconstruction. Four‐dimensional cone beam CT (4D‐CBCT) FDK reconstructions were registered to either mean or respiratory phase reference frame images to model respiratory motion. The resulting 4D transformation was used to deform projection data during the FDK backprojection operation to create a motion‐compensated reconstruction. To simulate clinically realistic conditions, superior quality projection datasets were sampled using a phase‐binned striding method. Tissue interface sharpness (TIS) was defined as the slope of a sigmoid curve fit to the lung‐diaphragm boundary or to the carina tissue‐airway boundary when no diaphragm was discernable. Image quality improvement was assessed in 19 clinical cases by evaluating mitigation of view‐aliasing artifacts, tissue interface sharpness recovery, and noise reduction.Results:For clinical datasets, evaluated average TIS recovery relative to base 4D‐CBCT reconstructions was observed to be 87% using fixed‐frame registration alone; 87% using fixed‐frame with motion‐compensated reconstruction; 92% using mean‐frame registration alone; and 90% using mean‐frame with motion‐compensated reconstruction. Soft tissue noise was reduced on average by 43% and 44% for the fixed‐frame registration and registration with motion‐compensation methods, respectively, and by 40% and 42% for the corresponding mean‐frame methods. Considerable reductions in view aliasing artifacts were observed for each method.Conclusion:Data‐driven groupwise registration and motion‐compensated reconstruction have the potential to improve the quality of 4D‐CBCT images acquired under clinical conditions. For clinical image datasets, the addition of motion compensation after groupwise registration visibly reduced artifact impact.This work was supported by the National Cancer Institute of the National Institutes of Health under Award Number R01CA166119. Hugo and Weiss hold a research agreement with Philips Healthcare and license agreement with Varian Medical Systems. Weiss receives royalties from UpToDate. Christensen receives funds from Roger Koch to support research.

Reconstruction Algorithms for a Class of Restricted Ray Transforms Without Added Singularities

Let X and \(X^*\) denote a restricted ray transform along curves and a corresponding backprojection operator, respectively. Theoretical analysis of reconstruction from the data Xf is usually based on a study of the composition \(X^* D X\), where D is some local operator (usually a derivative). If \(X^*\) is chosen appropriately, then \(X^* D X\) is a Fourier integral operator (FIO) with singular symbol. The singularity of the symbol leads to the appearance of artifacts (added singularities) that can be as strong as the original (or, useful) singularities. By choosing D in a special way one can reduce the strength of added singularities, but it is impossible to get rid of them completely. In the paper we follow a similar approach, but make two changes. First, we replace D with a nonlocal operator \(\tilde{D}\) that integrates Xf along a curve in the data space. The result \(\tilde{D} Xf\) resembles the generalized Radon transform R of f. The function \(\tilde{D} Xf\) is defined on pairs \((x_0,\Theta )\in U\times S^2\), where \(U\subset {\mathbb R}^3\) is an open set containing the support of f, and \(S^2\) is the unit sphere in \({\mathbb R}^3\). Second, we replace \(X^*\) with a backprojection operator \(R^*\) that integrates with respect to \(\Theta \) over \(S^2\). It turns out that if \(\tilde{D}\) and \(R^*\) are appropriately selected, then the composition \(R^* \tilde{D} X\) is an elliptic pseudodifferential operator of order zero with principal symbol 1. Thus, we obtain an approximate reconstruction formula that recovers all the singularities correctly and does not produce artifacts. The advantage of our approach is that by inserting \(\tilde{D}\) we get access to the frequency variable \(\Theta \). In particular, we can incorporate suitable cut-offs in \(R^*\) to eliminate bad directions \(\Theta \), which lead to added singularities.

GPU-accelerated iterative reconstruction from Compton scattered data using a matched pair of conic projector and backprojector

Background and ObjectiveIterative reconstruction from Compton scattered data is known to be computationally more challenging than that from conventional line-projection based emission data in that the gamma rays that undergo Compton scattering are modeled as conic projections rather than line projections. In conventional tomographic reconstruction, to parallelize the projection and backprojection operations using the graphics processing unit (GPU), approximated methods that use an unmatched pair of ray-tracing forward projector and voxel-driven backprojector have been widely used. In this work, we propose a new GPU-accelerated method for Compton camera reconstruction which is more accurate by using exactly matched pair of projector and backprojector. MethodsTo calculate conic forward projection, we first sample the cone surface into conic rays and accumulate the intersecting chord lengths of the conic rays passing through voxels using a fast ray-tracing method (RTM). For conic backprojection, to obtain the true adjoint of the conic forward projection, while retaining the computational efficiency of the GPU, we use a voxel-driven RTM which is essentially the same as the standard RTM used for the conic forward projector. ResultsOur simulation results show that, while the new method is about 3 times slower than the approximated method, it is still about 16 times faster than the CPU-based method without any loss of accuracy. ConclusionsThe net conclusion is that our proposed method is guaranteed to retain the reconstruction accuracy regardless of the number of iterations by providing a perfectly matched projector–backprojector pair, which makes iterative reconstruction methods for Compton imaging faster and more accurate.

View-dependent geometric calibration for offset flat-panel cone beam computed tomography systems

Geometric parameters that define the geometry of imaging systems are crucial for image reconstruction and image quality in x-ray computed tomography (CT). The problem of determining geometric parameters for an offset flat-panel cone beam CT (CBCT) system, a recently introduced modality with a large field of view, with the assumption of an unstable mechanism and geometric parameters that vary in each view, is considered. To accurately and rapidly find the geometric parameters for each projection view, we use the projection matrix method and design a dedicated phantom that is partially visible in all projection views. The phantom consists of balls distributed symmetrically in a cylinder to ensure the inclusion of the phantom in all views, and a large portion of the phantom is covered in the projection image. To efficiently use calibrated geometric information in the reconstruction process and get rid of approximation errors, instead of decomposing the projection matrix into actual geometric parameters that are manually corrected before being used in reconstruction, as in conventional methods, we directly use the projection matrix and its pseudo-inverse in projection and backprojection operations of reconstruction algorithms. The experiments illustrate the efficacy of the proposed method with a real offset flat-panel CBCT system in dental imaging.

Single remote sensing image scale-up combining modulation transform function compensation

Remote sensing images usually need scale-up for visualization or representation, using only one original image. According to the performance of detective sensors, a new and more applicable method is proposed here. To enhance the high-frequency components, the modulation transform function compensation (MTFC) part focuses on how to adjust the spatial response before and after launch, under signal-to-noise ratio control. This largely reduces the ring phenomenon from incorrect point spread function guesses. Then a contour stencil prior manages to limit edge artifacts in the upscaled image after MTFC. An iterative backprojection operation with fast convergence is also utilized to bring about intensity and contour consistency. We finally present our analysis based on real images with parallel design for full speed. Compared with existing algorithms, the operator illustrates its potential to keep geometric features and extend the visual and quantitative quality for further analysis.

Fast Algorithms and Efficient GPU Implementations for the Radon Transform and the Back-Projection Operator Represented as Convolution Operators

The Radon transform and its adjoint, the back-projection operator, can both be expressed as convolutions in log-polar coordinates. Hence, fast algorithms for the application of these operators can be constructed by using the FFT, if data is resampled at log-polar coordinates. Radon data is typically measured on an equally spaced grid in polar coordinates, and reconstructions are represented (as images) in Cartesian coordinates. Therefore, in addition to FFT, several steps of interpolation have to be conducted in order to apply the Radon transform and the back-projection operator by means of convolutions. However, in comparison to the interpolation conducted in Fourier-based gridding methods, the interpolation performed in the Radon and image domains will typically deal with functions that are substantially less oscillatory. Reasonable reconstruction results can thus be expected using interpolation schemes of moderate order. The approach also provides better control over the artifacts that can appear due t...

Inversion of the Attenuated Geodesic X-Ray Transform over Functions and Vector Fields on Simple Surfaces

We derive explicit reconstruction formulas for the attenuated geodesic X-ray transform over functions and, in the case of nonvanishing attenuation, vector fields, on a class of simple Riemannian surfaces with boundary. These formulas partly rely on new explicit approaches to construct continuous right-inverses for backprojection operators (and, in turn, holomorphic integrating factors), which were previously unavailable in a systematic form. The reconstruction of functions is presented in two ways, the latter being motivated by numerical considerations and successfully implemented in the last section of the paper (before the appendix). Constructing the aforementioned right-inverses requires that certain Fredholm equations, first appearing in [L. Pestov and G. Uhlmann, Int. Math. Res. Not., 2004 (2004), pp. 4331--4347], be invertible. Whether this last condition reduces the applicability of the overall approach to a strict subset of simple surfaces remains open at present.

On the computational implementation of forward and back-projection operations for cone-beam computed tomography.

Forward- and back-projection operations are the main computational burden in iterative image reconstruction in computed tomography. In addition, their implementation has to be accurate to ensure stable convergence to a high-quality image. This paper reviews and compares some of the variations in the implementation of these operations in cone-beam computed tomography. We compare four algorithms for computing the system matrix, including a distance-driven algorithm, an algorithm based on cubic basis functions, another based on spherically symmetric basis functions, and a voxel-driven algorithm. The focus of our study is on understanding how the choice of the implementation of the system matrix will influence the performance of iterative image reconstruction algorithms, including such factors as the noise strength and spatial resolution in the reconstructed image. Our experiments with simulated and real cone-beam data reveal the significance of the speed-accuracy trade-off in the implementation of the system matrix. Our results suggest that fast convergence of iterative image reconstruction methods requires accurate implementation of forward- and back-projection operations, involving a direct estimation of the convolution of the footprint of the voxel basis function with the surface of the detectors. The required accuracy decreases by increasing the resolution of the projection measurements beyond the resolution of the reconstructed image. Moreover, reconstruction of low-contrast objects needs more accurate implementation of these operations. Our results also show that, compared with regularized reconstruction methods, the behavior of iterative reconstruction algorithms that do not use a proper regularization is influenced more significantly by the implementation of the forward- and back-projection operations.

Reconstruction of the singular support of a tensor field given in a refracting medium by its ray transform

We propose some approaches for numerically solving the problem of reconstructing the singular support of a symmetric tensor field, given in a refracting medium, by its known ray transform. To solve the problem, we use the back-projection operators that act on the ray transforms and the tensor analysis methods on Riemannian manifolds. We construct the operators of the medium inhomogeneity indicator that allow us to identify the sets of points of the singular support of the scalar, vector, and tensor fields. We propose and implement algorithms for solving the problem under study.

Matched Backprojection Operator for Combined Scanning Transmission Electron Microscopy Tilt- and Focal Series.

Combined tilt- and focal series scanning transmission electron microscopy is a recently developed method to obtain nanoscale three-dimensional (3D) information of thin specimens. In this study, we formulate the forward projection in this acquisition scheme as a linear operator and prove that it is a generalization of the Ray transform for parallel illumination. We analytically derive the corresponding backprojection operator as the adjoint of the forward projection. We further demonstrate that the matched backprojection operator drastically improves the convergence rate of iterative 3D reconstruction compared to the case where a backprojection based on heuristic weighting is used. In addition, we show that the 3D reconstruction is of better quality.

TH-CD-303-07: Evaluation of Four Data-Driven Respiratory Motion-Compensation Methods for Four-Dimensional Cone-Beam CT Registration and Reconstruction

Purpose: To develop and compare four 4D-CBCT reconstruction methods that correct for respiratory motion and enhance image quality. Methods: Four motion-compensation workflows were developed employing a combination of deformable image registration (DIR) and image reconstruction. Groupwise registration was used to simultaneously register all frames of a base Four-Dimensional Cone-Beam CT (4D-CBCT) FDK reconstruction to a reference frame, providing a 4D transformation. This 4D transformation expresses a respiratory motion model and was used to apply a deformation during the backprojection operation to create a motion-compensated reconstruction. Two variables were assessed: reference image of registration (fixed or mean frame) and inclusion or exclusion of the motion-compensated reconstruction step. Fixed frame registration refers to using an actual frame from the 4D image (end of inhalation) as the target image. Mean frame registration refers to using a per iteration-estimated average frame image instead. Image quality was assessed in one clinical case through tissue boundary sharpness, noise, and presence of view-aliasing artifact. Results: Diaphragm edge sharpness (DES) was defined as the slope of the intensity gradient along the lung-diaphragm boundary. DES relative to the base 4D-CBCT reconstruction improved by 26.8% using fixed-frame registration alone; 15.5% using fixed-frame with motion-compensated reconstruction; decreased by 2.9% using mean-frame registration alone; and improved by 12.2% using mean-frame with motion-compensated reconstruction. Noise was reduced in soft tissue by 8.7% and 75.8% for the fixed-frame registration and registration with motion-compensation methods, respectively, and by 8.8% and 77.7% for the corresponding mean-frame methods. Reductions in undersampling artifacts were visible for all four methods relative to the base reconstruction while retaining higher frequency spatial details (i.e. blood vessels). Conclusion: A data-driven approach combining groupwise registration and motion-compensation has the potential to improve 4D-CBCT image quality. Adding motion compensation after groupwise registration reduced sharpness somewhat but resulted in reduced noise and structural artifacts. This work was supported by the National Cancer Institute of the National Institutes of Health under Award Number R01CA166119. The authors have no conflicts of interest.

Time-of-flight PET image reconstruction using origin ensembles

The origin ensemble (OE) algorithm is a novel statistical method for minimum-mean-square-error (MMSE) reconstruction of emission tomography data. This method allows one to perform reconstruction entirely in the image domain, i.e. without the use of forward and backprojection operations. We have investigated the OE algorithm in the context of list-mode (LM) time-of-flight (TOF) PET reconstruction. In this paper, we provide a general introduction to MMSE reconstruction, and a statistically rigorous derivation of the OE algorithm. We show how to efficiently incorporate TOF information into the reconstruction process, and how to correct for random coincidences and scattered events. To examine the feasibility of LM-TOF MMSE reconstruction with the OE algorithm, we applied MMSE-OE and standard maximum-likelihood expectation-maximization (ML-EM) reconstruction to LM-TOF phantom data with a count number typically registered in clinical PET examinations. We analyzed the convergence behavior of the OE algorithm, and compared reconstruction time and image quality to that of the EM algorithm. In summary, during the reconstruction process, MMSE-OE contrast recovery (CRV) remained approximately the same, while background variability (BV) gradually decreased with an increasing number of OE iterations. The final MMSE-OE images exhibited lower BV and a slightly lower CRV than the corresponding ML-EM images. The reconstruction time of the OE algorithm was approximately 1.3 times longer. At the same time, the OE algorithm can inherently provide a comprehensive statistical characterization of the acquired data. This characterization can be utilized for further data processing, e.g. in kinetic analysis and image registration, making the OE algorithm a promising approach in a variety of applications.

Parallelizing a Matched Pair of Ray-Tracing Projector and Backprojector for Iterative Cone-Beam CT Reconstruction

Iterative reconstruction methods in X-ray CT can provide better image quality than analytical methods but their applications in practice are still limited due to computationally expensive calculations of repeated projection and backprojection operations. In the past decade, GPU-accelerated methods have been successfully used to reduce the computation time for projection and backprojection. However, it has been of a difficult problem to overcome a trade-off between the accuracy of reconstructed images and the efficiency of parallel computations. For example, when the size of the voxel in the reconstructed volume is larger than that of the detector bin, the use of the conventional unmatched projector-backprojector pair can lower the accuracy of reconstructed images due to the error caused by the mismatch between the projector and backprojector. In this paper, we propose a new GPU-accelerated scheme for the most widely used ray-tracing method (RTM) to perform projection and backprojection operations. Unlike the previous works that accelerate the computation of backprojection by using approximations, our method does not use any approximations for parallelizing the projection and backprojection operations. Since our method is exact, the results are as accurate as those obtained from the nonaccelerated method. We apply our method to iterative reconstruction for dental cone-beam CT systems and test its performance using both the simulated data using a 3-D digital phantom and the real data acquired from an offset flat-panel X-ray CT system. Our experimental results show that, the proposed method achieves a substantially high acceleration rate while retaining the accuracy of the RTM for both projection and backprojection.

Fast parallel algorithm for three-dimensional distance-driven model in iterative computed tomography reconstruction**Projected supported by the National High Technology Research and Development Program of China (Grant No. 2012AA011603) and the National Natural Science Foundation of China (Grant No. 61372172).

The projection matrix model is used to describe the physical relationship between reconstructed object and projection. Such a model has a strong influence on projection and backprojection, two vital operations in iterative computed tomographic reconstruction. The distance-driven model (DDM) is a state-of-the-art technology that simulates forward and back projections. This model has a low computational complexity and a relatively high spatial resolution; however, it includes only a few methods in a parallel operation with a matched model scheme. This study introduces a fast and parallelizable algorithm to improve the traditional DDM for computing the parallel projection and backprojection operations. Our proposed model has been implemented on a GPU (graphic processing unit) platform and has achieved satisfactory computational efficiency with no approximation. The runtime for the projection and backprojection operations with our model is approximately 4.5 s and 10.5 s per loop, respectively, with an image size of 256×256×256 and 360 projections with a size of 512×512. We compare several general algorithms that have been proposed for maximizing GPU efficiency by using the unmatched projection/backprojection models in a parallel computation. The imaging resolution is not sacrificed and remains accurate during computed tomographic reconstruction.

Approximate reconstruction from ray integrals of a function on a domain with low refraction

An approach is proposed for reconstruction from ray integrals of a function defined on a Riemannian domain with low refraction. Using the back-projection operator and the fast Fourier transform, an inversion algorithm is constructed for the ray transform and its numerical study is carried out.

Proton computed tomography reconstruction using a backprojection-then-filtering approach

A novel approach to proton CT reconstruction using backprojection-then-filtering (BPF) is proposed. A list-mode algorithm is formulated accommodating non-linear proton paths. The analytical form is derived for the deblurring kernel necessary for the filtering step. Further, a finite matrix correction is derived to correct for the limited size of the backprojection matrix. High quantitative accuracy in relative stopping power is demonstrated (⩽0.1%) using Monte Carlo simulations. This accuracy makes the algorithm a promising candidate for future proton CT systems in proton therapy applications. For the purposes of reconstruction, each proton path in the object-of-interest was estimated based on a cubic spline fit to the proton entry and exit vectors. The superior spatial-resolution of the BPF method over the standard filtering-then-backprojection approach is demonstrated. As the BPF algorithm requires only one backprojection and filtering operation on a scan data set, it also offers computational advantages over an iterative reconstruction approach.

Inversion without Explicit Jacobian Calculations in Electrical Impedance Tomography

Electrical impedance tomography (EIT) is the inverse problem of finding the internal conductivity distribution of a medium given boundary electrical measurements performed via an electrode array onto its surface. Conventional inversion schemes adopt Tikhonov regularized Newton-type methods. Following a transport back-transport approach, we develop in this work an adjoint approach which allows reducing computational burden in enabling inversion without explicit Jacobian calculation. Forward and back-projection operators are defined from potential gradients, along with their efficient implementation. These derivations allow the transparent use of inversion algorithms. We first check the implementation of operators. We then evaluate how reconstructions perform on simulated noisy data using a preconditioned conjugate gradient. We eventually practice our inversion framework on experimental data acquired in vitro from a saline phantom.