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

Block-iterative Fisher scoring algorithms for maximum penalized likelihood image reconstruction in emission tomography

This paper introduces and evaluates a block-iterative Fisher scoring (BFS) algorithm. The algorithm provides regularized estimation in tomographic models of projection data with Poisson variability. Regularization is achieved by penalized likelihood with a general quadratic penalty. Local convergence of the block-iterative algorithm is proven under conditions that do not require iteration dependent relaxation. We show that, when the algorithm converges, it converges to the unconstrained maximum penalized likelihood (MPL) solution. Simulation studies demonstrate that, with suitable choice of relaxation parameter and restriction of the algorithm to respect nonnegative constraints, the BFS algorithm provides convergence to the constrained MPL solution. Constrained BFS often attains a maximum penalized likelihood faster than other block-iterative algorithms which are designed for nonnegatively constrained penalized reconstruction.

Advances in PET Image Reconstruction

Until recently, the most widely used methods for image reconstruction were direct analytic techniques. Iterative techniques, although computationally much more intensive, produce improved images (principally arising from more accurate modeling of the acquired projection data), enabling these techniques to replace analytic techniques not only in research settings but also in the clinic. This article offers an overview of image reconstruction theory and algorithms for PET, with a particular emphasis on statistical iterative reconstruction techniques. Future directions for image reconstruction in PET are considered, which concern mainly improving the modeling of the data acquisition process and task-specific specification of the parameters to be estimated in image reconstruction.

Noise and resolution in images reconstructed with FBP and OSC algorithms for CT

This paper presents a comparison between an analytical and a statistical iterative reconstruction algorithm for computed transmission tomography concerning their noise and resolution performance. The reconstruction of two-dimensional images from simulated fan-beam transmission data is performed with a filtered back-projection (FBP) type reconstruction and an iterative ordered subsets convex (OSC) maximum-likelihood method. A special software phantom, which allows measuring the resolution and noise in a nonambiguous way, is used to simulate transmission tomography scans with different signal-to-noise ratios (SNR). The noise and modulation transfer function is calculated for FBP and OSC reconstruction at several positions, distributed over the field-of-view (FOV). The reconstruction with OSC using different numbers of subsets shows an inverse linear relation to the number of iterations that are necessary to reach a certain resolution and SNR, i.e., increasing the number of subsets by a factor x reduces the number of required iterations by the same factor. The OSC algorithm is able to achieve a nearly homogeneous high resolution over the whole FOV, which is not achieved with FBP. The OSC method achieves a lower level of noise compared with FBP at the same resolution. The reconstruction with OSC can save a factor of up to nine of x-ray dose compared with FBP in the investigated range of noise levels.

An overview of fast convergent ordered-subsets reconstruction methods for emission tomography based on the incremental EM algorithm

Statistical reconstruction has become popular in emission computed tomography but suffers slow convergence (to the MAP or ML solution). Methods proposed to address this problem include the fast but non-convergent OSEM and the convergent RAMLA [1] for the ML case, and the convergent BSREM [2], relaxed OS-SPS and modified BSREM [3] for the MAP case. The convergent algorithms required a user-determined relaxation schedule. We proposed fast convergent OS reconstruction algorithms for both ML and MAP cases, called COSEM (Complete-data OSEM), which avoid the use of a relaxation schedule while maintaining convergence. COSEM is a form of incremental EM algorithm. Here, we provide a derivation of our COSEM algorithms and demonstrate COSEM using simulations. At early iterations, COSEM-ML is typically slower than RAMLA, and COSEM-MAP is typically slower than optimized BSREM while remaining much faster than conventional MAP-EM. We discuss how COSEM may be modified to overcome these limitations.

Optimization algorithms and weighting factors for analysis of dynamic PET studies

Positron emission tomography (PET) pharmacokinetic analysis involves fitting of measured PET data to a PET pharmacokinetic model. The fitted parameters may, however, suffer from bias or be unrealistic, especially in the case of noisy data. There are many optimization algorithms, each having different characteristics. The purpose of the present study was to evaluate (1) the performance of different optimization algorithms and (2) the effects of using incorrect weighting factors during optimization in terms of both accuracy and reproducibility of fitted PET pharmacokinetic parameters. In this study, the performance of commonly used optimization algorithms (i.e. interior-reflective Newton methods) and a simulated annealing (SA) method was evaluated. This SA algorithm, known as basin hopping, was modified for the present application. In addition, optimization was performed using various weighting factors. Algorithms and effects of using incorrect weighting factors were studied using both simulated and clinical time-activity curves (TACs). Input data, taken from [(15)O]H(2)O, [(11)C]flumazenil and [(11)C](R)-PK11195 studies, were used to simulate time-activity curves at various variance levels (0-15% COV). Clinical evaluation was based on studies with the same three tracers. SA was able to produce accurate results without the need for selecting appropriate starting values for (kinetic) parameters, in contrast to the interior-reflective Newton method. The latter gave biased results unless it was modified to allow for a range of starting values for the different parameters. For patient studies, where large variability is expected, both SA and the extended Newton method provided accurate results. Simulations and clinical assessment showed similar results for the evaluation of different weighting models in that small to intermediate mismatches between data variance and weighting factors did not significantly affect the outcome of the fits. Large errors were observed only when the mismatch between weighting model and data variance was large. It is concluded that selection of specific optimization algorithms and weighting factors can have a large effect on the accuracy and precision of PET pharmacokinetic analysis. Apart from carefully selecting appropriate algorithms and variance models, further improvement in accuracy might be obtained by using noise reducing strategies, such as wavelet filtering, provided that these methods do not introduce significant bias.

Penalized maximum-likelihood image reconstruction for lesion detection

Detecting cancerous lesions is one major application in emission tomography. In this paper, we study penalized maximum-likelihood image reconstruction for this important clinical task. Compared to analytical reconstruction methods, statistical approaches can improve the image quality by accurately modelling the photon detection process and measurement noise in imaging systems. To explore the full potential of penalized maximum-likelihood image reconstruction for lesion detection, we derived simplified theoretical expressions that allow fast evaluation of the detectability of a random lesion. The theoretical results are used to design the regularization parameters to improve lesion detectability. We conducted computer-based Monte Carlo simulations to compare the proposed penalty function, conventional penalty function, and a penalty function for isotropic point spread function. The lesion detectability is measured by a channelized Hotelling observer. The results show that the proposed penalty function outperforms the other penalty functions for lesion detection. The relative improvement is dependent on the size of the lesion. However, we found that the penalty function optimized for a 5 mm lesion still outperforms the other two penalty functions for detecting a 14 mm lesion. Therefore, it is feasible to use the penalty function designed for small lesions in image reconstruction, because detection of large lesions is relatively easy.

Iterative reconstruction techniques in emission computed tomography

In emission tomography statistically based iterative methods can improve image quality relative to analytic image reconstruction through more accurate physical and statistical modelling of high-energy photon production and detection processes. Continued exponential improvements in computing power, coupled with the development of fast algorithms, have made routine use of iterative techniques practical, resulting in their increasing popularity in both clinical and research environments. Here we review recent progress in developing statistically based iterative techniques for emission computed tomography. We describe the different formulations of the emission image reconstruction problem and their properties. We then describe the numerical algorithms that are used for optimizing these functions and illustrate their behaviour using small scale simulations.

Convergent incremental optimization transfer algorithms: application to tomography

No convergent ordered subsets (OS) type image reconstruction algorithms for transmission tomography have been proposed to date. In contrast, in emission tomography, there are two known families of convergent OS algorithms: methods that use relaxation parameters, and methods based on the incremental expectation-maximization (EM) approach. This paper generalizes the incremental EM approach by introducing a general framework, "incremental optimization transfer." The proposed algorithms accelerate convergence speeds and ensure global convergence without requiring relaxation parameters. The general optimization transfer framework allows the use of a very broad family of surrogate functions, enabling the development of new algorithms. This paper provides the first convergent OS-type algorithm for (nonconcave) penalized-likelihood (PL) transmission image reconstruction by using separable paraboloidal surrogates (SPS) which yield closed-form maximization steps. We found it is very effective to achieve fast convergence rates by starting with an OS algorithm with a large number of subsets and switching to the new "transmission incremental optimization transfer (TRIOT)" algorithm. Results show that TRIOT is faster in increasing the PL objective than nonincremental ordinary SPS and even OS-SPS yet is convergent.

A modified OSEM algorithm for PET reconstruction using wavelet processing

Ordered subset expectation-maximization (OSEM) method in positron emission tomography (PET) has been very popular recently. It is an iterative algorithm and provides images with superior noise characteristics compared to conventional filtered backprojection (FBP) algorithms. Due to the lack of smoothness in images in OSEM iterations, however, some type of inter-smoothing is required. For this purpose, the smoothing based on the convolution with the Gaussian kernel has been used in clinical PET practices. In this paper, we incorporated a robust wavelet de-noising method into OSEM iterations as an inter-smoothing tool. The proposed wavelet method is based on a hybrid use of the standard wavelet shrinkage and the robust wavelet shrinkage to have edge preserving and robust de-noising simultaneously. The performances of the proposed method were compared with those of the smoothing methods based on the convolution with Gaussian kernel using software phantoms, physical phantoms, and human PET studies. The results demonstrated that the proposed wavelet method provided better spatial resolution characteristic than the smoothing methods based on the Gaussian convolution, while having comparable performance in noise removal.

A fast image reconstruction algorithm based on penalized-likelihood estimate

Statistical iterative methods for image reconstruction like maximum likelihood expectation maximization (ML-EM) are more robust and flexible than analytical inversion methods and allow for accurately modeling the counting statistics and the photon transport during acquisition. They are rapidly becoming the standard for image reconstruction in emission computed tomography. The maximum likelihood approach provides images with superior noise characteristics compared to the conventional filtered back projection algorithm. But a major drawback of the statistical iterative image reconstruction is its high computational cost. In this paper, a fast algorithm is proposed as a modified OS-EM (MOS-EM) using a penalized function, which is applied to the least squares merit function to accelerate image reconstruction and to achieve better convergence. The experimental results show that the algorithm can provide high quality reconstructed images with a small number of iterations.

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.

Statistical list-mode image reconstruction for the high resolution research tomograph

We have investigated statistical list-mode reconstruction applicable to a depth-encoding high resolution research tomograph. An image non-negativity constraint has been employed in the reconstructions and is shown to effectively remove the overestimation bias introduced by the sinogram non-negativity constraint. We have furthermore implemented a convergent subsetized (CS) list-mode reconstruction algorithm, based on previous work (Hsiao et al 2002 Conf. Rec. SPIE Med. Imaging 4684 10-19; Hsiao et al 2002 Conf. Rec. IEEE Int. Symp. Biomed. Imaging 409-12) on convergent histogram OSEM reconstruction. We have demonstrated that the first step of the convergent algorithm is exactly equivalent (unlike the histogram-mode case) to the regular subsetized list-mode EM algorithm, while the second and final step takes the form of additive updates in image space. We have shown that in terms of contrast, noise as well as FWHM width behaviour, the CS algorithm is robust and does not result in limit cycles. A hybrid algorithm based on the ordinary and the convergent algorithms is also proposed, and is shown to combine the advantages of the two algorithms (i.e. it is able to reach a higher image quality in fewer iterations while maintaining the convergent behaviour), making the hybrid approach a good alternative to the ordinary subsetized list-mode EM algorithm.

Regularized image reconstruction algorithms for positron emission tomography

We develop algorithms for obtaining regularized estimates of emission means in positron emission tomography. The first algorithm iteratively minimizes a penalized maximum-likelihood (PML) objective function. It is based on standard de-coupled surrogate functions for the ML objective function and de-coupled surrogate functions for a certain class of penalty functions. As desired, the PML algorithm guarantees nonnegative estimates and monotonically decreases the PML objective function with increasing iterations. The second algorithm is based on an iteration dependent, de-coupled penalty function that introduces smoothing while preserving edges. For the purpose of making comparisons, the MLEM algorithm and a penalized weighted least-squares algorithm were implemented. In experiments using synthetic data and real phantom data, it was found that, for a fixed level of background noise, the contrast in the images produced by the proposed algorithms was the most accurate.

An accelerated convergent ordered subsets algorithm for emission tomography

We propose an algorithm, E-COSEM (enhanced complete-data ordered subsets expectation-maximization), for fast maximum likelihood (ML) reconstruction in emission tomography. E-COSEM is founded on an incremental EM approach. Unlike the familiar OSEM (ordered subsets EM) algorithm which is not convergent, we show that E-COSEM converges to the ML solution. Alternatives to the OSEM include RAMLA, and for the related maximum a posteriori (MAP) problem, the BSREM and OS-SPS algorithms. These are fast and convergent, but require ajudicious choice of a user-specified relaxation schedule. E-COSEM itself uses a sequence of iteration-dependent parameters (very roughly akin to relaxation parameters) to control a tradeoff between a greedy, fast but non-convergent update and a slower but convergent update. These parameters are computed automatically at each iteration and require no user specification. For the ML case, our simulations show that E-COSEM is nearly as fast as RAMLA.

Analysis of lesion detectability in Bayesian emission reconstruction with nonstationary object variability

Bayesian methods based on the maximum a posteriori principle (also called penalized maximum-likelihood methods) have been developed to improve image quality in emission tomography. To explore the full potential of Bayesian reconstruction for lesion detection, we derive simplified theoretical expressions that allow fast evaluation of the detectability of a lesion in Bayesian reconstruction. This work is builded on the recent progress on the theoretical analysis of image properties of statistical reconstructions and the development of numerical observers. We explicitly model the nonstationary variation of the lesion and background without assuming that they are locally stationary. The results can be used to choose the optimum prior parameters for the maximum lesion detectability. The theoretical results are validated using Monte Carlo simulations. The comparisons show good agreement between the theoretical predictions and the Monte Carlo results. We also demonstrate that the lesion detectability can be reliably estimated using one noisy data set.

A unified noise analysis for iterative image estimation

Iterative image estimation methods have been widely used in emission tomography. Accurate estimation of the uncertainty of the reconstructed images is essential for quantitative applications. While both iteration-based noise analysis and fixed-point noise analysis have been developed, current iteration-based results are limited to only a few algorithms that have an explicit multiplicative update equation and some may not converge to the fixed-point result. This paper presents a theoretical noise analysis that is applicable to a wide range of preconditioned gradient-type algorithms. Under a certain condition, the proposed method does not require an explicit expression of the preconditioner. By deriving the fixed-point expression from the iteration-based result, we show that the proposed iteration-based noise analysis is consistent with fixed-point analysis. Examples in emission tomography and transmission tomography are shown. The results are validated using Monte Carlo simulations.

Subset-dependent relaxation in block-iterative algorithms for image reconstruction in emission tomography

This paper presents a row-action maximum likelihood algorithm (RAMLA), in which the relaxation parameter is controlled in such a way that the noise propagation from projection data to the reconstructed image is substantially independent of the access order of the input data (subsets) in each cycle of the sub-iterations. The 'subset-dependent' relaxation parameter lambda(k) (q) is expressed as lambda(k)(q) = beta0/(beta0 + q + gamma k M), where M is the number of angular views, q (0 < or = q < or = M - 1) is the access order of the angular view, k is the iteration number and beta0 and gamma are constants. The constant beta0 deals with the balance of the noise propagation and the constant gamma controls the convergence of iterations. The value of beta0 is determined from the geometrical correlation coefficients among lines of coincidence response. The proposed RAMLA using the subset-dependent (dynamic) relaxation 'dynamic RAMLA (DRAMA)' provides a reasonable signal-to-noise ratio with a satisfactory spatial resolution by a few iterations in the two-dimensional image reconstruction for PET. Dynamic OS-EM (DOSEM) has also been developed, which allows the use of a larger number of subsets (OS level) Msub without loss of signal-to-noise ratio as compared to the conventional OS-EM. DRAMA is a special case of DOSEM, where Msub = M, and it is no more profitable to use DOSEM with a smaller Msub (< M), because DRAMA provides similar performance with the fastest convergence and smallest computer burden. This paper describes the theory, algorithm and the results of the simulation studies on the performance of DRAMA and DOSEM.

Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms

We present two types of globally convergent relaxed ordered subsets (OS) algorithms for penalized-likelihood image reconstruction in emission tomography: modified block sequential regularized expectation-maximization (BSREM) and relaxed OS separable paraboloidal surrogates (OS-SPS). The global convergence proof of the existing BSREM (De Pierro and Yamagishi, 2001) required a few a posteriori assumptions. By modifying the scaling functions of BSREM, we are able to prove the convergence of the modified BSREM under realistic assumptions. Our modification also makes stepsize selection more convenient. In addition, we introduce relaxation into the OS-SPS algorithm (Erdoğan and Fessler, 1999) that otherwise would converge to a limit cycle. We prove the global convergence of diagonally scaled incremental gradient methods of which the relaxed OS-SPS is a special case; main results of the proofs are from (Nedić and Bertsekas, 2001) and (Correa and Lemaréchal, 1993). Simulation results showed that both new algorithms achieve global convergence yet retain the fast initial convergence speed of conventional unrelaxed ordered subsets algorithms.

Relaxed ordered-subset algorithm for penalized-likelihood image restoration

The expectation-maximization (EM) algorithm for maximum-likelihood image recovery is guaranteed to converge, but it converges slowly. Its ordered-subset version (OS-EM) is used widely in tomographic image reconstruction because of its order-of-magnitude acceleration compared with the EM algorithm, but it does not guarantee convergence. Recently the ordered-subset, separable-paraboloidal-surrogate (OS-SPS) algorithm with relaxation has been shown to converge to the optimal point while providing fast convergence. We adapt the relaxed OS-SPS algorithm to the problem of image restoration. Because data acquisition in image restoration is different from that in tomography, we employ a different strategy for choosing subsets, using pixel locations rather than projection angles. Simulation results show that the relaxed OS-SPS algorithm can provide an order-of-magnitude acceleration over the EM algorithm for image restoration. This new algorithm now provides the speed and guaranteed convergence necessary for efficient image restoration.