Robert M. Lewitt

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.

Alternatives to voxels for image representation in iterative reconstruction algorithms

Spherically symmetric volume elements are alternatives to the more conventional voxels for the construction of volume images in the computer. The image representation, and the calculation of projections of it, are essential components of iterative algorithms for image reconstruction from projection data. A two-parameter family of spherical volume elements is described that allows control of the smoothness properties of the represented image, whereas conventional voxels are discontinuous. The rotational symmetry of the spherical elements leads to efficient calculation of projections of the represented image, as required in iterative reconstruction algorithms. For volume elements whose shape is ellipsoidal (rather than spherical) it is shown that efficient calculation of the projections is also possible by means of an image space transformation.

Practical considerations for 3-D image reconstruction using spherically symmetric volume elements

Spherically symmetric volume elements with smooth tapering of the values near their boundaries are alternatives to the more conventional voxels for the construction of volume images in the computer. Their use, instead of voxels, introduces additional parameters which enable the user to control the shape of the volume element (blob) and consequently to control the characteristics of the images produced by iterative methods for reconstruction from projection data. For images composed of blobs, efficient algorithms have been designed for the projection and discrete back-projection operations, which are the crucial parts of iterative reconstruction methods. The authors have investigated the relationship between the values of the blob parameters and the properties of images represented by the blobs. Experiments show that using blobs in iterative reconstruction methods leads to substantial improvement in the reconstruction performance, based on visual quality and on quantitative measures, in comparison with the voxel case. The images reconstructed using appropriately chosen blobs are characterized by less image noise for both noiseless data and noisy data, without loss of image resolution.

Multidimensional digital image representations using generalized Kaiser-Bessel window functions.

Inverse problems that require the solution of integral equations are inherent in a number of indirect imaging applications, such as computerized tomography. Numerical solutions based on discretization of the mathematical model of the imaging process, or on discretization of analytic formulas for iterative inversion of the integral equations, require a discrete representation of an underlying continuous image. This paper describes discrete image representations, in n-dimensional space, that are constructed by the superposition of shifted copies of a rotationally symmetric basis function. The basis function is constructed using a generalization of the Kaiser-Bessel window function of digital signal processing. The generalization of the window function involves going from one dimension to a rotationally symmetric function in n dimensions and going from the zero-order modified Bessel function of the standard window to a function involving the modified Bessel function of order m. Three methods are given for the construction, in n-dimensional space, of basis functions having a specified (finite) number of continuous derivatives, and formulas are derived for the Fourier transform, the x-ray transform, the gradient, and the Laplacian of these basis functions. Properties of the new image representations using these basis functions are discussed, primarily in the context of two-dimensional and three-dimensional image reconstruction from line-integral data by iterative inversion of the x-ray transform. Potential applications to three-dimensional image display are also mentioned.

Accelerated Iterative Reconstruction for Positron Emission Tomography Based on the EM Algorithm for Maximum Likelihood Estimation

The EM method that was originally developed for maximum likelihood estimation in the context of mathematical statistics may be applied to a stochastic model of positron emission tomography (PET). The result is an iterative algorithm for image reconstruction that is finding increasing use in PET, due to its attractive theoretical and practical properties. Its major disadvantage is the large amount of computation that is often required, due to the algorithms slow rate of convergence. This paper presents an accelerated form of the EM algorithm for PET in which the changes to the image, as calculated by the standard algorithm, are multiplied at each iteration by an overrelaxation parameter. The accelerated algorithm retains two of the important practical properties of the standard algorithm, namely the selfnormalization and nonnegativity of the reconstructed images. Experimental results are presented using measured data obtained from a hexagonal detector system for PET. The likelihood function and the norm of the data residual were monitored during the iterative process. According to both of these measures, the images reconstructed at iterations 7 and 11 of the accelerated algorithm are similar to those at iterations 15 and 30 of the standard algorithm, for two different sets of data. Important theoretical properties remain to be investigated, namely the convergence of the accelerated algorithm and its performance as a maximum likelihood estimator.

Fourier correction for spatially variant collimator blurring in SPECT

In single-photon emission computed tomography (SPECT), projection data are acquired by rotating the photon detector around a patient, either in a circular orbit or in a noncircular orbit. The projection data of the desired spatial distribution of emission activity is blurred by the point-response function of the collimator that is used to define the range of directions of gamma-ray photons reaching the detector. The point-response function of the collimator is not spatially stationary, but depends on the distance from the collimator to the point. Conventional methods for deblurring collimator projection data are based on approximating the actual distance-dependent point-response function by a spatially invariant blurring function, so that deconvolution methods can be applied independently to the data at each angle of view. A method is described here for distance-dependent preprocessing of SPECT projection data prior to image reconstruction. Based on the special distance-dependent characteristics of the Fourier coefficients of the sinogram, a spatially variant inverse filter can be developed to process the projection data in all views simultaneously. The algorithm is first derived from Fourier analysis of the projection data from the circular orbit geometry. For circular orbit projection data, experimental results from both simulated data and real phantom data indicate the potential of this method. It is shown that the spatial filtering method can be extended to the projection data from the noncircular orbit geometry. Experiments on simulated projection data from an elliptical orbit demonstrate correction of the spatially variant blurring and distortion in the reconstructed image caused by the noncircular orbit geometry.

Novel Properties Of The Fourier Decomposition Of The Sinogram

The double Fourier decomposition of the sinogram is obtained by first taking the Fourier transform of each parallel-ray projection and then calculating the coefficients of a Fourier series with respect to angle for each frequency component of the transformed projections. The values of these coefficients may be plotted on a two-dimensional map whose coordinates are spatial frequency w (continuous) and angular harmonic number n (discrete). For |w| large, the Fourier coefficients on the line n=kw of slope k through the origin of the coefficient space are found to depend strongly on the contributions to the projection data that, for each view, come from a certain distance to the detector plane, where the distance is a linear function of k. The values of these coefficients depend only weakly on contributions from other distances from the detector. The theoretical basis of this property is presented in this paper and a potential application to emission computerized tomography is discussed.

Three-dimensional image reconstruction for PET by multi-slice rebinning and axial image filtering

A fast method is described for reconstructing volume images from three-dimensional (3D) coincidence data in positron emission tomography (PET). The reconstruction method makes use of all coincidence data acquired by high-sensitivity PET systems that do not have inter-slice absorbers (septa) to restrict the axial acceptance angle. The reconstruction method requires only a small amount of storage and computation, making it well suited for dynamic and whole-body studies. The method consists of three steps: (i) rebinning of coincidence data into a stack of 2D sinograms; (ii) slice-by-slice reconstruction of the sinogram associated with each slice to produce a preliminary 3D image having strong blurring in the axial (z) direction, but with different blurring at different z positions; and (iii) spatially variant filtering of the 3D image in the axial direction (i.e. 1D filtering in z for each x-y column) to produce the final image. The first step involves a new form of the rebinning operation in which multiple sinograms are incremented for each oblique coincidence line (multi-slice rebinning). The axial filtering step is formulated and implemented using the singular value decomposition (SVD). The method has been applied successfully to simulated data and to measured data for different kinds of phantom (multiple point sources, multiple discs, a cylinder with cold spheres, and a 3D brain phantom).

Fourier Method For Correction Of Depth-Dependent Collimator Blurring

A method is described for preprocessing projection data prior to image reconstruction in single-photon emission computed tomography. The projection data of the desired spatial distribution of emission activity is blurred by the point-response function of the collimator that is used to define the range of directions of gamma-ray photons reaching the detector. The point-response function of the collimator is not stationary, but depends on the distance from the collimator to the point. Conventional methods for deblurring collimator projection data are based on approximating the actual depth-dependent point-response function by a spatially-invariant blurring function, so that deconvolution methods can be applied independently to the data at each angle of view. The method described in this paper is based on Fourier analysis of the multi-angular data set as a whole, using special depth-dependent characteristics of the Fourier coefficients to achieve spatially-variant inverse filtering of the data in all views simultaneously. Preliminary results are presented for simulated data with a simple collimator model.

A methodology for testing for statistically significant differences between fully 3D PET reconstruction algorithms

We present a practical methodology for evaluating 3D PET reconstruction methods. It includes generation of random samples from a statistically described ensemble of 3D images resembling those to which PET would be applied in a medical situation, generation of corresponding projection data with noise and detector point spread function simulating those of a 3D PET scanner, assignment of figures of merit appropriate for the intended medical applications, optimization of the reconstruction algorithms on a training set of data, and statistical testing of the validity of hypotheses that say that two reconstruction algorithms perform equally well (from the point of view of a particular figure of merit) as compared to the alternative hypotheses that say that one of the algorithms outperforms the other. Although the methodology was developed with the 3D PET in mind, it can be used, with minor changes, for other 3D data collection methods, such as fully 3D cr or SPECT.