Ordered subsets Bayesian tomographic reconstruction using 2-D smoothing splines as priors

https://doi.org/10.1016/S0169-2607(02)00112-8Get rights and content

Abstract

The ordered subsets expectation maximization (OS-EM) algorithm has enjoyed considerable interest for accelerating the well-known EM algorithm for emission tomography. The OS principle has also been applied to several regularized EM algorithms, such as nonquadratic convex minimization-based maximum a posteriori (MAP) algorithms. However, most of these methods have not been as practical as OS-EM due to their complex optimization methods and difficulties in hyperparameter estimation. We note here that, by relaxing the requirement of imposing sharp edges and using instead useful quadratic spline priors, solutions are much easier to compute, and hyperparameter calculation becomes less of a problem. In this work, we use two-dimensional smoothing splines as priors and apply a method of iterated conditional modes for the optimization. In this case, step sizes or line-search algorithms necessary for gradient-based descent methods are avoided. We also accelerate the resulting algorithm using the OS approach and propose a principled way of scaling smoothing parameters to retain the strength of smoothing for different subset numbers. Our experimental results show that the OS approach applied to our quadratic MAP algorithms provides a considerable acceleration while retaining the advantages of quadratic spline 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 asPr(F=f|G=g)=Pr(G=g|F=f)Pr(F=f)Pr(G=g),where 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,f̂ijk+1=f̂ijkHtθ,ijHtθ,ijgmnHtθ,mnf̂mnk,uses the estimation of all projections (∑mnHtθ,mnf̂mnk) and calculates the ratios (r̂k=defg/∑mnHtθ,mnf̂mnk) between measured and estimated values for all projections in the backprojection step (∑Htθ,ijr̂k), 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)

  • D. Terzopoulos

    Multilevel computational processes for visual surface reconstruction

    Comput. Vis. Graph. Image Processing

    (1983)
  • Y. Vardi et al.

    A statistical model for positron emission tomography

    J. Am. Stat. Assoc.

    (1985)
  • K. Lange et al.

    EM reconstruction algorithms for emission and transmission tomography

    J. Comput. Assist. Tomogr.

    (1984)
  • H.M. Hudson et al.

    Accelerated image reconstruction using ordered subsets of projection data

    IEEE Trans. Med. Imaging

    (1994)
  • A.D. Pierro

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

    IEEE Trans. Med. Imaging

    (1995)
  • J.A. Browne et al.

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

    IEEE Trans. Image Processing

    (1996)
  • J.A. Fessler et al.

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

    IEEE Trans. Image Processing

    (1996)
  • C. Bouman et al.

    A unified approach to statistical tomography using coordinate descent optimization

    IEEE Trans. Image Processing

    (1996)
  • P.J. Green

    Bayesian reconstructions from emission tomography data using a modified EM algorithm

    IEEE Trans. Med. Imaging

    (1990)
  • S. Geman et al.

    Statistical methods for tomographic image reconstruction

    Bull. Int. Stat. Inst.

    (1987)
  • T. Hebert et al.

    A generalized em algorithm for 3-D Bayesian reconstruction for poisson data using Gibbs priors

    IEEE Trans. Med. Imaging

    (1989)
  • V.E. Johnson et al.

    Bayesian reconstruction of PET images using Gibbs priors

  • D.S. Lalush et al.

    Simulation evaluation of Gibbs prior distributions for use in maximum a posteriori SPECT reconstructions

    IEEE Trans. Med. Imaging

    (1992)
  • S.J. Lee et al.

    Bayesian image reconstruction in SPECT using higher order mechanical models as priors

    IEEE Trans. Med. Imaging

    (1995)
  • J. Besag

    Spatial interaction and the statistical analysis of lattice systems

    J. R. Stat. Soc. Ser. B

    (1974)
  • Cited by (0)

    View full text