Ordered subsets Bayesian tomographic reconstruction using 2-D smoothing splines as priors
Introduction
Statistical image reconstruction methods have played an important role in emission tomography since they accurately model the Poisson noise associated with gamma-ray projection data. A significant advance for these methods was made with the introduction of the maximum likelihood (ML) approach using the famous EM algorithm [1], [2]. Since the ML-EM algorithm provides a number of potential benefits, its application has led to the introduction of many related reconstruction algorithms and its value is now well established. Although the application of ML-EM in clinical practice in early days was a difficult problem due to high per-iteration costs and large numbers of iterations, with the dramatically improved speed of computers, it is now inevitably becoming more practical and is indeed commonly available from nuclear medicine equipment suppliers. However, it is important to point out that, while the main limitation in applying EM due to processing time is becoming less of a problem, there is also continuing need for more complex models to characterize the actual physics involved more accurately, which requires more computational load.
Recently, there has been increasing emphasis on accelerating iterative reconstruction methods which can produce an image with much less iterations than the ML-EM algorithm. Some accelerated methods demonstrate a possibility of permitting reconstruction in a clinically acceptable time with a few iterations. For example, the ordered subsets expectation maximization (OS-EM) algorithm proposed by Hudson and Larkin [3] has been successful in accelerating the existing ML-EM algorithm and is continually increasing in its popularity. This is presumably due to the fact that, while retaining the advantages of ML-EM such as accurate modeling of any type of system, it provides an order-of-magnitude acceleration over ML-EM.
Besides OS-EM, a number of other fast algorithms for finding ML or maximum a posteriori (MAP) estimation for emission tomography have also been introduced in the literature [4], [5], [6], [7]. To avoid the slow convergence problem in the traditional EM algorithm, some algorithms are based on direct maximization [7] of the objective function rather than on the EM algorithm. Nevertheless, the OS-EM algorithm has been more extensively used than others and has recently found widespread popularity in clinical practice. In this paper, we are concerned exclusively with the OS approach.
Although the OS-EM algorithm provides good-quality reconstructions after only a few iterations, beyond a critical number of iterations, the noise artifact of ML solutions is magnified due to the inherent instability problem of ML-EM. One can alleviate the instability of OS-EM by applying the OS principle to regularized EM in the context of a Bayesian MAP framework, which allows the incorporation of suitable prior models to regularize the ill-posed nature of the tomographic inversion problem. A significant body of work for Bayesian tomographic reconstruction has been developed over the last decade. Early methods focused on the regularization of the unstable ML-EM algorithm. Later, the Bayesian approach was also viewed as a means of incorporating actual, known information regarding the local spatial character of the source. A host of different formulations to model the Bayesian priors have been proposed in the literature [8], [9], [10], [11], [12], [13]; some of these implicitly model the underlying radionuclide density as globally smooth, and others extend the smoothness model by allowing for spatial discontinuities. While the former is associated with a smoothing penalty that is a quadratic function of nearby pixel, the latter uses a nonquadratic penalty function.
Hudson and Larkin [3] applied the ordered subsets principle to the one-step-late (OSL) algorithm proposed by Green [8] using a nonquadratic log-cosh penalty function and demonstrated that the instability of OS-EM could be reduced. Unfortunately, however, this method seems to be less than ideal due to the following reasons: (i) although the OSL approach provides an easy way of implementing complex nonquadratic priors as well as simple quadratic priors, it is known to be unstable when the smoothing parameter that weights the prior relative to the likelihood is relatively large; (ii) the use of nonquadratic priors requires estimation of additional ‘hyperparameters’ to adjust penalty functions for edge preservation. With regards to point (i), a number of other gradient-based descent algorithms exist, which can provide a stable solution for nonquadratic convex priors as well as quadratic priors. However, the performance of such algorithms depends highly on the choice of step sizes or line-search algorithms. Despite the good performance under properly chosen optimization algorithms and hyperparameters, any enthusiasm for using nonquadratic priors is usually tempered in clinical practice by the fact that most methods for estimating the hyperparameters involved in such priors are characterized by severe computational cost.
We note here that, by relaxing the requirement of preserving sharp edges, which are hardly found in clinical emission data, and using instead useful quadratic smoothing spline priors, the difficulties in Bayesian methods described above may be avoided, and the OS algorithm may become more practical. The remainder of this paper formulates MAP-EM algorithms with two-dimensional (2-D) smoothing spline priors, develops an optimization algorithm using a method of iterated conditional modes (ICM), derives an accelerated algorithm using the OS approach, and presents our experimental results and conclusion.
Section snippets
Formulation of MAP-EM with smoothing spline priors
The MAP approach in the context of a Bayesian framework is to estimate the underlying source field f by maximizing the posterior probability, given aswhere f and g are the 2-D vector fields for the source intensities and projection data, respectively, and F and G are the associated random fields. In Eq. (1) Pr(G=g ∣ F=f) is the likelihood which is Poisson distributed, and Pr(F=f) is the prior distribution.
In Bayesian approaches the prior is usually modeled as
Ordered subsets
In order to accelerate the convergence speed of the ICM algorithm in Eq. (13), we apply the OS idea [3] which was used for accelerating the EM algorithm.
While the standard EM algorithm, which is given by,uses the estimation of all projections and calculates the ratios between measured and estimated values for all projections in the backprojection step , the OS-EM algorithm groups
Simulations and results
We first performed 2-D simulation studies with projection data from two 128×128 digital phantoms, designated A and B for convenience, with 128 projection angles over 180° and 128 detector bins at each projection. Phantom A (Fig. 2(a)) is a digital Hoffman brain phantom with activities of 4:1:0 in grey matter, white matter, and cerebro-spinal fluid (CSF), respectively. Phantom B (Fig. 2(b)) was derived from a digitized rhesus monkey autoradiograph [29].
The motivation of using the autoradiograph
Conclusion
We proposed a method to apply the ordered subsets principle to the Bayesian reconstruction problem for emission tomography. Although the OS principle can be applied to the OSL algorithm so that OS-EM can be extended to a Bayesian MAP approach, our method applied to a version of the ICM method provides more robust results in that the convergence is not affected by a large value of the smoothing parameter.
One could also consider applying the OS algorithm to other MAP-EM algorithms which offer
Acknowledgements
This work was supported by the Nuclear Energy Research Grant from Korea Ministry of Science and Technology.
References (30)
Multilevel computational processes for visual surface reconstruction
Comput. Vis. Graph. Image Processing
(1983)- et al.
A statistical model for positron emission tomography
J. Am. Stat. Assoc.
(1985) - et al.
EM reconstruction algorithms for emission and transmission tomography
J. Comput. Assist. Tomogr.
(1984) - et al.
Accelerated image reconstruction using ordered subsets of projection data
IEEE Trans. Med. Imaging
(1994) A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography
IEEE Trans. Med. Imaging
(1995)- et al.
A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography
IEEE Trans. Image Processing
(1996) - et al.
Penalized maximum-likelihood image reconstruction using space-alternating generalized EM algorithms
IEEE Trans. Image Processing
(1996) - et al.
A unified approach to statistical tomography using coordinate descent optimization
IEEE Trans. Image Processing
(1996) Bayesian reconstructions from emission tomography data using a modified EM algorithm
IEEE Trans. Med. Imaging
(1990)- et al.
Statistical methods for tomographic image reconstruction
Bull. Int. Stat. Inst.
(1987)