Comparison of Lesion Detection and Quantification in MAP Reconstruction with Gaussian and Non-Gaussian Priors

Statistical image reconstruction methods based on maximum a posteriori (MAP) principle have been developed for emission tomography. The prior distribution of the unknown image plays an important role in MAP reconstruction. The most commonly used prior are Gaussian priors, whose logarithm has a quadratic form. Gaussian priors are relatively easy to analyze. It has been shown that the effect of a Gaussian prior can be approximated by linear filtering a maximum likelihood (ML) reconstruction. As a result, sharp edges in reconstructed images are not preserved. To preserve sharp transitions, non-Gaussian priors have been proposed. However, their effect on clinical tasks is less obvious. In this paper, we compare MAP reconstruction with Gaussian and non-Gaussian priors for lesion detection and region of interest quantification using computer simulation. We evaluate three representative priors: Gaussian prior, Huber prior, and Geman-McClure prior. We simulate imaging a prostate tumor using positron emission tomography (PET). The detectability of a known tumor in either a fixed background or a random background is measured using a channelized Hotelling observer. The bias-variance tradeoff curves are calculated for quantification of the total tumor activity. The results show that for the detection and quantification tasks, the Gaussian prior is as effective as non-Gaussian priors.


INTRODUCTION
Statistical image reconstruction methods have been developed for emission tomography to improve the signal-tonoise ratio (SNR), for example, [1][2][3]. The most popular maximum-likelihood (ML) reconstruction algorithm is the expectation-maximization (EM) algorithm [1,4]. However, the ML estimate can be very noisy because emission tomography is an ill-posed problem. Hence some form of regularization is needed to obtain a useful image. Bayesian methods regularize the solution by using a prior probability distribution on the image. The prior probability distribution plays an important role in Bayesian image reconstruction. The most commonly used prior is the Gaussian prior which strongly discourages sharp transitions in images. To preserve edges, non-Gaussian priors have been proposed [3,5,6]. Empirical results showed that images reconstructed with edge-preserving non-Gaussian priors have less mean-squared error than those reconstructed with the Gaussian prior. However, the effect of edge-preserving priors on clinical tasks is not obvious.
Gifford et al. [7] compared a quadratic penalty function with Huber penalty functions in a penalized-EM algo-rithm for tumor detection. However, the comparison was done as a function of iteration, so the result is algorithmdependent. Nuyts and Michel [8] compared maximum a posteriori (MAP) reconstruction with a relative difference prior to post-smoothed ML reconstruction for hot lesion detection and found similar performance between MAP and post-smoothed ML reconstructions. In this paper we evaluate the performance of Bayesian image reconstructions with the Gaussian and non-Gaussian priors for hot lesion detection and region of interest quantification, which are two major clinical applications of emission tomography. Some preliminary results were reported in [9]. This paper is organized as follows. We briefly review Bayesian image reconstruction and different prior functions in Section 2. The methods for evaluating image quality are described in Section 3. The results are presented in Section 4. Finally, the conclusion is drawn in Section 5.

BAYESIAN IMAGE RECONSTRUCTION
In emission tomography data are well modeled as a collection of independent Poisson random variables with the meanȳ ∈ R M×1 related to the unknown image, x ∈ R N×1 , through an 2 International Journal of Biomedical Imaging where P ∈ R M×N is the detection probability matrix with element (i, j) equal to the probability of an event produced in voxel j being detected by sinogram bin i, and r ∈ R M×1 accounts for the presence of scatters and randoms in the data. The appropriate log-likelihood function is given by where y is the measured data. The ML reconstruction can be obtained by maximizing (2). However, the ML estimate can be very noisy because emission tomography is an ill-posed problem. Bayesian methods regularize the noise by using a prior probability distribution on the image. Most image priors use a Markov random field with a Gibbs distribution of the form where U(x) is the energy function, β is the hyperparameter that controls the resolution of the reconstructed image, and Z is a normalization constant. The Markovian properties of these distributions make them theoretically attractive as formalism for describing empirical local image properties, as well as computationally appealing. The energy function U(x) often contains potentials defined on pair-wise cliques of neighboring voxels: where N j denotes the set of neighboring voxels of voxel j, κ jk are weighting factors, and V (·) is the potential function. A wide range of potential functions has been studied in the literature that attempt to produce local smoothing while not removing or blurring true boundaries or edges in the image. All have the basic property that they are monotonic nondecreasing functions of the absolute intensity difference The potential function of a Gaussian prior is a quadratic function It produces smooth images with very low probability of sharp transitions in intensity. In an attempt to increase the probability of sharp transitions, Bouman and Sauer [3] propose using the generalized p-Gaussian model where V (x) = |x| p , 1 < p < 2. An alternative function is the Huber prior in which V (·) is defined as [10] V When δ is small, the Huber function approximates the absolute value function. It is probably the most edge-preserving prior with a convex potential function. Other potential functions with similar behavior are V (x) = log cosh(δx) [6] and [11]. In an attempt to produce even sharper intensity transitions, nonconvex functions have also been proposed. One example that we will study is the Geman-McClure prior [5], of which Figure 1 shows the potential functions of the Gaussian, Huber, and Geman-McClure priors. Note that in practice both the Huber prior and Geman-McClure prior can approach performance of the Gaussian prior by setting δ to be sufficiently large. To demonstrate the difference between these potential functions, we show two images in Figure 2. For the Huber prior (δ 1) the two images are equally probable, while the Gaussian prior strongly favors the cone image (with an energy ratio of 50 : 1) and the Geman-McClure prior favors the disk image. Combine the likelihood function and prior distribution, the MAP reconstruction x is found by maximizing the logposterior density function: For priors with a convex potential function, (8) generally has a unique solution. When nonconvex potential functions are used, there may exist multiple local optima and the solution of most deterministic optimization algorithms will depend on the initial image.

Computer simulation
We conduct computer simulations to study the effect of non-Gaussian priors on lesion detection. We simulate imaging of a prostate tumor using C-11 choline [12]. The simulated PET system has similar parameters as an ECAT HR+ clinical scanner (CPS, Knoxville, TN). It has 576 detectors forming a ring with radius 41.3 cm. The phantom has a body shape that is obtained from a patient image. The background has nearly uniform uptake of the radiotracer and has an attenuation coefficient of 0.0095 mm −1 . We place a round hot spot of different diameters (5 mm and 15 mm) at the center of the image to simulate a prostate tumor. Different pixel sizes are used in data generation and image reconstruction to introduce some model mismatch. For data generation, phantoms are represented by 256×256 2-mm square pixels, whereas 128×128 4.5-mm square pixels are used in reconstruction. Figure 3 shows various phantom images that we used. Figure 3(a) is a fixed uniform background. Figures 3 are three random backgrounds obtained by superimposing the background image in Figure 3(a) with a realization of lumpy backgrounds [13]. The lumpy backgrounds are modeled as where G(b, σ 2 , r k ) is a Gaussian blob with variance σ 2 and height b centering at a random location r k , and K is a Poisson random variable. The mean of K is set to 100. Two sets of b and σ are used: b = 0.02 and σ = 32 mm in Figure 3 we also add a hot region with activity to background ratio of 4 : 1 to mimic possible bladder uptake. In all four cases the mean activity of the background is about 0.20. For each type of the background we generate three groups of data: background only, background with the 5mm lesion, and background with the 15-mm lesion. Each group consists of 1000 independent identically distributed data sets. The expected total number of detected events is about 200 000. All data sets are independently reconstructed using a preconditioned conjugate gradient algorithm with the Gaussian prior, Huber prior, and Geman-McClure prior and different β and δ values. Five hundred iterations are used to ensure effective convergence of the algorithm. All reconstructions start from a uniform image. For Geman-McClure prior, the reconstructed image may correspond to a local optimum of the log-posterior density function because the objective function is nonconvex.

Lesion detection
Detection of cancerous lesions is one major task of emission tomography. A standard methodology for studying lesion detectability is the receiver operating characteristic (ROC) study that compares true positive versus false positive rates for human observers for the task of lesion detection. Numerical observers based on signal-detection theory have been developed to mimic human performance [14]. For a given reconstructed image x, a linear numerical observer computes a test statistic (a scalar-valued decision variable), η( x), by where t is the observer template. The detection performance can be measured by the SNR of η( x), which is defined as where H 0 is the null hypothesis representing lesion absent, H 1 is the hypothesis representing lesion present, Σ x|H1 and Σ x|H0 are the conditional covariance matrices of x under hypotheses of H 1 and H 0 , respectively, and is the difference between the mean reconstructions under the two hypotheses. When η( x) is normally distributed, the area under the ROC curve (AUC) is related to the SNR by where erf(x) is the error function. We use a channelized Hotelling observer (CHO), which has been shown to correlate with human performance [15][16][17][18][19]. The test statistic of CHO is where U denotes frequency-selective channels that mimic the human visual system, n is the internal channel noise that models the uncertainty in the human detection process [20] with zero mean and covariance K N , and K is the covariance of the channel outputs, that is, In this work the channel functions are the differences of four Gaussian functions with standard deviations σ = 2.653, 1.592, 0.995, and 0.573, respectively [21]. The internal noise is modeled as uncorrelated noise with where σ 2 i is the data variance in the ith channel output [22]. Monte Carlo reconstructions are used to calculate the expectation of the reconstruction and covariance matrices. The SNR of CHO is then calculated by because the SNR calculated from (15) is meaningful only when η( x) is normally distributed, we also calculate AUCs from empirical ROC curves using numerical integration.

Quantification performance
Another clinical task in emission tomography is to quantify the uptake of radioactive tracer in a region of interest (ROI). This can be written as where t is the indicator function of the ROI, that is, t j = 1 if voxel j is inside the ROI, and t j = 0 otherwise. The accuracy of the quantification can be measured by the bias and variance of η Q ( x) as where x denotes the true tracer uptake. Note that the ROI quantification is only performed on images in which a lesion is known to be present. We use the bias versus variance tradeoff curve to evaluate the quantification performance. The above equations are calculated from Monte Carlo reconstructions.

Figures 4 and 5 show examples of reconstructed images with different priors.
For all priors the reconstructed images become less noisy as we increase the hyperparameter β. The images reconstructed with the Gaussian prior have blurred edges, whereas the images reconstructed with the Huber prior and Geman-McClure prior tend to form piece-wise constant regions. Figure 6 shows the variation of AUC as a function of β and δ for detecting the 15-mm lesion in the lumpy background shown in Figure 3(c). In this case, the optimum parameters for the Gaussian, Huber, and Geman-McClure priors are β = 30, (β = 10, δ = 0.1), and (β = 0.3, δ = 0.1), respectively. Comparing to the results of the fixed background (not shown), we found that lumpy backgrounds reduce lesion detectability at low resolution (large β) and slightly decrease the optimum β value for lesion detection.
To give a fair comparison between different priors, we choose to compare the maximum SNR and maximum AUC of each prior. The results are shown in Figure 7. The error bars are computed using a bootstrap method. In all detection studies, the contrasts of the 5-mm and 15-mm lesions are 3 and 0.9, respectively, which are selected to obtain a reasonable detectability (AUC ≈ 0.9). In each case, SNR and AUC give similar ranking for the optimum performances of the three priors. No statistically significant advantage is found for non-Gaussian priors. For lesions with lower contrast, we expect the performances of the three priors to be even closer. We notice that some differences in AUC do not correspond well to the differences in SNR (see Figures 7(c) and 7(d)), which indicates that the test statistics of the numerical observer does not follow a Gaussian distribution.   Figure 8 shows the bias versus variance tradeoff curves for quantification. The contrast of the lesion is 3.0 for the 5-mm lesion and 0.9 for the 15-mm lesion, respectively. The ROIs were obtained from the original phantom image. Bias and standard deviation are normalized to the total activity inside each ROI. Here we plot all the cases. It is interesting to see that Gaussian prior seems to set a lower bound for the Huber and Geman-McClure priors (except in Figure 8(h)). Considering that both Huber and Geman-McClure priors include Gaussian prior as a special case, no improvement in ROI quantification is found by using the Huber or Geman-McClure prior. In many cases, the performance of the non-Gaussian priors is much worse than that of the Gaussian prior, indicating that hyperparameter selection is more important for non-Gaussian priors. In Figures 8(b), 8(d), 8(f) there is a kink in the Gaussian biasvariance curve at low noise levels. This is because at such low resolution the reconstructed background is no longer uniform, but forms a dome, which artificially increases the ac-tivity inside the ROI and, hence, reduces the bias. The decrease in bias at low noise levels in Figures 8(g) and 8(h) is mostly caused by the spill-over effect of the nearby hot region ("bladder").
We also studied lesions with higher contrast (9 for the 5-mm lesion and 2.7 for the 15-mm lesion, respectively) and found very similar results [9]. To investigate the bias-variance tradeoff for large regions, we quantify the total activity in the hot "bladder" in Figure 3(d). The results are shown in Figure 9. Even for this large region with 4 : 1 activity ratio, we do not see any advantage of using the Huber and Geman-McClure priors.

CONCLUSION AND DISCUSSION
We have compared the performance of three representative priors for lesion detection and ROI quantification. The Gaussian prior is the most commonly used prior in emission reconstruction; the Huber prior is probably the most  edge-preserving prior among all priors with a convex potential function; and the Geman-McClure is a typical edgepreserving prior with a nonconvex potential function. Note that both the Huber and Geman-McClure priors can approach the Gaussian prior by setting δ to be sufficiently large. Thus we focus on whether the Huber and Geman-McClure priors can outperform the Gaussian prior. In all the cases that we have tested, we have not observed any significant improvement. The results show that for the detection and quantification tasks that are considered here, the Gaussian prior is as effective as the more complex non-Gaussian priors.
We should note that while we have investigated each prior with a range of β and δ values, it is still possible that the results may not reflect the best performance for non-Gaussian priors because of the lack of theoretical guidance on the hyperparameter selection for non-Gaussian priors. Nonetheless, the results indicate that hyperparameter selection is ex-tremely important for non-Gaussian priors. For the Geman-McClure prior, the simulation results presented here may not correspond to the global maximum of the log posterior distribution because the objective function is nonconcave. It is possible that the performance of Geman-McClure prior may be improved by using deterministic or stochastic annealing techniques at an expense of increased computational cost [23].
In this paper we used a channelized Hotelling observer for lesion detection and a simple ROI estimator for quantification because of their popularity. It is possible that some results might change if different observers or estimators were used. In ROI quantification we defined ROI using the original phantom data. If an ROI is to be delineated on the reconstructed image, the error in ROI delineation will also affect the ROI quantification [24]. Such effect is not included in this study. We expect that with the recent development of combined PET/CT and SPECT/CT scanners, high-resolution   anatomical images will help to reduce the error in ROI definition.

ACKNOWLEDGMENT
This work is supported by National Institutes of Health under Grant R01 EB000194.