Evaluation of accelerated iterative x-ray CT image reconstruction using floating point graphics hardware.

Statistical reconstruction methods offer possibilities to improve image quality as compared with analytical methods, but current reconstruction times prohibit routine application in clinical and micro-CT. In particular, for cone-beam x-ray CT, the use of graphics hardware has been proposed to accelerate the forward and back-projection operations, in order to reduce reconstruction times. In the past, wide application of this texture hardware mapping approach was hampered owing to limited intrinsic accuracy. Recently, however, floating point precision has become available in the latest generation commodity graphics cards. In this paper, we utilize this feature to construct a graphics hardware accelerated version of the ordered subset convex reconstruction algorithm. The aims of this paper are (i) to study the impact of using graphics hardware acceleration for statistical reconstruction on the reconstructed image accuracy and (ii) to measure the speed increase one can obtain by using graphics hardware acceleration. We compare the unaccelerated algorithm with the graphics hardware accelerated version, and for the latter we consider two different interpolation techniques. A simulation study of a micro-CT scanner with a mathematical phantom shows that at almost preserved reconstructed image accuracy, speed-ups of a factor 40 to 222 can be achieved, compared with the unaccelerated algorithm, and depending on the phantom and detector sizes. Reconstruction from physical phantom data reconfirms the usability of the accelerated algorithm for practical cases.

Parallel statistical image reconstruction for cone-beam x-ray CT on a shared memory computation platform.

Statistical reconstruction methods offer possibilities of improving image quality as compared to analytical methods, but current reconstruction times prohibit routine clinical applications. To reduce reconstruction times we have parallelized a statistical reconstruction algorithm for cone-beam x-ray CT, the ordered subset convex algorithm (OSC), and evaluated it on a shared memory computer. Two different parallelization strategies were developed: one that employs parallelism by computing the work for all projections within a subset in parallel, and one that divides the total volume into parts and processes the work for each sub-volume in parallel. Both methods are used to reconstruct a three-dimensional mathematical phantom on two different grid densities. The reconstructed images are binary identical to the result of the serial (non-parallelized) algorithm. The speed-up factor equals approximately 30 when using 32 to 40 processors, and scales almost linearly with the number of cpus for both methods. The huge reduction in computation time allows us to apply statistical reconstruction to clinically relevant studies for the first time.

Statistical image reconstruction for transmission tomography using relaxed ordered subset algorithms.

Statistical reconstruction methods offer possibilities for improving image quality as compared to analytical methods, but current reconstruction times prohibit routine clinical applications in x-ray computed tomography (CT). To reduce reconstruction times, we have applied (under) relaxation to ordered subset algorithms. This enables us to use subsets consisting of only single projection angle, effectively increasing the number of image updates within an entire iteration. A second advantage of applying relaxation is that it can help improve convergence by removing the limit cycle behaviour of ordered subset algorithms, which normally do not converge to an optimal solution but rather a suboptimal limit cycle consisting of as many points as there are subsets. Relaxation suppresses the limit cycle behaviour by decreasing the stepsize for approaching the solution. A simulation study for a 2D mathematical phantom and three different ordered subset algorithms shows that all three algorithms benefit from relaxation: equal noise-to-resolution trade-off can be achieved using fewer iterations than the conventional algorithms, while a lower minimal normalized mean square error (NMSE) clearly indicates a better convergence. Two different schemes for setting the relaxation parameter are studied, and both schemes yield approximately the same minimal NMSE.

Evaluation of the ordered subset convex algorithm for cone-beam CT.

Statistical methods for image reconstruction such as maximum likelihood expectation maximization (ML-EM) are more robust and flexible than analytical inversion methods and allow for accurate modelling of the photon transport and noise. Statistical reconstruction is prohibitively slow when applied to clinical x-ray cone-beam CT due to the large data sets and the high number of iterations required for reconstructing high resolution images. One way to reduce the reconstruction time is to use ordered subsets of projections during the iterations, which has been successfully applied to fan-beam x-ray CT. In this paper, we quantitatively analyse the use of ordered subsets in concert with the convex algorithm for cone-beam x-ray CT reconstruction, for the case of circular acquisition orbits. We focus on the reconstructed image accuracy of a 3D head phantom. Acceleration factors larger than 300 were obtained with errors smaller than 1%, with the preservation of signal-to-noise ratio. Pushing the acceleration factor towards 600 by using an increasing number of subsets increases the reconstruction error up to 5% and significantly increases noise. The results indicate that the use of ordered subsets can be extremely useful for cone-beam x-ray CT.

One-step finite-difference time-domain algorithm to solve the Maxwell equations.

We present a one-step algorithm to solve the time-dependent Maxwell equations for systems with spatially varying permittivity and permeability. We compare the results of this algorithm with those obtained from the Yee algorithm and from unconditionally stable algorithms. We demonstrate that for a range of applications the one-step algorithm may be orders of magnitude more efficient than multiple time-step, finite-difference time-domain algorithms. We discuss both the virtues and limitations of this one-step approach.

Higher-order unconditionally stable algorithms to solve the time-dependent Maxwell equations.

For the recently introduced algorithms to solve the time-dependent Maxwell equations [J. S. Kole, M. T. Figge, and H. De Raedt, Phys. Rev. E 64, 066705 (2001)], we construct a variable grid implementation and an improved spatial discretization implementation that preserve the exceptional property of the algorithms to be unconditionally stable by construction. We find that the performance and accuracy of the corresponding algorithms are significant and illustrate their practical relevance by simulating various physical model systems.

Unconditionally stable algorithms to solve the time-dependent Maxwell equations.

Based on the Suzuki product-formula approach, we construct a family of unconditionally stable algorithms to solve the time-dependent Maxwell equations. We describe a practical implementation of these algorithms for one-, two-, and three-dimensional systems with spatially varying permittivity and permeability. The salient features of the algorithms are illustrated by computing selected eigenmodes and the full density of states of one-, two-, and three-dimensional models and by simulating the propagation of light in slabs of photonic band-gap materials.

Quantum Monte Carlo method for attractive Coulomb potentials.

Starting from an exact lower bound on the imaginary-time propagator, we present a path-integral quantum Monte Carlo method that can handle singular attractive potentials. We illustrate the basic ideas of this quantum Monte Carlo algorithm by simulating the ground state of hydrogen and helium.

Morphological image analysis of quantum motion in billiards.

Morphological image analysis is applied to the time evolution of the probability distribution of a quantum particle moving in two- and three-dimensional billiards. It is shown that the time-averaged Euler characteristic of the probability distribution provides a well defined quantity to distinguish between classically integrable and nonintegrable billiards. In three dimensions the time-averaged mean breadth of the probability distribution may also be used for this purpose.