Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): applications to tomography

Theoretical study of lesion detectability of MAP reconstruction using computer observers

The low signal-to-noise ratio (SNR) in emission data has stimulated the development of statistical image reconstruction methods based on the maximum a posteriori (MAP) principle. Experimental examples have shown that statistical methods improve image quality compared to the conventional filtered backprojection (FBP) method. However, these results depend on isolated data sets. Here we study the lesion detectability of MAP reconstruction theoretically, using computer observers. These theoretical results can be applied to different object structures. They show that for a quadratic smoothing prior, the lesion detectability using the prewhitening observer is independent of the smoothing parameter and the neighborhood of the prior, while the nonprewhitening observer exhibits an optimum smoothing point. We also compare the results to those of FBP reconstruction. The comparison shows that for ideal positron emission tomography (PET) systems (where data are true line integrals of the tracer distribution) the MAP reconstruction has a higher SNR for lesion detection than FBP reconstruction due to the modeling of the Poisson noise. For realistic systems, MAP reconstruction further benefits from accurately modeling the physical photon detection process in PET.

Maximum-likelihood transmission image reconstruction for overlapping transmission beams

In many transmission imaging geometries, the transmitted "beams" of photons overlap on the detector, such that a detector element may record photons that originated in different sources or source locations and thus traversed different paths through the object. Examples include systems based on scanning line sources or on multiple parallel rod sources. The overlap of these beams has been disregarded by both conventional analytical reconstruction methods as well as by previous statistical reconstruction methods. We propose a new algorithm for statistical image reconstruction of attenuation maps that explicitly accounts for overlapping beams in transmission scans. The algorithm is guaranteed to monotonically increase the objective function at each iteration. The availability of this algorithm enables the possibility of deliberately increasing the beam overlap so as to increase count rates. Simulated single photon emission tomography transmission scans based on a multiple line source array demonstrate that the proposed method yields improved resolution/noise tradeoffs relative to "conventional" reconstruction algorithms, both statistical and nonstatistical.

Nonparametric regression sinogram smoothing using a roughness-penalized Poisson likelihood objective function

We develop and investigate an approach to tomographic image reconstruction in which nonparametric regression using a roughness-penalized Poisson likelihood objective function is used to smooth each projection independently prior to reconstruction by unapodized filtered backprojection (FBP). As an added generalization, the roughness penalty is expressed in terms of a monotonic transform, known as the link function, of the projections. The approach is compared to shift-invariant projection filtering through the use of a Hanning window as well as to a related nonparametric regression approach that makes use of an objective function based on weighted least squares (WLS) rather than the Poisson likelihood. The approach is found to lead to improvements in resolution-noise tradeoffs over the Hanning filter as well as over the WLS approach. We also investigate the resolution and noise effects of three different link functions: the identity, square root, and logarithm links. The choice of link function is found to influence the resolution uniformity and isotropy properties of the reconstructed images. In particular, in the case of an idealized imaging system with intrinsically uniform and isotropic resolution, the choice of a square root link function yields the desirable outcome of essentially uniform and isotropic resolution in reconstructed images, with noise performance still superior to that of the Hanning filter as well as that of the WLS approach.

Resolution and noise properties of MAP reconstruction for fully 3-D PET

We derive approximate analytical expressions for the local impulse response and covariance of images reconstructed from fully three-dimensional (3-D) positron emission tomography (PET) data using maximum a posteriori (MAP) estimation. These expressions explicitly account for the spatially variant detector response and sensitivity of a 3-D tomograph. The resulting spatially variant impulse response and covariance are computed using 3-D Fourier transforms. A truncated Gaussian distribution is used to account for the effect on the variance of the nonnegativity constraint used in MAP reconstruction. Using Monte Carlo simulations and phantom data from the microPET small animal scanner, we show that the approximations provide reasonably accurate estimates of contrast recovery and covariance of MAP reconstruction for priors with quadratic energy functions. We also describe how these analytical results can be used to achieve near-uniform contrast recovery throughout the reconstructed volume.

Noise characterization of block-iterative reconstruction algorithms: I. Theory

Researchers have shown increasing interest in block-iterative image reconstruction algorithms due to the computational and modeling advantages they provide. Although their convergence properties have been well documented, little is known about how they behave in the presence of noise. In this work, we fully characterize the ensemble statistical properties of the rescaled block-iterative expectation-maximization (RBI-EM) reconstruction algorithm and the rescaled block-iterative simultaneous multiplicative algebraic reconstruction technique (RBI-SMART). Also included in the analysis are the special cases of RBI-EM, maximum-likelihood EM (ML-EM) and ordered-subset EM (OS-EM), and the special case of RBI-SMART, SMART. A theoretical formulation strategy similar to that previously outlined for ML-EM is followed for the RBI methods. The theoretical formulations in this paper rely on one approximation, namely, that the noise in the reconstructed image is small compared to the mean image. In a second paper, the approximation will be justified through Monte Carlo simulations covering a range of noise levels, iteration points, and subset orderings. The ensemble statistical parameters could then be used to evaluate objective measures of image quality.

Exact distribution of edge-preserving MAP estimators for linear signal models with Gaussian measurement noise

We derive the exact statistical distribution of maximum a posteriori (MAP) estimators having edge-preserving nonGaussian priors. Such estimators have been widely advocated for image restoration and reconstruction problems. Previous investigations of these image recovery methods have been primarily empirical; the distribution we derive enables theoretical analysis. The signal model is linear with Gaussian measurement noise. We assume that the energy function of the prior distribution is chosen to ensure a unimodal posterior distribution (for which convexity of the energy function is sufficient), and that the energy function satisfies a uniform Lipschitz regularity condition. The regularity conditions are sufficiently general to encompass popular priors such as the generalized Gaussian Markov random field prior and the Huber prior, even though those priors are not everywhere twice continuously differentiable.

Noise propagation in wave-front sensing with phase diversity

The phase diversity technique is studied as a wave-front sensor to be implemented with widely extended sources. The wave-front phase expanded on the Zernike polynomials is estimated from a pair of images (in focus and out of focus) by use of a maximum-likelihood approach. The propagation of the photon noise in the images on the estimated phase is derived from a theoretical analysis. The covariance matrix of the phase estimator is calculated, and the optimal distance between the observation planes that minimizes the noise propagation is determined. The phase error is inversely proportional to the number of photons in the images. The noise variance on the Zernike polynomials increases with the order of the polynomial. These results are confirmed with both numerical and experimental validations. The influence of the spectral bandwidth on the phase estimator is also studied with simulations.

Penalized-likelihood estimators and noise analysis for randoms-precorrected PET transmission scans

This paper analyzes and compares image reconstruction methods based on practical approximations to the exact log likelihood of randoms-precorrected positron emission tomography (PET) measurements. The methods apply to both emission and transmission tomography, however, in this paper we focus on transmission tomography. The results of experimental PET transmission scans and variance approximations demonstrate that the shifted Poisson (SP) method avoids the systematic bias of the conventional data-weighted least squares (WLS) method and leads to significantly lower variance than conventional statistical methods based on the log likelihood of the ordinary Poisson (OP) model. We develop covariance approximations to analyze the propagation of noise from attenuation maps into emission images via the attenuation correction factors (ACF's). Empirical pixel and region variances from real transmission data agree closely with the analytical predictions. Both the approximations and the empirical results show that the performance differences between the OP model and SP model are even larger, when considering noise propagation from the transmission images into the final emission images, than the differences in the attenuation maps themselves.

A theoretical study of the contrast recovery and variance of MAP reconstructions from PET data

We examine the spatial resolution and variance properties of PET images reconstructed using maximum a posteriori (MAP) or penalized-likelihood methods. Resolution is characterized by the contrast recovery coefficient (CRC) of the local impulse response. Simplified approximate expressions are derived for the local impulse response CRC's and variances for each voxel. Using these results we propose a practical scheme for selecting spatially variant smoothing parameters to optimize lesion detectability through maximization of the local CRC-to-noise ratio in the reconstructed image.

Conjugate-gradient preconditioning methods for shift-variant PET image reconstruction

Gradient-based iterative methods often converge slowly for tomographic image reconstruction and image restoration problems, but can be accelerated by suitable preconditioners. Diagonal preconditioners offer some improvement in convergence rate, but do not incorporate the structure of the Hessian matrices in imaging problems. Circulant preconditioners can provide remarkable acceleration for inverse problems that are approximately shift-invariant, i.e., for those with approximately block-Toeplitz or block-circulant Hessians. However, in applications with nonuniform noise variance, such as arises from Poisson statistics in emission tomography and in quantum-limited optical imaging, the Hessian of the weighted least-squares objective function is quite shift-variant, and circulant preconditioners perform poorly. Additional shift-variance is caused by edge-preserving regularization methods based on nonquadratic penalty functions. This paper describes new preconditioners that approximate more accurately the Hessian matrices of shift-variant imaging problems. Compared to diagonal or circulant preconditioning, the new preconditioners lead to significantly faster convergence rates for the unconstrained conjugate-gradient (CG) iteration. We also propose a new efficient method for the line-search step required by CG methods. Applications to positron emission tomography (PET) illustrate the method.

Errors in biased estimators for parametric ultrasonic imaging

Maximum likelihood (ML) methods are widely used in acoustic parameter estimation. Although ML methods are often unbiased, the variance is unacceptably large for many applications, including medical imaging. For such cases, Bayesian estimators can reduce variance and preserve contrast at the cost of an increased bias. Consequently, including prior knowledge about object and noise properties in the estimator can improve low-contrast target detectability of parametric ultrasound images by improving the precision of the estimates. In this paper, errors introduced by biased estimators are analyzed and approximate closed-form expressions are developed. The task-specific nature of the estimator performance is demonstrated through analysis, simulation, and experimentation. A strategy for selecting object priors is proposed. Acoustic scattering from kidney tissue is the emphasis of this paper, although the results are more generally applicable.

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.

Grouped-coordinate ascent algorithms for penalized-likelihood transmission image reconstruction

This paper presents a new class of algorithms for penalized-likelihood reconstruction of attenuation maps from low-count transmission scans. We derive the algorithms by applying to the transmission log-likelihood a version of the convexity technique developed by De Pierro for emission tomography. The new class includes the single-coordinate ascent (SCA) algorithm and Lange's convex algorithm for transmission tomography as special cases. The new grouped-coordinate ascent (GCA) algorithms in the class overcome several limitations associated with previous algorithms. 1) Fewer exponentiations are required than in the transmission maximum likelihood-expectation maximization (ML-EM) algorithm or in the SCA algorithm. 2) The algorithms intrinsically accommodate nonnegativity constraints, unlike many gradient-based methods. 3) The algorithms are easily parallelizable, unlike the SCA algorithm and perhaps line-search algorithms. We show that the GCA algorithms converge faster than the SCA algorithm, even on conventional workstations. An example from a low-count positron emission tomography (PET) transmission scan illustrates the method.

Spatial resolution properties of penalized-likelihood image reconstruction: space-invariant tomographs

This paper examines the spatial resolution properties of penalized-likelihood image reconstruction methods by analyzing the local impulse response. The analysis shows that standard regularization penalties induce space-variant local impulse response functions, even for space-invariant tomographic systems. Paradoxically, for emission image reconstruction, the local resolution is generally poorest in high-count regions. We show that the linearized local impulse response induced by quadratic roughness penalties depends on the object only through its projections. This analysis leads naturally to a modified regularization penalty that yields reconstructed images with nearly uniform resolution. The modified penalty also provides a very practical method for choosing the regularization parameter to obtain a specified resolution in images reconstructed by penalized-likelihood methods.