There are many practical problems where the observed data are not drawn directly from the density g of real interest, but rather from another distribution derived from g by the application of an integral operator. The estimation of g then entails both statistical and numerical difficulties. A natural statistical approach is by maximum likelihood, conveniently implemented using the EM algorithm, but this provides unsatisfactory reconstructions of g. In this paper, we modify the maximum likelihood-EM approach by introducing a simple smoothing step at each EM iteration. In our experience, this algorithm converges in relatively few iterations to good estimates of g that do not depend on the choice of starting configuration. Some theoretical background is given that relates this smoothed EM algorithm to a maximum penalized likelihood approach. Two applications are considered in detail. The first is the classical stereology problem of determining particle size distributions from data collected on a plane section through a composite medium. The second concerns the recovery of the structure of a section of the human body from external observations obtained by positron emission tomography; for this problem, we also suggest several technical improvements on existing methodology.

https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1990.tb01788.x

A smoothed EM approach to indirect estimation problems, with particular reference to stereology and emission tomography

Extending the work of A.W. McCarthy et al. (1988) and M.I. Miller and B. Roysam (1991), the authors demonstrate that a fully parallel implementation of the maximum-likelihood method for single-photon emission computed tomography (SPECT) can be accomplished in clinical time frames on massively parallel systolic array processors. The authors show that for SPECT imaging on 64*64 image grids, with 96 view angles, the single-instruction, multiple data (SIMD) distributed array processor containing 64/sup 2/ processors performs the expectation-maximization (EM) algorithm with Good's smoothing at a rate of 1 iteration/1.5 s. This promises for emission tomography fully Bayesian reconstructions including regularization in clinical computation times which are on the order of 1 min/slice. The most important result of the implementations is that the scaling rules for computation times are roughly linear in the number of processors. >

Maximum likelihood SPECT in clinical computation times using mesh-connected parallel computers

A fully three-dimensional (3-D) implementation of the maximum a posteriori (MAP) method for single photon emission computed tomography (SPECT) is demonstrated. The 3-D reconstruction exhibits a major increase in resolution when compared to the generation of the series of separate 2-D slice reconstructions. As has been noted, the iterative EM algorithm for 2-D reconstruction is highly computational; the 3-D algorithm is far worse. To accommodate the computational complexity, previous work in the 2-D arena is extended, and an implementation on the class of massively parallel processors of the 3-D algorithm is demonstrated. Using a 16000- (4000-) processor MasPar/DECmpp-Sx machine, the algorithm is demonstrated to execute at 2.5

3-D maximum a posteriori estimation for single photon emission computed tomography on massively-parallel computers

Single photon emission computed tomography (SPECT) reconstructions performed using maximum a posteriori (penalized likelihood) estimation with the expectation maximization algorithm are discussed. Due to the large number of computations, the algorithms were performed on a massively parallel single-instruction multiple-data computer. Computation times for 200 iterations, using I.J. Good and R.A. Gaskins's (1971) roughness as a rotationally invariant roughness penalty, are shown to be on the order of 5 min for a 64*64 image with 96 view angles on an AMT-DAP 4096 processor machine and 1 min on a MasPar 4096 processor machine. Computer simulations performed using parameters for the Siemens gamma camera and clinical brain scan parameters are presented to compare two regularization techniques-regularization by kernel sieves and penalized likelihood with Good's rotationally invariant roughness measure-to filtered backprojection. Twenty-five independent sets of data are reconstructed for the pie and Hoffman brain phantoms. The average variance and average deviation are examined in various areas of the brain phantom. It is shown that while the geometry of the area examined greatly affects the observed results, in all cases the reconstructions using Good's roughness give superior variance and bias results to the two alternative methods. >

Maximum a posteriori estimation for SPECT using regularization techniques on massively parallel computers

Single photon emission computed tomography (SPECT) reconstructions were performed using maximum a posteriori (penalized likelihood) estimation via the expectation maximization algorithm Due to the large number of computations, the algorithms were performed on a massively parallel single-instruction multiple-data computer Computation times for 200 iterations using Good's rotationally invariant roughness penalty were on the order of 5 min for a 64*64 image with 96 view angles on a 4096 processor machine, and 40 s on a 16384 processor machine Computer simulations were performed using parameters for the Siemens gamma camera and clinical brain scan parameters comparing two regularization techniques to conventional reconstructions Regularization by kernel sieves and penalized likelihood with Good's rotationally invariant roughness measure are compared to filtered backprojection Twenty-five independent sets of data are reconstructed for the pie and Hoffman brain phantoms >

https://ieeexplore.ieee.org/iel2/623/6551/00259267.pdf

Maximum a posteriori estimation for SPECT using regularization techniques on massively-parallel computers

A constrained maximum-likelihood estimator is derived by incorporating a rotationally invariant roughness penalty proposed by I.J. Good (1981) into the likelihood functional. This leads to a set of nonlinear differential equations the solution of which is a spline-smoothing of the data. The nonlinear partial differential equations are mapped onto a grid via finite differences, and it is shown that the resulting computations possess a high degree of parallelism as well as locality in the data-passage, which allows an efficient implementation on a 48-by-48 mesh-connected array of NCR GAPP processors. The smooth reconstruction of the intensity functions of Poisson point processes is demonstrated in two dimensions. >

J.A. Shrauner

Papers

Bayesian imaging using Good's roughness measure-implementation on a massively parallel processor