Monotonic algorithms for transmission tomography

Statistical projection completion in X-ray CT using consistency conditions

Projection data incompleteness arises in many situations relevant to X-ray computed tomography (CT) imaging. We propose a penalized maximum likelihood statistical sinogram restoration approach that incorporates the Helgason-Ludwig (HL) consistency conditions to accommodate projection data incompleteness. Image reconstruction is performed by the filtered-backprojection (FBP) in a second step. In our problem formulation, the objective function consists of the log-likelihood of the X-ray CT data and a penalty term; the HL condition poses a linear constraint on the restored sinogram and can be implemented efficiently via fast Fourier transform (FFT) and inverse FFT. We derive an iterative algorithm that increases the objective function monotonically. The proposed algorithm is applied to both computer simulated data and real patient data. We study different factors in the problem formulation that affect the properties of the final FBP reconstructed images, including the data truncation level, the amount of prior knowledge on the object support, as well as different approximations of the statistical distribution of the available projection data. We also compare its performance with an analytical truncation artifacts reduction method. The proposed method greatly improves both the accuracy and the precision of the reconstructed images within the scan field-of-view, and to a certain extent recovers the truncated peripheral region of the object. The proposed method may also be applied in areas such as limited angle tomography, metal artifacts reduction, and sparse sampling imaging.

Generalized algorithms for direct reconstruction of parametric images from dynamic PET data

Indirect and direct methods have been developed for reconstructing parametric images from dynamic positron emission tomography (PET) data. Indirect methods are simple and easy to implement because reconstruction and kinetic modeling are performed in two separate steps. Direct methods estimate parametric images directly from dynamic PET sinograms and, in theory, can be statistically more efficient, but the algorithms are often difficult to implement and are very specific to the kinetic model being used. This paper presents a class of generalized algorithms for direct reconstruction of parametric images that are relatively easy to implement and can be adapted to different kinetic models. The proposed algorithms use optimization transfer principle to convert the maximization of a penalized likelihood into a pixel-wise weighted least squares (WLS) kinetic fitting problem at each iteration. Thus, it can employ existing WLS algorithms developed for kinetic models. The proposed algorithms resemble the empirical iterative implementation of the indirect approach, but converge to a solution of the direct formulation. Computer simulations showed that the proposed direct reconstruction algorithms are flexible and achieve a better bias-variance tradeoff than indirect reconstruction methods.

PET transmission tomography using a novel nonlocal MRF prior

In positron emission tomography, transmission scans can be performed to estimate attenuation correction factors (ACFs) which are in turn used to correct the emission scans. And such an attenuation correction is crucial for quantitatively accurate PET reconstructions. The prior model used in this work was based on our assumption that the attenuation values vary smoothly, with occasional discontinuities at anatomical borders. And on the other hand, long acquisition or scan times, although alleviating the noise effect of the count-limited scans, are blamed for patient uncomfortableness and movements. So, transmission tomography often suffers from the noise effect because of the short scan time. Thus reconstruction which is capable of overcoming the noise effect is highly needed. In this article, we apply the nonlocal prior Bayesian reconstruction method in PET transmission tomography. Resulting experimentations validate that the reconstructions using the nonlocal prior can reconstruct better transmission images and overcome noise effect even when the scan time is relatively short.

Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior

How to reduce the radiation dose delivered to the patients has always been a important concern since the introduction of computed tomography (CT). Though clinically desired, low-dose CT images can be severely degraded by the excessive quantum noise under extremely low X-ray dose circumstances. Bayesian statistical reconstructions outperform the traditional filtered back-projection (FBP) reconstructions by accurately expressing the system models of physical effects and the statistical character of the measurement data. This work aims to improve the image quality of low-dose CT images using a novel AW nonlocal (adaptive-weighting nonlocal) prior statistical reconstruction approach. Compared to traditional local priors, the proposed prior can adaptively and selectively exploit the global image information. It imposes an effective resolution-preserving and noise-removing regularization for reconstructions. Experimentation validates that the reconstructions using the proposed prior have excellent performance for X-ray CT with low-dose scans.

Electronic noise modeling in statistical iterative reconstruction

We consider electronic noise modeling in tomographic image reconstruction when the measured signal is the sum of a Gaussian distributed electronic noise component and another random variable whose log-likelihood function satisfies a certain linearity condition. Examples of such likelihood functions include the Poisson distribution and an exponential dispersion (ED) model that can approximate the signal statistics in integration mode X-ray detectors. We formulate the image reconstruction problem as a maximum-likelihood estimation problem. Using an expectation-maximization approach, we demonstrate that a reconstruction algorithm can be obtained following a simple substitution rule from the one previously derived without electronic noise considerations. To illustrate the applicability of the substitution rule, we present examples of a fully iterative reconstruction algorithm and a sinogram smoothing algorithm both in transmission CT reconstruction when the measured signal contains additive electronic noise. Our simulation studies show the potential usefulness of accurate electronic noise modeling in low-dose CT applications.

Development of novel statistical reconstruction algorithms for poly-energetic X-ray computed tomography

A beam-hardening effect is a common problem affecting the quantitative aspects of X-ray computed tomography (CT). We have developed two statistical reconstruction algorithms for poly-energetic X-ray CT that can effectively reduce the beam-hardening effect. Phantom tests were used to evaluate our approach in comparison with traditional correction methods. Unlike previous methods, our algorithm utilizes multiple energy-corresponding blank scans to estimate the attenuation map for a particular energy spectrum. Therefore, our algorithm is an energy-selective reconstruction. In addition to benefits over other statistical algorithms for poly-energetic reconstruction, our algorithm has the advantage of not requiring prior knowledge of the object material, the energy spectrum of the source and the energy sensitivity of the detector. The results showed an improvement in coefficient of variation, uniformity and signal-to-noise ratio; overall, this novel approach produces a better beam-hardening correction.

Analysis of Resolution and Noise Properties of Nonquadratically Regularized Image Reconstruction Methods for PET

We present accurate and efficient methods for estimating the spatial resolution and noise properties of nonquadratically regularized image reconstruction for positron emission tomography (PET). It is well known that quadratic regularization tends to over-smooth sharp edges. Many types of edge-preserving nonquadratic penalties have been proposed to overcome this problem. However, there has been little research on the quantitative analysis of nonquadratic regularization due to its nonlinearity. In contrast, quadratically regularized estimators are approximately linear and are well understood in terms of resolution and variance properties. We derive new approximate expressions for the linearized local perturbation response (LLPR) and variance using the Taylor expansion with the remainder term. Although the expressions are implicit, we can use them to accurately predict resolution and variance for nonquadratic regularization where the conventional expressions based on the first-order Taylor truncation fail. They also motivate us to extend the use of a certainty-based modified penalty to nonquadratic regularization cases in order to achieve spatially uniform perturbation responses, analogous to uniform spatial resolution in quadratic regularization. Finally, we develop computationally efficient methods for predicting resolution and variance of nonquadratically regularized reconstruction and present simulations that illustrate the validity of these methods.

Majorization-minimization algorithms for wavelet-based image restoration

Standard formulations of image/signal deconvolution under wavelet-based priors/regularizers lead to very high-dimensional optimization problems involving the following difficulties: the non-Gaussian (heavy-tailed) wavelet priors lead to objective functions which are nonquadratic, usually nondifferentiable, and sometimes even nonconvex; the presence of the convolution operator destroys the separability which underlies the simplicity of wavelet-based denoising. This paper presents a unified view of several recently proposed algorithms for handling this class of optimization problems, placing them in a common majorization-minimization (MM) framework. One of the classes of algorithms considered (when using quadratic bounds on nondifferentiable log-priors) shares the infamous "singularity issue" (SI) of "iteratively reweighted least squares" (IRLS) algorithms: the possibility of having to handle infinite weights, which may cause both numerical and convergence issues. In this paper, we prove several new results which strongly support the claim that the SI does not compromise the usefulness of this class of algorithms. Exploiting the unified MM perspective, we introduce a new algorithm, resulting from using l1 bounds for nonconvex regularizers; the experiments confirm the superior performance of this method, when compared to the one based on quadratic majorization. Finally, an experimental comparison of the several algorithms, reveals their relative merits for different standard types of scenarios.

An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms

The majorize-minimize (MM) optimization technique has received considerable attention in signal and image processing applications, as well as in statistics literature. At each iteration of an MM algorithm, one constructs a tangent majorant function that majorizes the given cost function and is equal to it at the current iterate. The next iterate is obtained by minimizing this tangent majorant function, resulting in a sequence of iterates that reduces the cost function monotonically. A well-known special case of MM methods are expectation-maximization algorithms. In this paper, we expand on previous analyses of MM, due to Fessler and Hero, that allowed the tangent majorants to be constructed in iteration-dependent ways. Also, this paper overcomes an error in one of those earlier analyses. There are three main aspects in which our analysis builds upon previous work. First, our treatment relaxes many assumptions related to the structure of the cost function, feasible set, and tangent majorants. For example, the cost function can be nonconvex and the feasible set for the problem can be any convex set. Second, we propose convergence conditions, based on upper curvature bounds, that can be easier to verify than more standard continuity conditions. Furthermore, these conditions allow for considerable design freedom in the iteration-dependent behavior of the algorithm. Finally, we give an original characterization of the local region of convergence of MM algorithms based on connected (e.g., convex) tangent majorants. For such algorithms, cost function minimizers will locally attract the iterates over larger neighborhoods than typically is guaranteed with other methods. This expanded treatment widens the scope of the MM algorithm designs that can be considered for signal and image processing applications, allows us to verify the convergent behavior of previously published algorithms, and gives a fuller understanding overall of how these algorithms behave.

SAR image autofocus by sharpness optimization: a theoretical study

Synthetic aperture radar (SAR) autofocus techniques that optimize sharpness metrics can produce excellent restorations in comparison with conventional autofocus approaches. To help formalize the understanding of metric-based SAR autofocus methods, and to gain more insight into their performance, we present a theoretical analysis of these techniques using simple image models. Specifically, we consider the intensity-squared metric, and a dominant point-targets image model, and derive expressions for the resulting objective function. We examine the conditions under which the perfectly focused image models correspond to stationary points of the objective function. A key contribution is that we demonstrate formally, for the specific case of intensity-squared minimization autofocus, the mechanism by which metric-based methods utilize the multichannel defocusing model of SAR autofocus to enforce the stationary point property for multiple image columns. Furthermore, our analysis shows that the objective function has a special separble property through which it can be well approximated locally by a sum of 1-D functions of each phase error component. This allows fast performance through solving a sequence of 1-D optimization problems for each phase component simultaneously. Simulation results using the proposed models and actual SAR imagery confirm that the analysis extends well to realistic situations.

Alternating minimization algorithms for transmission tomography

A family of alternating minimization algorithms for finding maximum-likelihood estimates of attenuation functions in transmission X-ray tomography is described. The model from which the algorithms are derived includes polyenergetic photon spectra, background events, and nonideal point spread functions. The maximum-likelihood image reconstruction problem is reformulated as a double minimization of the I-divergence. A novel application of the convex decomposition lemma results in an alternating minimization algorithm that monotonically decreases the objective function. Each step of the minimization is in closed form. The family of algorithms includes variations that use ordered subset techniques for increasing the speed of convergence. Simulations demonstrate the ability to correct the cupping artifact due to beam hardening and the ability to reduce streaking artifacts that arise from beam hardening and background events.

Estimating 3-D respiratory motion from orbiting views by tomographic image registration

Respiratory motion remains a significant source of errors in treatment planning for the thorax and upper abdomen. Recently, we proposed a method to estimate two-dimensional (2-D) object motion from a sequence of slowly rotating X-ray projection views, which we called deformation from orbiting views (DOVs). In this method, we model the motion as a time varying deformation of a static prior of the anatomy. We then optimize the parameters of the motion model by maximizing the similarity between the modeled and actual projection views. This paper extends the method to full three-dimensional (3-D) motion and cone-beam projection views. We address several practical issues for using a cone-beam computed tomography (CBCT) scanner that is integrated in a radiotherapy system, such as the effects of Compton scatter and the limited gantry rotation for one breathing cycle. We also present simulation and phantom results to illustrate the performance of this method.

Pose estimation of known objects during transmission tomographic image reconstruction

We address the problem of image formation in transmission tomography when metal objects of known composition and shape, but unknown pose, are present in the scan subject. Using an alternating minimization (AM) algorithm, derived from a model in which the detected data are viewed as Poisson-distributed photon counts, we seek to eliminate the streaking artifacts commonly seen in filtered back projection images containing high-contrast objects. We show that this algorithm, which minimizes the I-divergence (or equivalently, maximizes the log-likelihood) between the measured data and model-based estimates of the means of the data, converges much faster when knowledge of the high-density materials (such as brachytherapy applicators or prosthetic implants) is exploited. The algorithm incorporates a steepest descent-based method to find the position and orientation (collectively called the pose) of the known objects. This pose is then used to constrain the image pixels to their known attenuation values, or, for example, to form a mask on the "missing" projection data in the shadow of the objects. Results from two-dimensional simulations are shown in this paper. The extension of the model and methods used to three dimensions is outlined.

Monotonic penalized-likelihood image reconstruction for X-ray fluorescence computed tomography

In this paper, we derive a monotonic penalized-likelihood algorithm for image reconstruction in X-ray fluorescence computed tomography (XFCT) when the attenuation maps at the energies of the fluorescence X-rays are unknown. In XFCT, a sample is irradiated with pencil beams of monochromatic synchrotron radiation that stimulate the emission of fluorescence X-rays from atoms of elements whose K- or L-edges lie below the energy of the stimulating beam. Scanning and rotating the object through the beam allows for acquisition of a tomographic dataset that can be used to reconstruct images of the distribution of the elements in question. XFCT is a stimulated emission tomography modality, and it is thus necessary to correct for attenuation of the incident and fluorescence photons. The attenuation map is, however, generally known only at the stimulating beam energy and not at the energies of the various fluorescence X-rays of interest. We have developed a penalized-likelihood image reconstruction strategy for this problem. The approach alternates between updating the distribution of a given element and updating the attenuation map for that element's fluorescence X-rays. The approach is guaranteed to increase the penalized likelihood at each iteration. Because the joint objective function is not necessarily concave, the approach may drive the solution to a local maximum. To encourage the algorithm to seek out a reasonable local maximum, we include in the objective function a prior that encourages a relationship, based on physical considerations, between the fluorescence attenuation map and the distribution of the element being reconstructed.

Penalized-likelihood sinogram restoration for computed tomography

We formulate computed tomography (CT) sinogram preprocessing as a statistical restoration problem in which the goal is to obtain the best estimate of the line integrals needed for reconstruction from the set of noisy, degraded measurements. CT measurement data are degraded by a number of factors-including beam hardening and off-focal radiation-that produce artifacts in reconstructed images unless properly corrected. Currently, such effects are addressed by a sequence of sinogram-preprocessing steps, including deconvolution corrections for off-focal radiation, that have the potential to amplify noise. Noise itself is generally mitigated through apodization of the reconstruction kernel, which effectively ignores the measurement statistics, although in high-noise situations adaptive filtering methods that loosely model data statistics are sometimes applied. As an alternative, we present a general imaging model relating the degraded measurements to the sinogram of ideal line integrals and propose to estimate these line integrals by iteratively optimizing a statistically based objective function. We consider three different strategies for estimating the set of ideal line integrals, one based on direct estimation of ideal "monochromatic" line integrals that have been corrected for single-material beam hardening, one based on estimation of ideal "polychromatic" line integrals that can be readily mapped to monochromatic line integrals, and one based on estimation of ideal transmitted intensities, from which ideal, monochromatic line integrals can be readily estimated. The first two approaches involve maximization of a penalized Poisson-likelihood objective function while the third involves minimization of a quadratic penalized weighted least squares (PWLS) objective applied in the transmitted intensity domain. We find that at low exposure levels typical of those being considered for screening CT, the Poisson-likelihood based approaches outperform the PWLS objective as well as a standard approach based on adaptive filtering followed by deconvolution. At higher exposure levels, the approaches all perform similarly.

Incorporation of system resolution compensation (RC) in the ordered-subset transmission (OSTR) algorithm for transmission imaging in SPECT

In order to reconstruct attenuation maps with improved spatial resolution and quantitative accuracy, we developed an approximate method of incorporating system resolution compensation (RC) in the ordered-subset transmission (OSTR) algorithm for transmission reconstruction. Our method approximately models the blur caused by the finite intrinsic detector resolution, the nonideal source collimation and detector collimation. We derived the formulation using the optimization transfer principle as in the derivation of the OSTR algorithm. The formulation includes one forward-blur step and one back-blur step, which do not severely slow down reconstruction. The formulation could be applicable to various transmission geometries, such as point-source, line-source, and sheet-source systems. Through computer simulations of the MCAT phantom and transmission measurements of the air-filled Data Spectrum Deluxe single photo emission computed tomography (SPECT) Phantom on a system which employed a cone-beam geometry and a system which employed a scanning-line-source geometry, we showed that incorporation of RC increased spatial resolution and improved the quantitative accuracy of reconstruction. In simulation studies, attenuation maps reconstructed with RC correction improved the quantitative accuracy of emission reconstruction.

On global and local convergence of half-quadratic algorithms

This paper provides original results on the global and local convergence properties of half-quadratic (HQ) algorithms resulting from the Geman and Yang (GY) and Geman and Reynolds (GR) primal-dual constructions. First, we show that the convergence domain of the GY algorithm can be extended with the benefit of an improved convergence rate. Second, we provide a precise comparison of the convergence rates for both algorithms. This analysis shows that the GR form does not benefit from a better convergence rate in general. Moreover, the GY iterates often take advantage of a low cost implementation. In this case, the GY form is usually faster than the GR form from the CPU time viewpoint.

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.

Comparison of quadratic- and median-based roughness penalties for penalized-likelihood sinogram restoration in computed tomography

We have compared the performance of two different penalty choices for a penalized-likelihood sinogram-restoration strategy we have been developing. One is a quadratic penalty we have employed previously and the other is a new median-based penalty. We compared the approaches to a noniterative adaptive filter that loosely but not explicitly models data statistics. We found that the two approaches produced similar resolution-variance tradeoffs to each other and that they outperformed the adaptive filter in the low-dose regime, which suggests that the particular choice of penalty in our approach may be less important than the fact that we are explicitly modeling data statistics at all. Since the quadratic penalty allows for derivation of an algorithm that is guaranteed to monotonically increase the penalized-likelihood objective function, we find it to be preferable to the median-based penalty.

A new simple iterative reconstruction algorithm for SPECT transmission measurement

This paper proposes a new iterative reconstruction algorithm for transmission tomography and compares this algorithm with several other methods. The new algorithm is simple and resembles the emission ML-EM algorithm in form. Due to its simplicity, it is easy to implement and fast to compute a new update at each iteration. The algorithm also always guarantees non-negative solutions. Evaluations are performed using simulation studies and real phantom data. Comparisons with other algorithms such as convex, gradient, and logMLEM show that the proposed algorithm is as good as others and performs better in some cases.

Penalized-likelihood sinogram smoothing for low-dose CT

We have developed a sinogram smoothing approach for low-dose computed tomography (CT) that seeks to estimate the line integrals needed for reconstruction from the noisy measurements by maximizing a penalized-likelihood objective function. The maximization is performed by an algorithm derived by use of the separable paraboloidal surrogates framework. The approach overcomes some of the computational limitations of a previously proposed spline-based penalized-likelihood sinogram smoothing approach, and it is found to yield better resolution-variance tradeoffs than this spline-based approach as well an existing adaptive filtering approach. Such sinogram smoothing approaches could be valuable when applied to the low-dose data acquired in CT screening exams, such as those being considered for lung-nodule detection.

Accelerated penalized weighted least-squares and maximum likelihood algorithms for reconstructing transmission images from PET transmission data

We present penalized weighted least-squares (PWLS) and penalized maximum-likelihood (PML) methods for reconstructing transmission images from positron emission tomography transmission data. First, we view the problem of minimizing the weighted least-squares (WLS) and maximum likelihood objective functions as a sequence of nonnegative least-squares minimization problems. This viewpoint follows from using certain quadratic functions as surrogate functions for the WLS and maximum likelihood objective functions. Second, we construct surrogate functions for a class of penalty functions that yield closed form expressions for the iterates of the PWLS and PML algorithms. Due to the slow convergence of the PWLS and PML algorithms, accelerated versions of them are developed that are theoretically guaranteed to monotonically decrease their respective objective functions. In experiments using real phantom data, the PML images produced the most accurate attenuation correction factors. On the other hand, the PWLS images produced images with the highest levels of contrast for low-count data.

Evaluation of an efficient compensation method for quantitative fan-beam brain SPECT reconstruction

Fan-beam collimators are designed to improve the system sensitivity and resolution for imaging small objects such as the human brain and breasts in single photon emission computed tomography (SPECT). Many reconstruction algorithms have been studied and applied to this geometry to deal with every kind of degradation factor. This paper presents a new reconstruction approach for SPECT with circular orbit, which demonstrated good performance in terms of both accuracy and efficiency. The new approach compensates for degradation factors including noise, scatter, attenuation, and spatially variant detector response. Its uniform attenuation approximation strategy avoids the additional transmission scan for the brain attenuation map, hence reducing the patient radiation dose and furthermore simplifying the imaging procedure. We evaluate and compare this new approach with the well-established ordered-subset expectation-maximization iterative method, using Monte Carlo simulations. We perform quantitative analysis with regional bias-variance, receiver operating characteristics, and Hotelling trace studies for both methods. The results demonstrate that our reconstruction strategy has comparable performance with a significant reduction of computing time.

Efficient calculation of resolution and covariance for penalized-likelihood reconstruction in fully 3-D SPECT

Resolution and covariance predictors have been derived previously for penalized-likelihood estimators. These predictors can provide accurate approximations to the local resolution properties and covariance functions for tomographic systems given a good estimate of the mean measurements. Although these predictors may be evaluated iteratively, circulant approximations are often made for practical computation times. However, when numerous evaluations are made repeatedly (as in penalty design or calculation of variance images), these predictors still require large amounts of computing time. In Stayman and Fessler (2000), we discussed methods for precomputing a large portion of the predictor for shift-invariant system geometries. In this paper, we generalize the efficient procedure discussed in Stayman and Fessler (2000) to shift-variant single photon emission computed tomography (SPECT) systems. This generalization relies on a new attenuation approximation and several observations on the symmetries in SPECT systems. These new general procedures apply to both two-dimensional and fully three-dimensional (3-D) SPECT models, that may be either precomputed and stored, or written in procedural form. We demonstrate the high accuracy of the predictions based on these methods using a simulated anthropomorphic phantom and fully 3-D SPECT system. The evaluation of these predictors requires significantly less computation time than traditional prediction techniques, once the system geometry specific precomputations have been made.

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.

Emission image reconstruction for randoms-precorrected PET allowing negative sinogram values

Most positron emission tomography (PET) emission scans are corrected for accidental coincidence (AC) events by real-time subtraction of delayed-window coincidences, leaving only the randoms-precorrected data available for image reconstruction. The real-time randoms precorrection compensates in mean for AC events but destroys the Poisson statistics. The exact log-likelihood for randoms-precorrected data is inconvenient, so practical approximations are needed for maximum likelihood or penalized-likelihood image reconstruction. Conventional approximations involve setting negative sinogram values to zero, which can induce positive systematic biases, particularly for scans with low counts per ray. We propose new likelihood approximations that allow negative sinogram values without requiring zero-thresholding. With negative sinogram values, the log-likelihood functions can be nonconcave, complicating maximization; nevertheless, we develop monotonic algorithms for the new models by modifying the separable paraboloidal surrogates and the maximum-likelihood expectation-maximization (ML-EM) methods. These algorithms ascend to local maximizers of the objective function. Analysis and simulation results show that the new shifted Poisson (SP) model is nearly free of systematic bias yet keeps low variance. Despite its simpler implementation, the new SP performs comparably to the saddle-point model which has shown the best performance (as to systematic bias and variance) in randoms-precorrected PET emission reconstruction.

Penalized-likelihood image reconstruction for digital holography

Conventional numerical reconstruction for digital holography using a filter applied in the spatial-frequency domain to extract the primary image may yield suboptimal image quality because of the loss in high-frequency components and interference from other undesirable terms of a hologram. We propose a new numerical reconstruction approach using a statistical technique. This approach reconstructs the complex field of the object from the real-valued hologram intensity data. Because holographic image reconstruction is an ill-posed problem, our statistical technique is based on penalized-likelihood estimation. We develop a Poisson statistical model for this problem and derive an optimization transfer algorithm that monotonically decreases the cost function at each iteration. Simulation results show that our statistical technique has the potential to improve image quality in digital holography relative to conventional reconstruction techniques.

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.

Segmentation-free statistical image reconstruction for polyenergetic x-ray computed tomography with experimental validation

This paper describes a statistical image reconstruction method for x-ray CT that is based on a physical model that accounts for the polyenergetic x-ray source spectrum and the measurement nonlinearities caused by energy-dependent attenuation. Unlike our earlier work, the proposed algorithm does not require pre-segmentation of the object into the various tissue classes (e.g., bone and soft tissue) and allows mixed pixels. The attenuation coefficient of each voxel is modelled as the product of its unknown density and a weighted sum of energy-dependent mass attenuation coefficients. We formulate a penalized-likelihood function for this polyenergetic model and develop an iterative algorithm for estimating the unknown density of each voxel. Applying this method to simulated x-ray CT measurements of objects containing both bone and soft tissue yields images with significantly reduced beam hardening artefacts relative to conventional beam hardening correction methods. We also apply the method to real data acquired from a phantom containing various concentrations of potassium phosphate solution. The algorithm reconstructs an image with accurate density values for the different concentrations, demonstrating its potential for quantitative CT applications.

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.

An iterative transmission algorithm incorporating cross-talk correction for SPECT

Simultaneous emission/transmission acquisitions in cardiac SPECT with a Tc99m/Gd153 source combination offer the capability for nonuniform attenuation correction. However, cross-talk of Tc99m photons downscattered into the Gd153 energy window contaminates the reconstructed transmission map used for attenuation correction. The estimated cross-talk contribution can be subtracted prior to transmission reconstruction or incorporated in the reconstruction algorithm itself. In this work, we propose an iterative transmission algorithm (MLTG-S) based on the maximum-likelihood gradient algorithm (MLTG) that explicitly accounts for this cross-talk estimate. Clinical images were acquired on a three-headed SPECT camera, acquiring Tc99m emission and Gd153 transmission images simultaneously. Subtracting the cross-talk estimate prior to transmission reconstruction can result in negative and zero values if the estimate is larger than or equal to the count in the transmission projection bin, especially with increased attenuator size or amount of cross-talk. This results in inaccurate attenuation coefficients for MLTG reconstructions with cross-talk subtraction. MLTG-S reconstructions on the other hand, yield better estimates of attenuation maps, by avoiding the subtraction of the cross-talk estimate. Comparison of emission slices corrected for nonuniform attenuation reveals that inaccuracies in the reconstructed attenuation map caused by cross-talk can artificially enhance the extra-cardiac activity, confounding the ability to visualize the left-ventricular walls.

Modeling the axial extension of a transmission line source within iterative reconstruction via multiple transmission sources

Reconstruction algorithms for transmission tomography have generally assumed that the photons reaching a particular detector bin at a particular angle originate from a single point source. In this paper, we highlight several cases of extended transmission sources, in which it may be useful to approach the estimation of attenuation coefficients as a problem involving multiple transmission point sources. Examined in detail is the case of a fixed transmission line source with a fan-beam collimator. This geometry can result in attenuation images that have significant axial blur. Herein it is also shown, empirically, that extended transmission sources can result in biased estimates of the average attenuation, and an explanation is proposed. The finite axial resolution of the transmission line source configuration is modeled within iterative reconstruction using an expectation-maximization algorithm that was previously derived for estimating attenuation coefficients from single photon emission computed tomography (SPECT) emission data. The same algorithm is applicable to both problems because both can be thought of as involving multiple transmission sources. It is shown that modeling axial blur within reconstruction removes the bias in the average estimated attenuation and substantially improves the axial resolution of attenuation images.

Statistical image reconstruction for polyenergetic X-ray computed tomography

This paper describes a statistical image reconstruction method for X-ray computed tomography (CT) that is based on a physical model that accounts for the polyenergetic X-ray source spectrum and the measurement nonlinearities caused by energy-dependent attenuation. We assume that the object consists of a given number of nonoverlapping materials, such as soft tissue and bone. The attenuation coefficient of each voxel is the product of its unknown density and a known energy-dependent mass attenuation coefficient. We formulate a penalized-likelihood function for this polyenergetic model and develop an ordered-subsets iterative algorithm for estimating the unknown densities in each voxel. The algorithm monotonically decreases the cost function at each iteration when one subset is used. Applying this method to simulated X-ray CT measurements of objects containing both bone and soft tissue yields images with significantly reduced beam hardening artifacts.

Edge-preserving tomographic reconstruction with nonlocal regularization

Tomographic image reconstruction using statistical methods can provide more accurate system modeling, statistical models, and physical constraints than the conventional filtered backprojection (FBP) method. Because of the ill posedness of the reconstruction problem, a roughness penalty is often imposed on the solution to control noise. To avoid smoothing of edges, which are important image attributes, various edge-preserving regularization methods have been proposed. Most of these schemes rely on information from local neighborhoods to determine the presence of edges. In this paper, we propose a cost function that incorporates nonlocal boundary information into the regularization method. We use an alternating minimization algorithm with deterministic annealing to minimize the proposed cost function, jointly estimating region boundaries and object pixel values. We apply variational techniques implemented using level-sets methods to update the boundary estimates; then, using the most recent boundary estimate, we minimize a space-variant quadratic cost function to update the image estimate. For the positron emission tomography transmission reconstruction application, we compare the bias-variance tradeoff of this method with that of a "conventional" penalized-likelihood algorithm with local Huber roughness penalty.

Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms

Iterative coordinate ascent algorithms have been shown to be useful for image recovery, but are poorly suited to parallel computing due to their sequential nature. This paper presents a new fast converging parallelizable algorithm for image recovery that can be applied to a very broad class of objective functions. This method is based on paraboloidal surrogate functions and a concavity technique. The paraboloidal surrogates simplify the optimization problem. The idea of the concavity technique is to partition pixels into subsets that can be updated in parallel to reduce the computation time. For fast convergence, pixels within each subset are updated sequentially using a coordinate ascent algorithm. The proposed algorithm is guaranteed to monotonically increase the objective function and intrinsically accommodates nonnegativity constraints. A global convergence proof is summarized. Simulation results show that the proposed algorithm requires less elapsed time for convergence than iterative coordinate ascent algorithms. With four parallel processors, the proposed algorithm yields a speedup factor of 3.77 relative to single processor coordinate ascent algorithms for a three-dimensional (3-D) confocal image restoration problem.

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.