Performance evaluation of an iterative image reconstruction algorithm for positron emission tomography

A new approach to overcome the inconsistency between SPECT and the anatomical map in maximum A-posterior expectation-maximization reconstruction algorithm

Noise reduction while preserving spatial resolution is one of the most important challenges in the reconstructing of emission tomography images. One of the resolving methods is the Bowsher maximum a-posteriori expectation-maximization reconstruction (MAPEM) algorithm. This method considers a binary selection of the neighbors of each voxel based on the prior anatomical values to use in the regularization function. This method is particularly susceptible to imposing the wrong data into the reconstructed image due to the spatial or functional inconsistencies between the anatomical image and the actual activity distribution. Because of the poor spatial resolution of single-photon emission tomography (SPECT) images and the different nature of emission and anatomical imaging, there is not enough certainty of inconsistency with anatomical images. Therefore, we proposed a new weighted Bowsher method that can overcome this weakness while the image quality indexes, especially the spatial resolution, are almost preserved. In the proposed method, each of the neighbors of a specific voxel takes a constant weight depending on the order of its value and independent of its intensity quantity. The proposed method was evaluated using some different physical phantoms and a patient scan. The results show that the proposed method has superiority in the presence of inconsistency; moreover, the proposed method gives nearly similar results to the regular Bowsher MAPEM in case of consistency. In conclusion, we show that using a suitable constant weighting factor in Bowsher MAPEM, one can operatively reduce the image noise while preserving the image quality parameters where the emission tomography images are either consistent or inconsistent with the prior anatomical map.

Image reconstruction with the Heaviside equation in photoacoustic tomography accounting for dispersive acoustic media

A challenging issue in photoacoustic biomedical imaging is to take into account the presence of dispersive acoustic media, since these are prone to induce amplitude attenuation and scattering of the photoacoustic frequency components. These perturbations are largely the cause for which the photoacoustic tomographic image reconstruction from projections lacks a plane-wave transport formalism. Attending this problem, we further develop an analytic formalism of the transport and its numerical implementation accounting for dispersive acoustic media. We differentiate three variations of an acoustically perturbing media. Our object of interest is a numerical description of the light absorption map of a coronal human breast image. Then, we analyze conditions for which the propagation of photoacoustic perturbations can obey the generalized Heaviside telegraph equation. In addition, we provide a study of the causality consistency of the wave propagation models. We observe transport implications due to the presence of dispersive acoustic media and derive model adjustments that include attenuation and diffusion approximations within the two-dimensional forward problem. Next, we restore the inverse problem description with the deduced perturbation components. Finally, we solve the nonlinear inverse problem with a numerical strategy for a filtered backprojection reconstruction. At a stage prior to the image reconstruction, we compensate for the effect of acoustic attenuation and diffusion to calculate the inversions of the wave perturbations located within the projections. In this way, we manage to significantly reduce reconstruction artifacts. In consequence, we prevent the use of some additional image processing of noise reduction. We demonstrate a feasible strategy on how to solve the stated nonlinear inverse problem of photoacoustic tomography accounting for dispersive acoustic media. In particular, we emphasize efforts to achieve an analytical description, and thus an algorithm is placed, for imaged sound perturbations to be cleaned from acoustic scattering in a simplified manner.

Image quality of Zr-89 PET imaging in the Siemens microPET Focus 220 preclinical scanner

PURPOSE

Zr-89 positron emission tomography (PET) is a valuable tool for understanding the biodistribution and pharmacokinetics of antibody-based therapeutics. We compared the image quality of Zr-89 PET and F-18 PET in the Siemens microPET Focus 220 preclinical scanner using different reconstruction methods.

PROCEDURES

Image quality metrics were measured in various Zr-89 and F-18 PET phantoms, including the NEMA NU 4-2008 image quality phantom. Images were reconstructed using various algorithms.

RESULTS

Zr-89 PET had greater image noise, inferior spatial resolution, and greater spillover than F-18 PET, but comparable recovery coefficients for cylinders of various diameters. Of the reconstruction methods, OSEM3D resulted in the lowest noise, highest recovery coefficients, best spatial resolution, but also the greatest spillover. Scatter correction results were found to be sensitive to varying object sizes.

CONCLUSIONS

Zr-89 PET image quality was inferior to that of F-18, and no single reconstruction method was superior in all aspects of image quality.

A Bayesian approach for suppression of limited angular sampling artifacts in single particle 3D reconstruction

In the single particle reconstruction, the initial 3D structure often suffers from the limited angular sampling artifact. Selecting 2D class averages of particle images generally improves the accuracy and efficiency of the reference-free 3D angle estimation, but causes an insufficient angular sampling to fill the information of the target object in the 3D frequency space. Similarly, the initial 3D structure by the random-conical tilt reconstruction has the well-known "missing cone" artifact. Here, we attempted to solve the limited angular sampling problem by sequentially applying maximum a posteriori estimate with expectation maximization algorithm (sMAP-EM). Using both simulated and experimental cryo-electron microscope images, the sMAP-EM was compared to the direct Fourier method on the basis of reconstruction error and resolution. To establish selection criteria of the final regularization weight for the sMAP-EM, the effects of noise level and sampling sparseness on the reconstructions were examined with evenly distributed sampling simulations. The frequency information filled in the missing cone of the conical tilt sampling simulations was assessed by developing new quantitative measurements. All the results of visual and numerical evaluations showed the sMAP-EM performed better than the direct Fourier method, regardless of the sampling method, noise level, and sampling sparseness. Furthermore, the frequency domain analysis demonstrated that the sMAP-EM can fill the meaningful information in the unmeasured angular space without detailed a priori knowledge of the objects. The current research demonstrated that the sMAP-EM has a high potential to facilitate the determination of 3D protein structures at near atomic-resolution.

Reconstruction of 4-D dynamic SPECT images from inconsistent projections using a Spline initialized FADS algorithm (SIFADS)

In this paper, we propose and validate an algorithm of extracting voxel-by-voxel time activity curves directly from inconsistent projections applied in dynamic cardiac SPECT. The algorithm was derived based on factor analysis of dynamic structures (FADS) approach and imposes prior information by applying several regularization functions with adaptively changing relative weighting. The anatomical information of the imaged subject was used to apply the proposed regularization functions adaptively in the spatial domain. The algorithm performance is validated by reconstructing dynamic datasets simulated using the NCAT phantom with a range of different input tissue time-activity curves. The results are compared to the spline-based and FADS methods. The validated algorithm is then applied to reconstruct pre-clinical cardiac SPECT data from canine and murine subjects. Images, generated from both simulated and experimentally acquired data confirm the ability of the new algorithm to solve the inverse problem of dynamic SPECT with slow gantry rotation.

Compensation of missing wedge effects with sequential statistical reconstruction in electron tomography

Electron tomography (ET) of biological samples is used to study the organization and the structure of the whole cell and subcellular complexes in great detail. However, projections cannot be acquired over full tilt angle range with biological samples in electron microscopy. ET image reconstruction can be considered an ill-posed problem because of this missing information. This results in artifacts, seen as the loss of three-dimensional (3D) resolution in the reconstructed images. The goal of this study was to achieve isotropic resolution with a statistical reconstruction method, sequential maximum a posteriori expectation maximization (sMAP-EM), using no prior morphological knowledge about the specimen. The missing wedge effects on sMAP-EM were examined with a synthetic cell phantom to assess the effects of noise. An experimental dataset of a multivesicular body was evaluated with a number of gold particles. An ellipsoid fitting based method was developed to realize the quantitative measures elongation and contrast in an automated, objective, and reliable way. The method statistically evaluates the sub-volumes containing gold particles randomly located in various parts of the whole volume, thus giving information about the robustness of the volume reconstruction. The quantitative results were also compared with reconstructions made with widely-used weighted backprojection and simultaneous iterative reconstruction technique methods. The results showed that the proposed sMAP-EM method significantly suppresses the effects of the missing information producing isotropic resolution. Furthermore, this method improves the contrast ratio, enhancing the applicability of further automatic and semi-automatic analysis. These improvements in ET reconstruction by sMAP-EM enable analysis of subcellular structures with higher three-dimensional resolution and contrast than conventional methods.

PET reconstruction artifact can be minimized by using sinogram correction and filtered back-projection technique

Filtered Back-Projection (FBP) has become an outdated image reconstruction technique in new-generation positron emission tomography (PET)/computed tomography (CT) scanners. Iterative reconstruction used in all new-generation PET scanners is a much improved reconstruction technique. Though a well-calibrated PET system can only be used for clinical imaging in few situations like ours, when compromised PET scanner with one PET module bypassed was used for PET acquisition, FBP with sinogram correction proved to be a better reconstruction technique to minimize streak artifact present in the image reconstructed by the iterative technique.

SNARK09 - a software package for reconstruction of 2D images from 1D projections

The problem of reconstruction of slices and volumes from 1D and 2D projections has arisen in a large number of scientific fields (including computerized tomography, electron microscopy, X-ray microscopy, radiology, radio astronomy and holography). Many different methods (algorithms) have been suggested for its solution. In this paper we present a software package, SNARK09, for reconstruction of 2D images from their 1D projections. In the area of image reconstruction, researchers often desire to compare two or more reconstruction techniques and assess their relative merits. SNARK09 provides a uniform framework to implement algorithms and evaluate their performance. It has been designed to treat both parallel and divergent projection geometries and can either create test data (with or without noise) for use by reconstruction algorithms or use data collected by another software or a physical device. A number of frequently-used classical reconstruction algorithms are incorporated. The package provides a means for easy incorporation of new algorithms for their testing, comparison and evaluation. It comes with tools for statistical analysis of the results and ten worked examples.

Distributed MLEM: an iterative tomographic image reconstruction algorithm for distributed memory architectures

The processing speed for positron emission tomography (PET) image reconstruction has been greatly improved in recent years by simply dividing the workload to multiple processors of a graphics processing unit (GPU). However, if this strategy is generalized to a multi-GPU cluster, the processing speed does not improve linearly with the number of GPUs. This is because large data transfer is required between the GPUs after each iteration, effectively reducing the parallelism. This paper proposes a novel approach to reformulate the maximum likelihood expectation maximization (MLEM) algorithm so that it can scale up to many GPU nodes with less frequent inter-node communication. While being mathematically different, the new algorithm maximizes the same convex likelihood function as MLEM, thus converges to the same solution. Experiments on a multi-GPU cluster demonstrate the effectiveness of the proposed approach.

EXACT GRAYSCALE IMAGE RECONSTRUCTION FROM PROJECTIONS

This paper proposes a new method for the exact reconstruction of gray-scale images from projections. The image projections construct an accumulator array, which is used afterwards to reconstruct the original grayscale image by applying the proposed decomposition algorithm. The proposed method determines the number of projections and the number of rays in each projection that are required in order to achieve the reconstruction. These two parameters also define the dimensions of the accumulator array. Using an accumulator array with proper dimensions ensures that there is always a unique characteristic sample for each pixel, which is used during the reconstruction process to extract the pixel's grayscale value. During the reconstruction phase, the sinusoidal contribution of each pixel is removed from the accumulator array. At the end of the decomposition process the accumulator array becomes empty and the original image is exactly reconstructed. The experimental results confirm the robustness and efficiency of the proposed method.

Automatic left ventricular contour extraction from cardiac magnetic resonance images using cantilever beam and random walk approach

Heart failure is a well-known debilitating disease. From clinical point of view, segmentation of left ventricle (LV) is important in a cardiac magnetic resonance (CMR) image. Accurate parameters are desired for better diagnosis. Proper and fast image segmentation of LV is of paramount importance prior to estimation of these parameters. We prefer random walk approach over other existing techniques due to two of its advantages: (1) robustness to noise and, (2) it does not require any special condition to work. Performance of the method solely depends on the selection of initial seed and parameter β. Problems arise while applying this method to different kind of CMR images bearing different ischemia. It is due due to their implicit geometry definitions unlike general images, where the boundary of LV in the image is not available in an explicit form. This type of images bear multi-labeled LV and the manual seed selection in these images introduces variability in the results. In view of this, the paper presents two modifications in the algorithm: (1) automatic seed selection and, (2) automatic estimation of β from the image. The highlight of our method is its ability to succeed with minimum number of initial seeds.

Model-based respiratory motion compensation for emission tomography image reconstruction

In emission tomography imaging, respiratory motion causes artifacts in lungs and cardiac reconstructed images, which lead to misinterpretations, imprecise diagnosis, impairing of fusion with other modalities, etc. Solutions like respiratory gating, correlated dynamic PET techniques, list-mode data based techniques and others have been tested, which lead to improvements over the spatial activity distribution in lungs lesions, but which have the disadvantages of requiring additional instrumentation or the need of discarding part of the projection data used for reconstruction. The objective of this study is to incorporate respiratory motion compensation directly into the image reconstruction process, without any additional acquisition protocol consideration. To this end, we propose an extension to the maximum likelihood expectation maximization (MLEM) algorithm that includes a respiratory motion model, which takes into account the displacements and volume deformations produced by the respiratory motion during the data acquisition process. We present results from synthetic simulations incorporating real respiratory motion as well as from phantom and patient data.

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.

Discretization of the radon transform and of its inverse by spline convolutions

We present an explicit formula for B-spline convolution kernels; these are defined as the convolution of several B-splines of variable widths h(i) and degrees n(i). We apply our results to derive spline-convolution-based algorithms for two closely related problems: the computation of the Radon transform and of its inverse. First, we present an efficient discrete implementation of the Radon transform that is optimal in the least-squares sense. We then consider the reverse problem and introduce a new spline-convolution version of the filtered back-projection algorithm for tomographic reconstruction. In both cases, our explicit kernel formula allows for the use of high-degree splines; these offer better approximation performance than the conventional lower-degree formulations (e.g., piecewise constant or piecewise linear models). We present multiple experiments to validate our approach and to find the parameters that give the best tradeoff between image quality and computational complexity. In particular, we find that it can be computationally more efficient to increase the approximation degree than to increase the sampling rate.

Fast EM-like methods for maximum "a posteriori" estimates in emission tomography

The maximum-likelihood (ML) approach in emission tomography provides images with superior noise characteristics compared to conventional filtered backprojection (FBP) algorithms. The expectation-maximization (EM) algorithm is an iterative algorithm for maximizing the Poisson likelihood in emission computed tomography that became very popular for solving the ML problem because of its attractive theoretical and practical properties. Recently, (Browne and DePierro, 1996 and Hudson and Larkin, 1994) block sequential versions of the EM algorithm that take advantage of the scanner's geometry have been proposed in order to accelerate its convergence. In Hudson and Larkin, 1994, the ordered subsets EM (OS-EM) method was applied to the ML problem and a modification (OS-GP) to the maximum a posteriori (MAP) regularized approach without showing convergence. In Browne and DePierro, 1996, we presented a relaxed version of OS-EM (RAMLA) that converges to an ML solution. In this paper, we present an extension of RAMLA for MAP reconstruction. We show that, if the sequence generated by this method converges, then it must converge to the true MAP solution. Experimental evidence of this convergence is also shown. To illustrate this behavior we apply the algorithm to positron emission tomography simulated data comparing its performance to OS-GP.

Interior-point methodology for 3-D PET reconstruction

Interior-point methods have been successfully applied to a wide variety of linear and nonlinear programming applications. This paper presents a class of algorithms, based on path-following interior-point methodology, for performing regularized maximum-likelihood (ML) reconstructions on three-dimensional (3-D) emission tomography data. The algorithms solve a sequence of subproblems that converge to the regularized maximum likelihood solution from the interior of the feasible region (the nonnegative orthant). We propose two methods, a primal method which updates only the primal image variables and a primal-dual method which simultaneously updates the primal variables and the Lagrange multipliers. A parallel implementation permits the interior-point methods to scale to very large reconstruction problems. Termination is based on well-defined convergence measures, namely, the Karush-Kuhn-Tucker first-order necessary conditions for optimality. We demonstrate the rapid convergence of the path-following interior-point methods using both data from a small animal scanner and Monte Carlo simulated data. The proposed methods can readily be applied to solve the regularized, weighted least squares reconstruction problem.

Prediction of human observer performance by numerical observers: an experimental study

Numerical observers are investigated for predicting the outcome of a free-response human observer study involving the detection of simulated pulmonary nodules in images reconstructed from low-dose computed tomography projection data by use of several reconstruction algorithms. A new way of calculating the figure of merit of a numerical observer is proposed wherein the detectability of signals in a particular image depends on the noise properties associated with that image and not the other images in the data set. The resulting variants of numerical observers are found to perform better than their traditional counterparts. In particular, the imagewise variant of the region-of-interest observer is found to predict best the rank ordering of algorithms by human observers for the free-response task.

Accelerated EM reconstruction in total-body PET: potential for improving tumour detectability

Total-body positron emission tomography (PET) is a useful diagnostic tool for evaluating malignant disease. However, tumour detection is limited by image artefacts due to the lack of attenuation correction and noise. Attenuation correction may be possible using transmission data acquired after or simultaneously with emission data. Despite the elimination of attenuation artefacts, however, tumour detection is still hampered by noise, which is amplified during image reconstruction by filtered backprojection (FBP). We have investigated, as an alternative to FBP, an accelerated expectation maximization (EM) algorithm for its potential to improve tumour detectability in total-body PET. Signal to noise ratio (SNR), calculated for a tumour with respect to the surrounding background, is used as a figure of merit. A software tumour phantom, with conditions typical of those encountered in a total-body PET study using simultaneous acquisition, is used to optimize and compare various reconstruction approaches. Accelerated EM reconstruction followed by two-dimensional filtering is shown to yield significantly higher SNR than FBP for a range of tumour sizes, concentrations and counting statistics (deltaSNR = 6.3 +/- 3.9, p < 0.001). The methods developed are illustrated by examples derived from physical phantom and patient data.

Evaluation of task-oriented performance of several fully 3D PET reconstruction algorithms

The relative performance of five fully 3D PET reconstruction algorithms is evaluated. The algorithms are a filtered backprojection (FBP) method and two variants each of the EM-ML and ART iterative methods. For each of the iterative methods, one variant makes use of voxels and the other makes use of 'blobs' (spherically symmetric functions smoothly decaying to zero at their boundaries) as basis functions in its discrete reconstruction model. The methods are evaluated from the point of view of the efficacy of the reconstructions produced by them for three typical medical tasks--estimation of the average activity inside specific regions of interest, detection of hot spots, and detection of cold spots. A free parameter is allowed in the description of each of the five algorithms; the parameters are determined by a training process during which a value of the free parameter is selected which (nearly) maximizes a technical figure of merit. Such training and the actual comparative evaluation is done by making use of randomly generated phantoms and their projection data. The methodology allows assignation of levels of statistical significance to claims of the relative superiority of one algorithm over another for a particular task. We find that using blobs as basis functions in the iterative algorithms is definitely advantageous over using voxels. This result has high statistical significance. (We also include a visual illustration of it.) Comparing FBP, EM-ML using blobs, and ART using blobs, we do not find a clear difference in the overall performance of the investigated variants of the methods. If anything, our results suggest that ART using blobs may be the most efficacious of the three.

Image restoration and reconstruction with a Bayesian approach

We have extended Johnson's Bayesian method for image restoration and reconstruction by introducing diagonal line sites, using symmetric neighborhood configurations, and employing an additional hyperparameter for estimation of line sites. A general formulation for arbitrary neighborhood configurations was derived. The major part of this paper deals with the conduct of computer simulations intended to examine the effect of the hyperparameters, the diagonal line sites, and the size of the neighborhood configuration on the performance of the proposed Bayesian method. We show that, for optimal performance, distinct hyperparameters should be used for the intensity sites and line sites. The results also suggest that a large neighborhood configuration should be used. By comparing the near-optimal restored images, we demonstrated that the use of diagonal line sites, along with the symmetric configurations thus made possible, can effectively remove the blocky edge artifacts and produce images of better quality. When the method was applied to positron emission tomography (PET) image reconstruction, our results showed that the quality of the reconstructed images was improved for both computer-simulated and real patient PET data.

Image reconstruction for dynamic PET based on low-order approximation and restoration of the sinogram

Many image-reconstruction methods have been proposed to improve the spatial resolution of positron emission tomography (PET) images and, thus, to produce better quantification. However, these techniques, which are designed for static images, may be inadequate for good reconstruction from dynamic data. We present a simple, but effective, reconstruction approach intended specifically for dynamic studies. First, the level of noise in dynamic PET data is reduced by smoothing along the time axis using a low-order approximation. Next, the denoised sinograms are restored spatially by the method of projections onto convex sets. Finally, images are reconstructed from the restored sinograms by ordinary filtered backprojection. We present experimental results that demonstrate substantial improvements in region-of-interest quantification in actual and simulated dopamine D-2 neuroreceptor-imaging studies of a monkey brain.

Noise analysis of MAP-EM algorithms for emission tomography

The ability to theoretically model the propagation of photon noise through PET and SPECT tomographic reconstruction algorithms is crucial in evaluating the reconstructed image quality as a function of parameters of the algorithm. In a previous approach for the important case of the iterative ML-EM (maximum-likelihood-expectation-maximization) algorithm, judicious linearizations were used to model theoretically the propagation of a mean image and a covariance matrix from one iteration to the next. Our analysis extends this approach to the case of MAP (maximum a posteriori)-EM algorithms, where the EM approach incorporates prior terms. We analyse in detail two cases: a MAP-EM algorithm incorporating an independent gamma prior, and a one-step-late (OSL) version of a MAP-EM algorithm incorporating a multivariate Gaussian prior, for which familiar smoothing priors are special cases. To validate our theoretical analyses, we use a Monte Carlo methodology to compare, at each iteration, theoretical estimates of mean and covariance with sample estimates, and show that the theory works well in practical situations where the noise and bias in the reconstructed images do not assume extreme values.

Bayesian image reconstruction for emission tomography based on median root prior

The aim of the present study was to investigate a new type of Bayesian one-step late reconstruction method which utilizes a median root prior (MRP). The method favours images which have locally monotonous radioactivity concentrations. The new reconstruction algorithm was applied to ideal simulated data, phantom data and some patient examinations with PET. The same projection data were reconstructed with filtered back-projection (FBP) and maximum likelihood-expectation maximization (ML-EM) methods for comparison. The MRP method provided good-quality images with a similar resolution to the FBP method with a ramp filter, and at the same time the noise properties were as good as with Hann-filtered FBP images. The typical artefacts seen in FBP reconstructed images outside of the object were completely removed, as was the grainy noise inside the object. Quantitatively, the resulting average regional radioactivity concentrations in a large region of interest in images produced by the MRP method corresponded to the FBP and ML-EM results but at the pixel by pixel level the MRP method proved to be the most accurate of the tested methods. In contrast to other iterative reconstruction methods, e.g. ML-EM, the MRP method was not sensitive to the number of iterations nor to the adjustment of reconstruction parameters. Only the Bayesian parameter beta had to be set. The proposed MRP method is much more simple to calculate than the methods described previously, both with regard to the parameter settings and in terms of general use. The new MRP reconstruction method was shown to produce high-quality quantitative emission images with only one parameter setting in addition to the number of iterations.

The EM algorithm in medical imaging

This article outlines the statistical developments that have taken place in the use of the EM algorithm in emission and transmission tomography during the past decade or so. We discuss the statistical aspects of the modelling of the projection data for both the emission and transmission cases and define the relevant probability models. This leads to the use of the method of maximum likelihood as a means of estimating the relevant unknown parameters within a given region of a patient's body and to the use of the EM algorithm to compute the reconstruction. Various different types of EM algorithm are discussed, including the SAGE algorithms of Fessler and Hero. The limitations of the EM algorithm, per se, are covered and the need for regularization is stressed. A number of different methods for penalizing the likelihood are described and a number of algorithms for the computation of the penalized EM reconstruction are discussed.

Bayesian reconstruction of PET images: methodology and performance analysis

We describe a practical statistical methodology for the reconstruction of PET images. Our approach is based on a Bayesian formulation of the imaging problem. The data are modelled as independent Poisson random variables and the image is modelled using a Markov random field smoothing prior. We describe a sequence of calibration procedures which are performed before reconstruction: (i) calculation of accurate attenuation correction factors from re-projected Bayesian reconstructions of the transmission image; (ii) estimation of the mean of the randoms component in the data; and (iii) computation of the scatter component in the data using a Klein-Nishina-based scatter estimation method. The Bayesian estimate of the PET image is then reconstructed using a pre-conditioned conjugate gradient method. We performed a quantitation study with a multi-compartment chest phantom in a Siemens/CTI ECAT931 system. Using 40 1 min frames, we computed the ensemble mean and variance over several regions of interest from images reconstructed using the Bayesian and a standard filtered backprojection (FBP) protocol. The values for the region of interest were compared with well counter data for each compartment. These results show that the Bayesian protocol can produce substantial improvements in relative quantitation over the standard FBP protocol, particularly when short transmission scans are used. An example showing the application of the method to a clinical chest study is also given.

Relevance of statistically significant differences between reconstruction algorithms

When comparing reconstruction algorithms, differences in figures of performance merit that are too small to be of any practical relevance may still be statistically significant. We formalize the notion of "relevance" and propose an evaluation methodology in which statistical significance is retained for relevant improvements, but not for irrelevant ones.

A row-action alternative to the EM algorithm for maximizing likelihood in emission tomography

The maximum likelihood (ML) approach to estimating the radioactive distribution in the body cross section has become very popular among researchers in emission computed tomography (ECT) since it has been shown to provide very good images compared to those produced with the conventional filtered backprojection (FBP) algorithm. The expectation maximization (EM) algorithm is an often-used iterative approach for maximizing the Poisson likelihood in ECT because of its attractive theoretical and practical properties. Its major disadvantage is that, due to its slow rate of convergence, a large amount of computation is often required to achieve an acceptable image. Here, the authors present a row-action maximum likelihood algorithm (RAMLA) as an alternative to the EM algorithm for maximizing the Poisson likelihood in ECT. The authors deduce the convergence properties of this algorithm and demonstrate by way of computer simulations that the early iterates of RAMLA increase the Poisson likelihood in ECT at an order of magnitude faster that the standard EM algorithm. Specifically, the authors show that, from the point of view of measuring total radionuclide uptake in simulated brain phantoms, iterations 1, 2, 3, and 4 of RAMLA perform at least as well as iterations 45, 60, 70, and 80, respectively, of EM. Moreover, the authors show that iterations 1, 2, 3, and 4 of RAMLA achieve comparable likelihood values as iterations 45, 60, 70, and 80, respectively, of EM. The authors also present a modified version of a recent fast ordered subsets EM (OS-EM) algorithm and show that RAMLA is a special case of this modified OS-EM. Furthermore, the authors show that their modification converges to a ML solution whereas the standard OS-EM does not.

A unified approach to statistical tomography using coordinate descent optimization

Over the past years there has been considerable interest in statistically optimal reconstruction of cross-sectional images from tomographic data. In particular, a variety of such algorithms have been proposed for maximum a posteriori (MAP) reconstruction from emission tomographic data. While MAP estimation requires the solution of an optimization problem, most existing reconstruction algorithms take an indirect approach based on the expectation maximization (EM) algorithm. We propose a new approach to statistically optimal image reconstruction based on direct optimization of the MAP criterion. The key to this direct optimization approach is greedy pixel-wise computations known as iterative coordinate decent (ICD). We propose a novel method for computing the ICD updates, which we call ICD/Newton-Raphson. We show that ICD/Newton-Raphson requires approximately the same amount of computation per iteration as EM-based approaches, but the new method converges much more rapidly (in our experiments, typically five to ten iterations). Other advantages of the ICD/Newton-Raphson method are that it is easily applied to MAP estimation of transmission tomograms, and typical convex constraints, such as positivity, are easily incorporated.

An expectation maximization reconstruction algorithm for emission tomography with non-uniform entropy prior

A Bayesian image reconstruction algorithm is proposed for emission tomography. It incorporates the Poisson nature of the noise in the projection data and uses a non-uniform entropy as an a priori probability distribution of the image in a maximum a posteriori (MAP) approach. The expectation maximization (EM) method was applied to find the MAP estimator. The Newton-Raphson numerical method whose convergence and positive solutions are proven, was used to solve the EM problem. The prior mean at iteration k was determined by smoothing the image obtained at iteration k-1. Comparisons between the ML and the MAP algorithm were carried out with a numerical phantom that contains a narrow valley region. The ML solution after 50 iterations was chosen as the initial solution for the MAP algorithm, since the global performance of the ML algorithm deteriorates with increasing number of iterations while its local performance in the valley region is always improving. The resulting algorithm is a compromise between ML who has the best local performance in the valley region and the MAP who has the best global performance.

Objective assessment of image quality. II. Fisher information, Fourier crosstalk, and figures of merit for task performance

Figures of merit for image quality are derived on the basis of the performance of mathematical observers on specific detection and estimation tasks. The tasks include detection of a known signal superimposed on a known background, detection of a known signal on a random background, estimation of Fourier coefficients of the object, and estimation of the integral of the object over a specified region of interest. The chosen observer for the detection tasks is the ideal linear discriminant, which we call the Hotelling observer. The figures of merit are based on the Fisher information matrix relevant to estimation of the Fourier coefficients and the closely related Fourier crosstalk matrix introduced earlier by Barrett and Gifford [Phys. Med. Biol. 39, 451 (1994)]. A finite submatrix of the infinite Fisher information matrix is used to set Cramer-Rao lower bounds on the variances of the estimates of the first N Fourier coefficients. The figures of merit for detection tasks are shown to be closely related to the concepts of noise-equivalent quanta (NEQ) and generalized NEQ, originally derived for linear, shift-invariant imaging systems and stationary noise. Application of these results to the design of imaging systems is discussed.

Penalized maximum-likelihood image reconstruction using space-alternating generalized EM algorithms

Most expectation-maximization (EM) type algorithms for penalized maximum-likelihood image reconstruction converge slowly, particularly when one incorporates additive background effects such as scatter, random coincidences, dark current, or cosmic radiation. In addition, regularizing smoothness penalties (or priors) introduce parameter coupling, rendering intractable the M-steps of most EM-type algorithms. This paper presents space-alternating generalized EM (SAGE) algorithms for image reconstruction, which update the parameters sequentially using a sequence of small "hidden" data spaces, rather than simultaneously using one large complete-data space. The sequential update decouples the M-step, so the maximization can typically be performed analytically. We introduce new hidden-data spaces that are less informative than the conventional complete-data space for Poisson data and that yield significant improvements in convergence rate. This acceleration is due to statistical considerations, not numerical overrelaxation methods, so monotonic increases in the objective function are guaranteed. We provide a general global convergence proof for SAGE methods with nonnegativity constraints.

A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography

The maximum likelihood (ML) expectation maximization (EM) approach in emission tomography has been very popular in medical imaging for several years. In spite of this, no satisfactory convergent modifications have been proposed for the regularized approach. Here, a modification of the EM algorithm is presented. The new method is a natural extension of the EM for maximizing likelihood with concave priors. Convergence proofs are given.

Statistical models for PET and SPECT data

This article outlines the statistical developments that have taken place in emission tomography during the past decade or so. We discuss the statistical aspects of the modelling of the projection data and define the additive Poisson regression model. This leads to the use of the method of maximum likelihood as a means of estimating the underlying isotope concentration within a given region of a patient's body, and to the use of the EM algorithm to compute the reconstruction. The need for the regulation of the maximum likelihood solution is tackled using Bayesian techniques. A number of algorithms for the computation of regularized solutions are outlined. The issue of parameter estimation is discussed and some open issues are mentioned.

On the relation between the ISRA and the EM algorithm for positron emission tomography

The image space reconstruction algorithm (ISRA) was proposed as a modification of the expectation maximization (EM) algorithm based on physical considerations for application in volume emission computered tomography. As a consequence of this modification, ISRA searches for least squares solutions instead of maximizing Poisson likelihoods as the EM algorithm. It is shown that both algorithms may be obtained from a common mathematical framework. This fact is used to extend ISRA for penalized likelihood estimates.

Acceleration of maximum likelihood reconstruction, using frequency amplification and attenuation compensation

Algorithms that calculate maximum likelihood (ML) and maximum a posteriori solutions using expectation-maximization have been successfully applied to SPECT and PET. These algorithms are appealing because of their solid theoretical basis and their guaranteed convergence. A major drawback is the slow convergence, which results in-long processing times. The authors present 2 new heuristic acceleration methods for maximum likelihood reconstruction of ECT images. The first method incorporates a frequency-dependent amplification in the calculations, to compensate for the low pass filtering of the backprojection operation. In the second method, an amplification factor is incorporated that suppresses the effect of attenuation on the updating factors. Both methods are compared to the 1-dimensional line search method proposed by Lewitt. All 3 methods accelerate the ML algorithm. On the authors' test images, Lewitt's method produced the strongest acceleration of the three individual methods. However, the combination of the frequency amplification with the line search method results in a new algorithm with still better performance. Under certain conditions, an effective frequency amplification can be already achieved by skipping some of the calculations required for ML.

Algebraic reconstruction techniques can be made computationally efficient [positron emission tomography application]

Algebraic reconstruction techniques (ART) are iterative procedures for recovering objects from their projections. It is claimed that by a careful adjustment of the order in which the collected data are accessed during the reconstruction procedure and of the so-called relaxation parameters that are to be chosen in an algebraic reconstruction technique, ART can produce high-quality reconstructions with excellent computational efficiency. This is demonstrated by an example based on a particular (but realistic) medical imaging task, showing that ART can match the performance of the standard expectation-maximization approach for maximizing likelihood (from the point of view of that particular medical task), but at an order of magnitude less computational cost.

Maximum likelihood, least squares, and penalized least squares for PET

The EM algorithm is the basic approach used to maximize the log likelihood objective function for the reconstruction problem in positron emission tomography (PET). The EM algorithm is a scaled steepest ascent algorithm that elegantly handles the nonnegativity constraints of the problem. It is shown that the same scaled steepest descent algorithm can be applied to the least squares merit function, and that it can be accelerated using the conjugate gradient approach. The experiments suggest that one can cut the computation by about a factor of 3 by using this technique. The results are applied to various penalized least squares functions which might be used to produce a smoother image.

Statistical analysis of maximum likelihood estimator images of human brain FDG PET studies

The work presented evaluates the statistical characteristics of regional bias and expected error in reconstructions of real positron emission tomography (PET) data of human brain fluoro-deoxiglucose (FDG) studies carried out by the maximum likelihood estimator (MLE) method with a robust stopping rule, and compares them with the results of filtered backprojection (FBP) reconstructions and with the method of sieves. The task of evaluating radioisotope uptake in regions-of-interest (ROIs) is investigated. An assessment of bias and variance in uptake measurements is carried out with simulated data. Then, by using three different transition matrices with different degrees of accuracy and a components of variance model for statistical analysis, it is shown that the characteristics obtained from real human FDG brain data are consistent with the results of the simulation studies.

The convergence of object dependent resolution in maximum likelihood based tomographic image reconstruction

Study of the maximum likelihood by EM algorithm (ML) with a reconstruction kernel equal to the intrinsic detector resolution and sieve regularization has demonstrated that any image improvements over filtered backprojection (FBP) are a function of image resolution. Comparing different reconstruction algorithms potentially requires measuring and matching the image resolution. Since there are no standard methods for describing the resolution of images from a nonlinear algorithm such as ML, we have defined measures of effective local Gaussian resolution (ELGR) and effective global Gaussian resolution (EGGR) and examined their behaviour in FBP images and in ML images using two different measurement techniques. For FBP these two resolution measures are equal and exhibit the standard convolution behaviour of linear systems. For ML, the FWHM of the ELGR monotonically increased with decreasing Gaussian object size due to slower convergence rates for smaller objects. For the simple simulated phantom used, this resolution dependence is independent of object position. With increasing object size, number of iterations and sieve size the object size dependence of the ELGR decreased. The FWHM of the EGGR converged after approximately 200 iterations, masking the fact that the ELGR for small objects was far from convergence. When FBP is compared to a nonlinear algorithm such as ML, it is recommended that at least the EGGR be matched; for ML this requires more than the number of iterations (e.g., < 100) that are typically run to minimize the mean square error or to satisfy a feasibility or similar stopping criterion. For many tasks, matching the EGGR of ML to FBP images may be insufficient and >> 200 iterations may be needed, particularly for small objects in the ML image because their ELGR has not yet converged.