Analysis of Resolution and Noise Properties of ... - IEEE Xplore
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
413
Analysis of Resolution and Noise Properties of Nonquadratically Regularized Image Reconstruction Methods for PET Sangtae Ahn, Member, IEEE, and Richard M. Leahy*, Fellow, IEEE
Abstract—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. Index Terms—Image resolution, noise, positron emission tomography (PET), regularization, statistical image reconstruction.
I. INTRODUCTION
S
TATISTICAL image reconstruction methods such as penalized-likelihood (PL) or maximum a posteriori (MAP) estimation are widely used for positron emission tomography (PET). The regularization or penalty function regularizes the ill-conditioned tomographic inverse problems. In other words, it controls a trade-off between resolution and noise in reconstructed images. The advantage of using penalty functions is that one can often predict accurately and/or efficiently the resolution and noise properties of reconstructed images and can design penalty functions for desired properties [1]–[8]. In contrast, other methods of stabilizing noise, such as early termination of
Manuscript received July 24, 2007; revised October 15, 2007. This work was supported by the National Institute of Biomedical Imaging and Bioengineering under Grant R01 EB00363. Asterisk indicates corresponding author. S. Ahn is with the Signal and Image Processing Institute, University of Southern California, Los Angeles, CA 90089 USA. *R. M. Leahy is with the Signal and Image Processing Institute, University of Southern California, Los Angeles, CA 90089 USA (e-mail: [email protected]. edu). Digital Object Identifier 10.1109/TMI.2007.911549
unregularized expectation maximization (EM) or ordered subsets EM (OSEM) [9], have only limited control over resolution and noise properties through choice of the number of iterations and subsets and post-reconstruction smoothing. Quadratically regularized estimators are approximately linear for reasonably high signal-to-noise ratio (SNR). Therefore, the resolution and variance of reconstructions can be estimated analytically without expensive Monte Carlo (MC) simulation [1], [2]. Furthermore, the use of Fourier transforms that take advantage of (approximate) Toeplitz structures in the system resolution matrices [10] enables one to efficiently predict resolution and variance as well as derived figures of merit for specific tasks such as detection and quantification and also to optimize regularization parameters for these tasks [4], [11]–[19]. However, the quadratic penalty tends to over-smooth sharp edges in reconstructed images. Many alternative nonquadratic penalty functions with edge-preserving properties have been proposed [20]–[26]. However, there has been little research on quantitative analysis of the resolution and noise properties for nonquadratic regularization in contrast to the rich literature on analysis and design of quadratic regularization. The lack of analytical results for nonquadratic regularization limits our ability to control the properties of reconstructed images through systematic selection of regularization parameters. Lalush and Tsui [27], [28] analyzed the smoothness in each iterate of the one-step-late (OSL) MAP-EM algorithm [22] in terms of the previous iterate by using the first derivative of a nonquadratic penalty function; however, the analysis does not apply to a case where iteration-independent and algorithm-independent MAP images are of interest. In this paper, we derive new expressions for estimating resolution and variance accurately and efficiently for nonquadratic regularization. In Section II, we define a local perturbation response (LPR), a generalization of the local impulse response (LIR) [1], to quantify resolution. We derive a new approximate expression in Section III-A for a linearized local perturbation response (LLPR), which approximates the LPR well for high SNR. While the LLPR is only implicitly expressed, it does accurately capture the nonlinear functional dependence of the perturbation response, which is intrinsic to nonquadratic regularization. In contrast, the conventional linearized local impulse response (LLIR) [1] is useful only when the superposition principle or linearity holds, which is not the case for nonquadratic regularization even for very high SNR. Although the expression for LLPR does not have a closed-form solution, we can solve the equation efficiently using a local estimation strategy as described in Section III-A. The expression allows us to achieve
0278-0062/$25.00 © 2007 IEEE
414
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
spatially uniform perturbation responses, akin to uniform spatial resolution in quadratic regularization, by using the certaintybased modified penalty [1], [5], as discussed in Section III-E. We also derive new expressions for the variance in the context of nonquadratic regularization in Section III-B. The implicit expression gives accurate estimates of the LPR and variance, even for low SNR, through efficient MC simulation. In Section III-C we describe a plug-in method that can be used when only a single noisy dataset is available and, in Section III-D, the use of a table lookup method to estimate resolution and variance. We perform simulation studies using a 2-D shift-invariant PET system to validate the methods in Section IV, and present concluding remarks in Section V. II. BACKGROUND
of horizontally or vertically neighboring pixels, for diagonally neighboring pixels, and 0 otherwise. Many types of nonconvex or nondifferentiable potential functions have been proposed [20], [21], [24], [32]. However, here we consider only convex and differentiable potential functions (e.g., [22], [23], [26]) because 1) our analysis requires the regularization function to be differentiable, 2) nonconvex potential functions may lead to unstable estimates in the sense that the estimator is not a continuous function of the data [24], and 3) it is more difficult to optimize nonconvex or nondifferentiable objective functions. Additionally, we assume the first derivative is absolutely continuous and therefore the second derivative exists almost everywhere [33, p. 109]. pairs of neighboring pixels Suppose there are a total of such that for . Define a differencing matrix as
A. Statistical Image Reconstruction is defined as A PL or MAP estimate a maximizer of an objective function as follows: (1) where denotes the PET emission data, is the number of pixels, is the number of detector pairs, and denotes matrix or vector transpose. The objective function is given by
and define for . The regularization function in (5) can then be rewritten as
The Hessian of
can be conveniently written as
(2) (6) where is the log-likelihood, is the regularization function or the energy function for a Gibbs prior, and is a global regularization parameter also known as a hyperparameter. Assuming the data follow the Poisson distribution [29], the log-likelihood can be written, neglecting constants, as
is the mean or noiseless data. The mean is affine in the object
denotes a diagwhich exists almost everywhere, where onal matrix that is formed by taking the arguments as diagonal elements. The simplest and most widely-used penalty is , which penalizes the squared difference of neighboring pixel values and consequently encourages smoothness in reconstructed images. However, the quadratic penalty tends to over-smooth sharp edges. We consider the following Huber penalty [34], [35] as a representative edge-preserving nonquadratic penalty, although the methods in this paper also apply to other penalties:
(4)
(7)
(3) where measurement
denotes the system or projection matrix charwhere acterizing the PET scanner and is the mean contribution of background events such as randoms and scatters. In this paper, we assume has already been estimated (for example, see [30] and [31]). For the regularization function, we focus on pairwise roughness penalty functions of the following form:
for some . The Huber penalty function quadratically penalizes small differences in neighboring pixels while applying a linear penalty to larger differences, with the transition determined by the parameter . Fig. 1 illustrates the quadratic and the Huber potential functions and their first and second derivatives. B. LPR
(5) where are nonnegative weights and is a potential function. and are symmetric, that is, We assume that and . Usually, in 2-D, is 1 for a pair
We define a LPR by which we quantify resolution properties. consists of a background Suppose that a true object and a signal of interest (8)
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
415
the resolution properties of nonquadratic regularization through the LPR as a function of perturbation . III. THEORY
Fig. 1. Quadratic and Huber potential functions and their first and second derivatives.
In Sections III-A and III-B, we develop methods which and require that the background and perturbation images, , are known. Such methods are useful for investigating the effects of the choice of penalty and its parameters on image quality metrics such as resolution, noise, lesion detectability, and quantification accuracy [4], [11]–[19] for various kinds of perturbations such as lesions of different sizes and contrasts. In Section III-C, we describe plug-in methods, which can be used when only a single noisy dataset is available. We then describe simplified methods for shift-invariant imaging systems in Section III-D and present a method for achieving spatially uniform perturbation responses in Section III-E. The analytical results in Sections III-A and III-E were partly presented in [40]. A. Asymptotic Resolution Analysis
For example, may represent a lesion or an organ of interest. We are interested in the effect of regularization upon recovery of the signal . In other words, we aim to study what would look like in a reconstructed image for a specific choice of regas the mean perturbation in the ularizer. We define a LPR reconstructed images caused by (9) where and are reconstructed images when the signal present and when absent, respectively. That is
(13) where is given in (1), and are defined in (12), and and in (10). By substituting (13) into (9), we define the linearized local perturbation response (LLPR) as
is (14) (10)
where when
1) LLPR: The LPR defined in (9) is intractable to analyze due to the mean operators. Fortunately, the mean of an estimator can be approximated by a noiseless estimate for high SNR1 [2]
is given in (1), and and represent the noisy data is present and absent, respectively. In other words (11)
where (12) ” In the equations above, is given in (4), and “ means that the elements of are independent Poisson random . Similar types of regional perturbation variables and responses have been used for analyzing and compensating for partial volume effects (e.g., see [36] and [37]). The LPR is a generalization of point, line, and edge spread is of a particular type. If is an functions where the signal impulse with an infinitesimal amplitude, then the LPR reduces to the LIR [1]. LIR is useful for studying resolution properties in quadratic regularization, which is approximately linear [5], [38] and thus the superposition principle holds. The LIR has been used as a key tool for quadratic penalty design methods [1], [4], [5], [39], [38]. However, in nonquadratic regularization, resoluand LIR tion properties depend heavily on the perturbation is of limited utility. For nonquadratic regularization, LIR based prediction of resolution and variance is inaccurate, particularly when SNR is low, as shown in Section IV. Therefore, we study
akin to the LLIR [1]. It is called a linearized response because the approximation is based on the first-order Taylor series truncation of the estimator [1]. That is, if , then . In Section IV-B, we show how well LLPR approximates LPR for various SNR levels. is to A straightforward method for computing the LLPR and using noiseless data and and then reconstruct to use (14); however, this gives little insight in the sense that it and . We is difficult to examine the relationship between instead derive an approximate expression that relates to the true perturbation . The new expression suggests a computationally efficient strategy and provides useful insights into resolution properties. Neglecting the nonnegativity constraint, which is reasonable for noiseless reconstructions [1], the noiseless estimates, and , satisfy the first-order optimality conditions (15) (16) 1Under mild conditions including convexity of and P y 6= 0, the objective function 8 in (1) is strictly concave [41, Lemma 1]. It follows that x ^( 1 ) in (1) is continuous [24, Theorem 1]. Note that if y ’s are iid random vectors, then y =K converges in probability to E [y ] by the weak law ^, x^( y =K ) of large numbers [42, p. 112]. Due to the continuity of x ^(E [y ]) [42, p. 124]. Therefore, as the photon converges in probability to x count increases (K ! 1), the estimate converges in probability to the noiseless estimate.
416
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
where is the column gradient operator with respect to the first argument. Consider the first-order Taylor series expansion about of
where . Note the second derivative term disappears in the above equation. Combining (4) and (15)–(19) yields
(17)
(22) (23)
is the Hessian operator with respect to the first where th element of the operator is argument, and the . The above approximation is reasonably accurate is linear in and varies slowly with because through its projection [see (3)]. However, in (16) may not be well approximated by a Taylor series truncation is nonlinear. For example, if is nonquadratic and thus whereas Fig. 1 shows the nonlinear behavior of is perfectly linear. Therefore, we take the zeroth-order Taylor expansion with the remainder term or equivalently the fundamental theorem of calculus [43, p. 292] (18)
Although (22) may apply to other system models and noise models, here we focus on the linear emission tomography system with the Poisson noise model. Using (3) and (4), one can compute the following derivative terms: (24) (25) where mating
is the Fisher information matrix for estigiven by (26)
Since
The approximation in (24) is reasonably accurate because 1) is a blurred version of, and hence very similar to, and 2) , and are all smoothing operators [1]. Substituting (4), (24), and (25) into (23) leads to (27) by (6), one can rewrite (18) as (19) where
This is a new expression that relates the LLPR to the true perturbation . To simplify the above expression, we consider the following factored form of the system matrix:
is defined as
(28) (20)
with (21) Interestingly, the above form of has appeared in different to design a contexts. For example, Huber has used quadratic surrogate for the penalty function with the aim of constructing monotonic optimization algorithms [44, p. 184]. Charbonnier et al. have described certain conditions for edge-pre[25]. Also, note that (20) serving properties in terms of and (21) are similar to (36) in [4], which is an ad hoc expression to predict the contrast recovery when the Huber penalty is used. If the penalty function is quadratic, then is a constant maof the penalty. By (20) and (21), trix, which is the Hessian the second term on the right hand of (19) can be calculated as
where denote ray-dependent factors including detector efficiencies, attenuation, and dead time, and is the geometric projection matrix possibly including space-variant sinogram blurring [1], [4]. Then the Fisher information matrix in (26) can be approximated as
(29) where (30) (31) because is concentrated about its diagonal and , which is interpreted as a measure of certainty, varies slowly with [1]. The approximation in (29) separates into an object-dependent and a precomputable object-independent part . part
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
Suppose the perturbation is localized around pixel . For a will be similarly localreasonable degree of regularization, ized. This localized perturbation assumption allows us to further simplify the approximation for
(32) is concentrated about its diagonal and because slowly [1]. Similarly
varies
(33) Substituting (32) and (33) into (27) yields a simplification of (27) (34) (cf. [1, (34)]). The approximate expression (27) for the LLPR and its simplified version (34) are the key results of this subsection, which can be viewed as a generalization of the approximations for the LLIR in [1]. We now briefly compare the LLPR and LLIR approximations. If we used the first-order Taylor series truncation instead of the Taylor expansion with remainder for in (19), we would obtain the following explicit expression for : (35) where (36) The th column of is the LLIR for pixel [1]. If in (27) is replaced with , then (27) becomes (35) and (36). Thus one can view (35) and (36) as being obtained in (27). The approxiby discarding the dependency of on mation (35) provides an explicit expression for , which is reasonably accurate for quadratic regularization. However, it does not predict the LLPR accurately for nonquadratic regularization particularly in a low SNR case whereas (27) does, as shown in Section IV-B. 2) Efficient Computation Using Local Estimation: One can in (27) by solving the following nonlinear least compute squares problem:
417
However, solving (37) is no less costly than performing noiseby (14). less image reconstructions and computing is localized about pixel , then If the perturbation and will also be concentrated about pixel . Based on this observation, one can downsize the problem (37) using local estimation. First, we define a set of region of support about pixel for some pixels within an such that . Next, we restrict the parameter space to by taking only the pixels in the set as optimization variables for the least squares problem. Equivalently, we use as a weight matrix in (37) where is an indicator function on . Choosing an appropriate depends and the regularization parameters such as on the support of . As increases, the LLPR computed by the local estimation method will become more accurate but the computational cost will increase. This downsizing strategy is somewhat similar to the local backprojection technique used in efficient quadratic penalty design [38] or the small volume support technique used in efficient calculation of resolution and covariance for fully 3-D SPECT [6]. An advantage of computing the LLPR by (37) rather than using (14) through noiseless image reconstructions is that one can use (37) under certain conditions even when only a single noisy dataset is available, as discussed in Section III-C. B. Efficient Monte Carlo Simulation Method for Variance Analysis To estimate the variance of noisy reconstruction , defined in (10), one can use a MC simulation method where a number of noisy datasets are generated and the sample variance is computed from noisy reconstructions for each pixel. However, we present a more efficient MC simulation method here. Neglecting again the nonnegativity constraint, the following first-order optimality conditions are satisfied: (39) (40) where and are defined in in (11), (12), and (13), respectively. As in Section III-A, combining the first-order Taylor around and the zeroth-order series truncation of around with the remainder term Taylor expansion of yields
(41) where the deviation
of
from
is defined as
(37) (42) where
is a diagonal weight matrix and Computing the derivatives as in (24) and (25) and substituting them into (41) leads to
Similarly, one can solve a least squares problem based on (34) by using (38)
(43)
418
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
where
Note that
is approximately white noise, that is, because and
and
. For the factored system model in (28), the right-hand side in (43) can be approximated as
will be . However, the nonlinin in (43) or (46), then earity in nonquadratic regularization manifests itself in nonzero , causing the LLPR to deviate from the LPR particularly for low SNR. If we used the first-order Taylor series truncation for in (40) in place of the Taylor expansion with the remainder term, and the covariwe would obtain an explicit expression for ance of , equal to the covariance of , could be written as
(48) (44) since the means and covariances of both sides above approximately match in view of (28) and (29). Combining (29), (43), and (44) yields
which is equivalent to the covariance expression derived in [2]. For nonquadratic regularization, the above expression becomes inaccurate as the photon count decreases, as shown in Section IV-B. C. Plug-In Method
(45) Again, if our region of interest is localized around pixel , (45) can be locally approximated as (46) The approximate expression (43) for the deviation and its simplified version (46) are the main results of this subsection. in the neighborhood One can estimate the variance of of pixel using (43) or (46) as follows. First, one generates a number of realizations of white noise , and for each reby taking a small alization one solves (43) or (46) for region of support about pixel and using the local estimation method, as described in Section III-A-2. The variance of is estimated by the sample variance of in view of (42). This MC simulation method with local estimation is much less expensive than the ordinary MC simulation method with full image reconstruction, as shown in Section IV-B. The efficient MC simulation method with local estimation can also provide an accurate estimate of the LPR. As the photon count decreases, the accuracy of the LLPR ’s approximation to LPR becomes worse. However, using the MC results, one can predict the LPR accurately even for low SNR. One can rewrite the LPR , defined in (9), as follows:
(47) where . Therefore, can be viewed as a correction term when one predicts the LPR from the LLPR given in (14). Consider a special case where the background image is locally uniform. In this case, since locally and therefore , one can compute by adding , estimated by the sample mean from the MC simulation using (43) or (46), to , estimated by solving (37). This locally uniform background assumption is similar to is not locally uniform, one would [18, Approximation 1]. If from another Monte Carlo simulation have to estimate is linear using instead of in (43) or (46). Note that if
In Sections III-A and III-B, we derived approximate expressions (34) and (46) for the LLPR and variance, respectively, for a factored system model. It is straightforward to use those expressions in a simulation study where the background and perturbation images, and , are known. However, in practice, only a single noisy dataset is available. We address three issues that arise when the methods given in Sections III-A and III-B are to be applied to this practical case. First, one can estimate from the noisy data and substitute in (34) and (46) [1], [2], [4], [5]. The plug-in esthem for varies slowly with the image timate is fairly accurate since . We use the generalized error lookup table and method [45] to estimate in Section IV. Second, the expression in (34) requires , which is unavailable in practice. It is not acceptable to substitute a noisy is sensitive to a noisy estimate reconstruction because of and moreover it is not possible to estimate only the background image in the presence of the signal . So we is locally uniform in the assume the background image region of interest, as is , and we restrict the analysis using the plug-in method to the locally uniform background case. in Based on this assumption, one can replace for some where (34) with and , consequently removing the dependency on the background image. Also, in (46) can be approximated by since under the locally uniform background assumption. is needed in (34); however, it Third, the true perturbation is unknown in practice. But one can investigate the resolution and variance properties for various typical perturbations as in a signal-known-exactly detection task (e.g., [12], [17]). One could also establish a statistical model of the perturbation and study the properties of reconstructed images as in signal-known-statistically tasks (e.g., [14]). D. Prediction of Resolution and Variance for Shift-Invariant Imaging Systems For a localized perturbation , one can use (34) to estimate and use (46) to estimate deviations , from the LLPR
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
which the variance of can be estimated. The expressions (34) and (46) are functions of 1) the true perturbation , 2) the regularization parameters, e.g., and for the Huber penalty, 3) measures of certainty , 4) (the mean reconstruction of the) , background image , and 5) the system response matrix which is data independent. As we add more assumptions, we can further simplify the situation. Assume that the background image , and hence , are locally uniform in the region of inwith in (34) terest. Then one can replace with in (46) as discussed and replace are no longer functions of in Section III-C, and now and the background image. Additionally, assume that the system is and will also be shift-invariant. Then (34) and (46) for shift-invariant, and the perturbation response and its statistical properties are not dependent on the location of the perturbation in the image domain. To summarize, if the background image is locally uniform and the system is shift-invariant, then for a given perturbation and regularization parameters and , one can precompute the LLPR (or LPR) and variance for a range of different ’s. Then the LLPR (or LPR) and variance in practice can be esregardless of the timated by table lookup using an estimate perturbation location and the background image. A similar table lookup strategy has been proposed in [3], although for quadratic regularization only.
419
Fig. 2. Digital phantom used in simulation studies. The background, the left hot disc and the right cold disc have relative emission activities of 2, 3, and 1, respectively. This figure also shows the circular perturbation of diameter 3 pixels in white that we used in the studies presented.
(34) in Section III-A, one can obtain the following approximawhen is used as a regularization function for the LLPR is localized about pixel : tion, assuming the perturbation
(50) where (51) with being defined as the th pair of neighboring is also localized about pixel , the following pixels. Since approximation is valid:
E. Spatially Uniform Perturbation Response
(52)
For nonquadratic regularization, the resolution is a function of the perturbation. When the LLPR , which approximates the LPR well for high SNR, is viewed as a function of is shift-variant even for a perturbation , the system shift-invariant imaging system partly due to the contribution of the background image through the ’s. In this subsection, we make the LLPR system shift-invariant using appropriate penalty weights. This is an extension of the quadratic penalty design for uniform spatial resolution [1], [5] to a nonquadratic penalty case. Note that the nonquadratic penalty is intended to yield nonuniform resolution and consequently what we aim to achieve is not uniform spatial resolution but rather a spatially uniform functional dependence of the LLPR on the perturbation. Here, we assume the imaging system is shift-invariant and the background is locally uniform. Using the implicit expression for the LLPR in (34), we show that the spatially uniform perturbation response of the LLPR can be achieved by using the following certainty-based modified penalty, which is proposed in [1]: (49) are estimated from noisy data, and and are the where same as given in (5). Following a similar approach to deriving
where
is defined in (20). Substituting (52) into (50) yields (53)
and we assume the background image is since locally uniform about pixel . Therefore, using the modified . penalty in (49) removes the dependence of LLPR on It can be seen in (53) that if one considers a system whose input is perturbation and whose output is LLPR , then the system is shift-invariant. Therefore, the modified penalty achieves the spatially uniform perturbation response properties. IV. RESULTS A. Two-Dimensional PET Simulation To validate the methods developed in Section III, 2-D PET emission scans were simulated. The synthetic emission phantom used in this study had relative emission activities of 1, 2, and 3 in the right cold disc, the warm elliptic background, and the left hot disc, respectively, as shown in Fig. 2. Since some of the methods we test are for shift-invariant imaging systems from Sections III-D and III-E, we generated a nearly shift-invariant 2-D PET system matrix using ASPIRE [46] to simplify the validation. The sinograms had 192 radial bins and 160 angles uniformly sampled over 180 . Nonuniform detector efficiencies were considered, and attenuation factors, from a uniform attenuation coefficient 0.095 cm , were included. The background
420
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
contribution factors corresponded to a uniform field of 10% of the total count. Reconstructed images were obtained, unless stated otherwise, using 20 iterations of fast ordered subset separable paraboloidal surrogates (OS-SPS) with 16 subsets [41] followed by 20 iterations of convergent SPS with optimal curvature [47] to ensure convergence, with filtered backprojection (FBP) initialization. The image had 128 128 pixels and a first-order neighborhood system was used. For a perturbation , a uniform circle of diameter 3 pixels was placed in the center of either the background, hot region or cold region, as shown in Fig. 2. For the local estimation described in Sections III, to solve the nonlinear least squares problems in (37), we used the function in the Optimization Toolbox of Matlab (Mathworks, Natick, MA). B. Accuracy of Prediction We evaluate the accuracy of the approximate expressions given in Sections III-A and III-B. The uniform circular activity located in the center of the background is taken as a perturbawith magnitude . Noisy datasets following tion the Poisson distribution were generated for various photon counts (0.25, 0.5, 1, 2, 4, and 8 millions) by scaling in (28) appropriately, in the presence and absence of the perturbation. was estimated by substituting the sample means, The LPR computed from 500 realizations, into the means in (9). The Huber penalty was used for reconstruction and the parameters and where is the total photon were set to count in millions in order to fix the resolution for different counts. The parameter was chosen to be 0.05, which is a reasonably small arbitrary value with respect to the SNR levels tested. Note that for larger values of the Huber penalty will exhibit approximately quadratic behavior. was computed through noiseless reconstrucThe LLPR tions using (14). For a quantitative comparison, we computed the recovery coefficient (RC) [48] and full-width half-maximum (FWHM) of the perturbation responses. The RC is defined as for the LPR , and as for the LLPR . The FWHM was calculated as an angularly averaged radius of a contour line corresponding to half the maximum function value of the perturbation response by using the of Image Reconstruction Toolbox [49]. Fig. 3 shows that the RC and FWHM predicted by the LLPR agrees well with those of the LPR (denoted by “measured”) for high SNR (total count M), and they become more inaccurate as SNR decreases M). The error bars represent the standard errors (total count estimated by a bootstrap method. We also estimate the LLPR based on the LLIR using (35), denoted by “conventional.” The LLIR does not predict the LPR accurately even for high SNR, as shown in Fig. 3. Finally, we estimated the LPR by generating 500 realizations of white noise and implementing the efficient MC simulation method using local estimation with an 11 11 pixel region of support as described in Section III-B. The efficient MC simulation method predicts LPR very accurately even for low SNR, as shown in Fig. 3. We also compared the sample standard deviation at the center pixel estimated from the full MC simulation with the predicted values when the perturbation is present and when absent. As shown in Fig. 4, the efficient MC simulation method using local
Fig. 3. Comparison of measured and predicted recovery coefficients in (a) and FWHMs of LPR in (b) for the circular perturbation in the center of the phantom.
estimation described in Section III-B predicts the measured values accurately, but the conventional method using (48) underestimates the standard deviation, that is, it overestimates the smoothing, particularly for low SNR. To summarize, the LLPR, which can be computed efficiently by the local estimation method, predicts the resolution accurately, particularly for high SNR. The conventional method based on first-order Taylor series truncation is inaccurate in predicting the resolution and standard deviation for low SNR. However, the efficient MC simulation method using local estimation predicts both resolution and standard deviation accurately. Next, we examine the degree of computational saving we obtain in the efficient MC method using local estimation instead of full image reconstruction. We compared computation times for a full image reconstruction using noisy sinogram data and for the local estimation method to solve (46) with an 11 11 pixel region of support using a realization of white noise. The total photon count was 1 M and the Huber penalty was used and . We make a fair comparison of the with
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
421
Fig. 5. Comparison of convergence rates of the (OS)-SPS algorithm for full image reconstruction and the trust region reflective Newton method by the lsqnonlin function of Matlab for local estimation. The solution space for full and that for local estimation was . image reconstruction was
Fig. 4. Comparison of measured and predicted standard deviations of the center pixel (a) when the perturbation is absent and (b) when it is present in the center of the phantom.
convergence rates as follows. First, a solution of each optimization problem is determined. For full image reconstruction, 5000 iterations of the convergent SPS algorithm were run. For and local estimation, the termination criterion parameters, , of the function were reduced until the solution did not change; then several additional iterations of the steepest descent algorithm with a bisection line-search were run to ensure convergence. Next, for each method, the normalized between an iterate and the sodistance lution was computed as a function of computation time. All computation was performed on an AMD Opteron 870 2.0-GHz computer. As shown in Fig. 5, the local estimation is much faster than full image reconstruction as a result of reducing the soluto . tion space from C. Plug-In Method We next examined the effects of using estimated from a single noisy dataset in the plug-in method described in
Section III-C. We compared the LLPR estimated using known and values and that estimated using estimates . The circular perturbation was located in the center of the background, the total photon count was 500 K, and the Huber and was used. The estimates penalty with were obtained by the generalized error lookup table method [45]. We also compared predicted standard deviations. The LLPR was estimated by the local estimation method with an 11 11 region of support and the standard deviation was estimated by the efficient MC method using local estimation. Fig. 6 shows profiles through the LLPR and standard deviation image estimated using known ’s and using noisy estimates ’s. In Fig. 6(a), the error bars are not shown because they are very small. It can be seen that the LLPR and standard deviation predicted using the ’s agree well with those estimated from true values. As we move far from the center, the degree of disagreement increases. This is because we use a finite region of support with local estimation. If we use a larger region of support, we can reduce the inaccuracy with increased computational cost. D. Efficient Prediction for Shift-Invariant Systems As discussed in Section III-D, for a shift-invariant imaging system under the locally uniform background assumption, given perturbation and regularization parameters and , the LPR and standard deviation can be precomputed and tabulated as a function of , regardless of the perturbation location and background value. To investigate this, we used a circular perturbawith . We first precomputed the LPR and tion standard deviation for a range of values (0.004–0.04) by the efficient MC method using local estimation with a 11 11 reand gion of support. The Huber penalty with was used. In order to evaluate the predictability of the precomputed values, we computed the LPR and standard deviation for different perturbation locations (the center of the background, hot disc, and cold disc) and for different photon counts (0.25, 0.5, 1, 2, 4, and 8 millions) by the MC method using full noisy image reconstructions. As shown in Fig. 7, the recovery coefficients and FWHMs from the precomputed LPRs agree well
422
Fig. 6. Horizontal profiles of (a) the LLPR and (b) standard deviation image ^. estimated using true values and using noisy estimates
with the estimated ones obtained from the full MC simulation. The error bars are not shown in the figure because they are very small. The precomputed standard deviations also agree reasonably well with the measured ones, as shown in Fig. 8. E. Spatially Uniform Perturbation Response We evaluated the performance of the certainty-based modified penalty in terms of its ability to produce a spatially uniform perturbation response as described in Section III-E. We comand pared the ordinary Huber penalty with and the modified Huber penalty, based on noisy estimates , given in (49). The total photon count was 1M. We tested for and difcircular perturbations with different magnitudes ferent locations (the center of the background, hot disc, and cold disc). We compared the recovery coefficients of the LPRs for different perturbation locations estimated from MC simulations using full image reconstruction from 500 noisy datasets in each , of the case. Fig. 9 shows plots of the RC, i.e., ; error bars are LPR versus the perturbation magnitude not shown because they are very small. As shown in Fig. 9(a),
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
Fig. 7. (a) Comparison of recovery coefficients predicted by the efficient MC method using local estimation for a range of values and measured ones for different locations of the perturbation and for various photon counts. (b) Same as (a) except that we show the FWHMs of the LPR.
the ordinary Huber penalty yields space-variant perturbation responses because of space-variant maps. However, the modified Huber penalty generates spatially more uniform responses as shown in Fig. 9(b) although some degree of nonuniformity remains. V. DISCUSSION We derived new approximate expressions for the LPR (or LLPR) and variance for nonquadratically regularized image reconstruction. The expressions can be used for predicting resolution and variance accurately and efficiently. Under certain assumptions such as shift-invariant imaging systems and locally uniform background, one can precompute and tabulate the resolution and standard deviation as a function of and regularization parameters (e.g., and ), given a perturbation signal. For shift-variant imaging systems, one could precompute such a lookup table for a subset of all the voxels as in [5], [38]. Although the MC method using local estimation is much more efficient than the conventional MC method
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
Fig. 8. Comparison of standard deviation of the center pixel of a perturbation predicted by the efficient MC method using local estimation for a range of values and measured ones for different perturbation locations and for various photon counts (a) when the perturbation is present and (b) when it is absent. The number of photons follows the same pattern as in Fig. 7.
using full image reconstruction, it is still costly to predict variance (and resolution) particularly for shift-variant imaging systems; further investigation is needed for accelerating the computation. For local estimation, we determined the region of support empirically. It would be possible to determine it, for a nonquadratic penalty , systematically from the correlation , map of a quadratic penalty with which gives an upper bound for the spread of the perturbation responses. We achieved spatially uniform perturbation responses by extending the certainty-based modified penalty proposed in [1]. Although the modified penalty does not yield perfect spatial uniformity, it improves the spatial uniformity of the perturbation re, sponses. If a look-up table for resolution as a function of and were available, then one would be able to choose local regularization parameters ( and ) depending on estimated values, as in [5], to improve uniformity.
423
Fig. 9. Comparison of recovery coefficients of LPR (a) for the Huber penalty and (b) for the modified Huber penalty. This figure shows the RC of the LPR when the circular perturbations are located in the center of the background, hot disc, and cold disc for different perturbation magnitudes.
Future work will include further investigation of other types of penalties and perturbations with respect to image quality, lesion detectability and quantification accuracy. Interestingly, Qi has recently reported from simulation studies that nonquadratic regularization does not work better than quadratic regularization in the context of lesion detection and region-of-interest quantification [50]. Further investigation into the performance of nonquadratic regularization is still needed. ACKNOWLEDGMENT The authors would like to thank Dr. J. Fessler for his insightful comments and constructive suggestions. REFERENCES [1] J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized-likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process., vol. 5, no. 9, pp. 1346–1358, Sep. 1996.
424
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
[2] J. A. Fessler, “Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography,” IEEE Trans. Image Process., vol. 5, no. 3, pp. 493–506, Mar. 1996. [3] J. A. Fessler, “Approximate variance images for penalized-likelihood image reconstruction,” in Proc. IEEE Nucl. Sci. Symp. Med. Imag. Conf., 1997, vol. 2, pp. 949–952. [4] J. Qi and R. M. Leahy, “A theoretical study of the contrast recovery and variance of MAP reconstructions for PET data,” IEEE Trans. Med. Imag., vol. 18, no. 4, pp. 293–305, Apr. 1999. [5] J. Qi and R. M. Leahy, “Resolution and noise properties MAP reconstruction for fully 3-D PET,” IEEE Trans. Med. Imag., vol. 19, no. 5, pp. 493–506, May 2000. [6] J. W. Stayman and J. A. Fessler, “Efficient calculation of resolution and covariance for penalized-likelihood reconstruction in fully 3-D SPECT,” IEEE Trans. Med. Imag., vol. 23, no. 12, pp. 1543–1556, Dec. 2004. [7] E. Asma and R. M. Leahy, “Mean and covariance properties of dynamic PET reconstructions from list-mode data,” IEEE Trans. Med. Imag., vol. 25, no. 1, pp. 42–54, Jan. 2006. [8] Y. Zhang-O’Connor and J. A. Fessler, “Fast predictions of variance images for fan-beam transmission tomography with quadratic regularization,” IEEE Trans. Med. Imag., vol. 26, no. 3, pp. 335–346, Mar. 2007. [9] E. Veklerov and J. Llacer, “Stopping rule for the MLE algorithm based on statistical hypothesis testing,” IEEE Trans. Med. Imag., vol. 6, no. 4, pp. 313–319, Dec. 1987. [10] J. A. Fessler and S. D. Booth, “Conjugate-gradient preconditioning methods for shift-variant PET image reconstruction,” IEEE Trans. Image Process., vol. 8, no. 5, pp. 688–99, May 1999. [11] P. Bonetto, J. Qi, and R. M. Leahy, “Covariance approximation for fast and accurate computation of channelized Hotelling observer statistics,” IEEE Trans. Nucl. Sci., vol. 47, no. 4, pp. 1567–1572, Aug. 2000. [12] J. Qi and R. H. Huesman, “Theoretical study of lesion detectability of MAP reconstruction using computer observers,” IEEE Trans. Med. Imag., vol. 20, no. 8, pp. 815–822, Aug. 2001. [13] Y. Xing, I.-T. Hsiao, and G. Gindi, “Rapid calculation of detectability in Bayesian single photon emission computed tomography,” Phys. Med. Biol., vol. 48, no. 22, pp. 3755–3773, Nov. 2003. [14] J. Qi, “Analysis of lesion detectability in Bayesian emission reconstruction with nonstationary object variability,” IEEE Trans. Med. Imag., vol. 23, no. 3, pp. 321–329, Mar. 2004. [15] P. Khurd and G. Gindi, “Fast LROC analysis of Bayesian reconstructed emission tomographic images using model observers,” Phys. Med. Biol., vol. 50, pp. 1519–1532, Apr. 2005. [16] P. Khurd and G. Gindi, “Rapid computation of LROC figures of merit using numerical observers (for SPECT/PET reconstruction),” IEEE Trans. Nucl. Sci., vol. 52, no. 3, pp. 618–626, Jun. 2005. [17] A. Yendiki and J. A. Fessler, “Analysis of observer performance in known-location tasks for tomographic image reconstruction,” IEEE Trans. Med. Imag., vol. 25, no. 1, pp. 28–41, Jan. 2006. [18] J. Qi and R. H. Huesman, “Theoretical study of penalized-likelihood image reconstruction for region of interest quantification,” IEEE Trans. Med. Imag., vol. 25, no. 5, pp. 640–648, May 2006. [19] J. Qi and R. H. Huesman, “Penalized maximum-likelihood image reconstruction for lesion detection,” Phys. Med. Biol., vol. 51, no. 16, pp. 4017–4029, Aug. 2006. [20] S. Geman and D. E. McClure, “Bayesian image analysis: An application to single photon emission tomography,” in Proc. Stat. Comp. Sect. Amer. Stat. Assoc., 1985, pp. 12–18. [21] T. Hebert and R. M. Leahy, “A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors,” IEEE Trans. Med. Imag., vol. 8, no. 2, pp. 194–202, Jun. 1989. [22] P. J. Green, “Bayesian reconstructions from emission tomography data using a modified EM algorithm,” IEEE Trans. Med. Imag., vol. 9, no. 1, pp. 84–93, Mar. 1990. [23] K. Lange, “Convergence of EM image reconstruction algorithms with Gibbs smoothing,” IEEE Trans. Med. Imag., vol. 9, no. 4, pp. 439–446, Dec. 1990. [24] C. Bouman and K. Sauer, “A generalized Gaussian image model for edge-preserving MAP estimation,” IEEE Trans. Image Process., vol. 2, no. 3, pp. 296–310, Jul. 1993. [25] P. Charbonnier, L. Blanc-Féraud, G. Aubert, and M. Barlaud, “Deterministic edge-preserving regularization in computed imaging,” IEEE Trans. Image Process., vol. 6, no. 2, pp. 298–311, Feb. 1997.
[26] J. Nuyts, D. Beque, P. Dupont, and L. Mortelmans, “A concave prior penalizing relative differences for maximum-a-posteriori reconstruction in emission tomography,” IEEE Trans. Nucl. Sci., vol. 49, no. 1, pp. 56–60, Feb. 2002. [27] D. S. Lalush and B. M. W. Tsui, “Simulation evaluation of Gibbs prior distribution for use in maximum a posteriori-EM SPECT reconstructions,” IEEE Trans. Med. Imag., vol. 11, no. 2, pp. 267–275, Jun. 1992. [28] D. S. Lalush and B. M. W. Tsui, “A generalized Gibbs prior for maximum a posteriori reconstruction in SPECT,” Phys. Med. Biol., vol. 38, no. 6, pp. 729–741, Jun. 1993. [29] A. J. Rockmore and A. Macovski, “A maximum likelihood approach to emission image reconstruction from projections,” IEEE Trans. Nucl. Sci., vol. 23, pp. 1428–1432, 1976. [30] M. E. Casey and E. J. Hoffman, “Quantitation in positron emission computed tomography: 7 a technique to reduce noise in accidental coincidence measurements and coincidence efficiency calibration,” J. Comput. Assist. Tomogr., vol. 10, no. 5, pp. 845–850, 1986. [31] J. M. Ollinger, “Model-based scatter correction for fully 3-D PET,” Phys. Med. Biol., vol. 41, no. 1, pp. 153–76, Jan. 1996. [32] A. H. Delaney and Y. Bresler, “Globally convergent edge-preserving regularized reconstruction: An application to limited-angle tomography,” IEEE Trans. Image Process., vol. 7, no. 2, pp. 204–221, Feb. 1998. [33] H. L. Royden, Real Analysis, 3rd ed. Englewood Cliffs, NJ: PrenticeHall, 1988. [34] I. A. Elbakri and J. A. Fessler, “Statistical image reconstruction for polyenergetic X-ray computed tomography,” IEEE Trans. Med. Imag., vol. 21, no. 2, pp. 89–99, Feb. 2002. [35] W. Chlewicki, F. Hermansen, and S. B. Hansen, “Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior,” Phys. Med. Biol., vol. 49, pp. 4717–4730, 2004. [36] Y. Du, B. M. W. Tsui, and E. C. Frey, “Partial volume effect compensation for quantitative brain SPECT imaging,” IEEE Trans. Med. Imag., vol. 24, no. 8, pp. 969–976, Aug. 2005. [37] G. Boening, P. H. Pretorius, and M. A. King, “Study of relative quantification of Tc-99 m with partial volume effect and spillover correction for SPECT oncology imaging,” IEEE Trans. Nucl. Sci., vol. 53, no. 3, pp. 1205–1212, Jun. 2006. [38] J. W. Stayman and J. A. Fessler, “Compensation for nonuniform resolution using penalized-likelihood reconstruction in space-variant imaging systems,” IEEE Trans. Med. Imag., vol. 23, no. 3, pp. 269–284, Mar. 2004. [39] J. W. Stayman and J. A. Fessler, “Regularization for uniform spatial resolution properties in penalized-likelihood image reconstruction,” IEEE Trans. Med. Imag., vol. 19, no. 6, pp. 601–615, June 2000. [40] S. Ahn and R. M. Leahy, “Spatial resolution properties of nonquadratically regularized image reconstruction for PET,” in Proc. IEEE Intl. Symp. Biomed. Imag., 2006, pp. 287–290. [41] S. Ahn and J. A. Fessler, “Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms,” IEEE Trans. Med. Imag., vol. 22, no. 5, pp. 613–626, May 2003. [42] C. R. Rao, Linear Statistical Inference and Its Applications. New York: Wiley, 1973. [43] P. Boxandall and H. Liebeck, Vector Calculus. Oxford, U.K.: Clarendon, 1986. [44] P. J. Huber, Robust Statistics. New York: Wiley, 1981. [45] Q. Li, E. Asma, J. Qi, J. R. Bading, and R. M. Leahy, “Accurate estimation of the Fisher information matrix for the PET image reconstruction problem,” IEEE Trans. Med. Imag., vol. 23, no. 9, pp. 1057–1064, Sep. 2004. [46] J. A. Fessler, ASPIRE 3.0 user’s guide: A sparse iterative reconstruction library Commun. Sign. Process. Lab., Dept. EECS, Univ. Michigan, Ann Arbor, Tech. Rep. 293, Jul. 1995. [47] H. Erdo˘gan and J. A. Fessler, “Monotonic algorithms for transmission tomography,” IEEE Trans. Med. Imag., vol. 18, no. 9, pp. 801–814, Sep. 1999. [48] E. J. Hoffman, S. C. Huang, and M. E. Phelps, “Quantitation in positron emission computed tomography: 1 Effect of object size,” J. Comput. Assist. Tomogr., vol. 3, no. 3, pp. 299–308, 1979. [49] J. A. Fessler, Image reconstruction toolbox [Online]. Available: http:// www.eecs.umich.edu/~fessler [50] J. Qi, “Comparison of lesion detection and quantification in MAP reconstruction with Gaussian and non-Gaussian priors,” Int. J. Biomed. Imag., 2006 [Online]. Available: http://www.hindawi.com/GetArticle. aspx?doi=10.1155/IJBI/2006/87567
413
Analysis of Resolution and Noise Properties of Nonquadratically Regularized Image Reconstruction Methods for PET Sangtae Ahn, Member, IEEE, and Richard M. Leahy*, Fellow, IEEE
Abstract—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. Index Terms—Image resolution, noise, positron emission tomography (PET), regularization, statistical image reconstruction.
I. INTRODUCTION
S
TATISTICAL image reconstruction methods such as penalized-likelihood (PL) or maximum a posteriori (MAP) estimation are widely used for positron emission tomography (PET). The regularization or penalty function regularizes the ill-conditioned tomographic inverse problems. In other words, it controls a trade-off between resolution and noise in reconstructed images. The advantage of using penalty functions is that one can often predict accurately and/or efficiently the resolution and noise properties of reconstructed images and can design penalty functions for desired properties [1]–[8]. In contrast, other methods of stabilizing noise, such as early termination of
Manuscript received July 24, 2007; revised October 15, 2007. This work was supported by the National Institute of Biomedical Imaging and Bioengineering under Grant R01 EB00363. Asterisk indicates corresponding author. S. Ahn is with the Signal and Image Processing Institute, University of Southern California, Los Angeles, CA 90089 USA. *R. M. Leahy is with the Signal and Image Processing Institute, University of Southern California, Los Angeles, CA 90089 USA (e-mail: [email protected]. edu). Digital Object Identifier 10.1109/TMI.2007.911549
unregularized expectation maximization (EM) or ordered subsets EM (OSEM) [9], have only limited control over resolution and noise properties through choice of the number of iterations and subsets and post-reconstruction smoothing. Quadratically regularized estimators are approximately linear for reasonably high signal-to-noise ratio (SNR). Therefore, the resolution and variance of reconstructions can be estimated analytically without expensive Monte Carlo (MC) simulation [1], [2]. Furthermore, the use of Fourier transforms that take advantage of (approximate) Toeplitz structures in the system resolution matrices [10] enables one to efficiently predict resolution and variance as well as derived figures of merit for specific tasks such as detection and quantification and also to optimize regularization parameters for these tasks [4], [11]–[19]. However, the quadratic penalty tends to over-smooth sharp edges in reconstructed images. Many alternative nonquadratic penalty functions with edge-preserving properties have been proposed [20]–[26]. However, there has been little research on quantitative analysis of the resolution and noise properties for nonquadratic regularization in contrast to the rich literature on analysis and design of quadratic regularization. The lack of analytical results for nonquadratic regularization limits our ability to control the properties of reconstructed images through systematic selection of regularization parameters. Lalush and Tsui [27], [28] analyzed the smoothness in each iterate of the one-step-late (OSL) MAP-EM algorithm [22] in terms of the previous iterate by using the first derivative of a nonquadratic penalty function; however, the analysis does not apply to a case where iteration-independent and algorithm-independent MAP images are of interest. In this paper, we derive new expressions for estimating resolution and variance accurately and efficiently for nonquadratic regularization. In Section II, we define a local perturbation response (LPR), a generalization of the local impulse response (LIR) [1], to quantify resolution. We derive a new approximate expression in Section III-A for a linearized local perturbation response (LLPR), which approximates the LPR well for high SNR. While the LLPR is only implicitly expressed, it does accurately capture the nonlinear functional dependence of the perturbation response, which is intrinsic to nonquadratic regularization. In contrast, the conventional linearized local impulse response (LLIR) [1] is useful only when the superposition principle or linearity holds, which is not the case for nonquadratic regularization even for very high SNR. Although the expression for LLPR does not have a closed-form solution, we can solve the equation efficiently using a local estimation strategy as described in Section III-A. The expression allows us to achieve
0278-0062/$25.00 © 2007 IEEE
414
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
spatially uniform perturbation responses, akin to uniform spatial resolution in quadratic regularization, by using the certaintybased modified penalty [1], [5], as discussed in Section III-E. We also derive new expressions for the variance in the context of nonquadratic regularization in Section III-B. The implicit expression gives accurate estimates of the LPR and variance, even for low SNR, through efficient MC simulation. In Section III-C we describe a plug-in method that can be used when only a single noisy dataset is available and, in Section III-D, the use of a table lookup method to estimate resolution and variance. We perform simulation studies using a 2-D shift-invariant PET system to validate the methods in Section IV, and present concluding remarks in Section V. II. BACKGROUND
of horizontally or vertically neighboring pixels, for diagonally neighboring pixels, and 0 otherwise. Many types of nonconvex or nondifferentiable potential functions have been proposed [20], [21], [24], [32]. However, here we consider only convex and differentiable potential functions (e.g., [22], [23], [26]) because 1) our analysis requires the regularization function to be differentiable, 2) nonconvex potential functions may lead to unstable estimates in the sense that the estimator is not a continuous function of the data [24], and 3) it is more difficult to optimize nonconvex or nondifferentiable objective functions. Additionally, we assume the first derivative is absolutely continuous and therefore the second derivative exists almost everywhere [33, p. 109]. pairs of neighboring pixels Suppose there are a total of such that for . Define a differencing matrix as
A. Statistical Image Reconstruction is defined as A PL or MAP estimate a maximizer of an objective function as follows: (1) where denotes the PET emission data, is the number of pixels, is the number of detector pairs, and denotes matrix or vector transpose. The objective function is given by
and define for . The regularization function in (5) can then be rewritten as
The Hessian of
can be conveniently written as
(2) (6) where is the log-likelihood, is the regularization function or the energy function for a Gibbs prior, and is a global regularization parameter also known as a hyperparameter. Assuming the data follow the Poisson distribution [29], the log-likelihood can be written, neglecting constants, as
is the mean or noiseless data. The mean is affine in the object
denotes a diagwhich exists almost everywhere, where onal matrix that is formed by taking the arguments as diagonal elements. The simplest and most widely-used penalty is , which penalizes the squared difference of neighboring pixel values and consequently encourages smoothness in reconstructed images. However, the quadratic penalty tends to over-smooth sharp edges. We consider the following Huber penalty [34], [35] as a representative edge-preserving nonquadratic penalty, although the methods in this paper also apply to other penalties:
(4)
(7)
(3) where measurement
denotes the system or projection matrix charwhere acterizing the PET scanner and is the mean contribution of background events such as randoms and scatters. In this paper, we assume has already been estimated (for example, see [30] and [31]). For the regularization function, we focus on pairwise roughness penalty functions of the following form:
for some . The Huber penalty function quadratically penalizes small differences in neighboring pixels while applying a linear penalty to larger differences, with the transition determined by the parameter . Fig. 1 illustrates the quadratic and the Huber potential functions and their first and second derivatives. B. LPR
(5) where are nonnegative weights and is a potential function. and are symmetric, that is, We assume that and . Usually, in 2-D, is 1 for a pair
We define a LPR by which we quantify resolution properties. consists of a background Suppose that a true object and a signal of interest (8)
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
415
the resolution properties of nonquadratic regularization through the LPR as a function of perturbation . III. THEORY
Fig. 1. Quadratic and Huber potential functions and their first and second derivatives.
In Sections III-A and III-B, we develop methods which and require that the background and perturbation images, , are known. Such methods are useful for investigating the effects of the choice of penalty and its parameters on image quality metrics such as resolution, noise, lesion detectability, and quantification accuracy [4], [11]–[19] for various kinds of perturbations such as lesions of different sizes and contrasts. In Section III-C, we describe plug-in methods, which can be used when only a single noisy dataset is available. We then describe simplified methods for shift-invariant imaging systems in Section III-D and present a method for achieving spatially uniform perturbation responses in Section III-E. The analytical results in Sections III-A and III-E were partly presented in [40]. A. Asymptotic Resolution Analysis
For example, may represent a lesion or an organ of interest. We are interested in the effect of regularization upon recovery of the signal . In other words, we aim to study what would look like in a reconstructed image for a specific choice of regas the mean perturbation in the ularizer. We define a LPR reconstructed images caused by (9) where and are reconstructed images when the signal present and when absent, respectively. That is
(13) where is given in (1), and are defined in (12), and and in (10). By substituting (13) into (9), we define the linearized local perturbation response (LLPR) as
is (14) (10)
where when
1) LLPR: The LPR defined in (9) is intractable to analyze due to the mean operators. Fortunately, the mean of an estimator can be approximated by a noiseless estimate for high SNR1 [2]
is given in (1), and and represent the noisy data is present and absent, respectively. In other words (11)
where (12) ” In the equations above, is given in (4), and “ means that the elements of are independent Poisson random . Similar types of regional perturbation variables and responses have been used for analyzing and compensating for partial volume effects (e.g., see [36] and [37]). The LPR is a generalization of point, line, and edge spread is of a particular type. If is an functions where the signal impulse with an infinitesimal amplitude, then the LPR reduces to the LIR [1]. LIR is useful for studying resolution properties in quadratic regularization, which is approximately linear [5], [38] and thus the superposition principle holds. The LIR has been used as a key tool for quadratic penalty design methods [1], [4], [5], [39], [38]. However, in nonquadratic regularization, resoluand LIR tion properties depend heavily on the perturbation is of limited utility. For nonquadratic regularization, LIR based prediction of resolution and variance is inaccurate, particularly when SNR is low, as shown in Section IV. Therefore, we study
akin to the LLIR [1]. It is called a linearized response because the approximation is based on the first-order Taylor series truncation of the estimator [1]. That is, if , then . In Section IV-B, we show how well LLPR approximates LPR for various SNR levels. is to A straightforward method for computing the LLPR and using noiseless data and and then reconstruct to use (14); however, this gives little insight in the sense that it and . We is difficult to examine the relationship between instead derive an approximate expression that relates to the true perturbation . The new expression suggests a computationally efficient strategy and provides useful insights into resolution properties. Neglecting the nonnegativity constraint, which is reasonable for noiseless reconstructions [1], the noiseless estimates, and , satisfy the first-order optimality conditions (15) (16) 1Under mild conditions including convexity of and P y 6= 0, the objective function 8 in (1) is strictly concave [41, Lemma 1]. It follows that x ^( 1 ) in (1) is continuous [24, Theorem 1]. Note that if y ’s are iid random vectors, then y =K converges in probability to E [y ] by the weak law ^, x^( y =K ) of large numbers [42, p. 112]. Due to the continuity of x ^(E [y ]) [42, p. 124]. Therefore, as the photon converges in probability to x count increases (K ! 1), the estimate converges in probability to the noiseless estimate.
416
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
where is the column gradient operator with respect to the first argument. Consider the first-order Taylor series expansion about of
where . Note the second derivative term disappears in the above equation. Combining (4) and (15)–(19) yields
(17)
(22) (23)
is the Hessian operator with respect to the first where th element of the operator is argument, and the . The above approximation is reasonably accurate is linear in and varies slowly with because through its projection [see (3)]. However, in (16) may not be well approximated by a Taylor series truncation is nonlinear. For example, if is nonquadratic and thus whereas Fig. 1 shows the nonlinear behavior of is perfectly linear. Therefore, we take the zeroth-order Taylor expansion with the remainder term or equivalently the fundamental theorem of calculus [43, p. 292] (18)
Although (22) may apply to other system models and noise models, here we focus on the linear emission tomography system with the Poisson noise model. Using (3) and (4), one can compute the following derivative terms: (24) (25) where mating
is the Fisher information matrix for estigiven by (26)
Since
The approximation in (24) is reasonably accurate because 1) is a blurred version of, and hence very similar to, and 2) , and are all smoothing operators [1]. Substituting (4), (24), and (25) into (23) leads to (27) by (6), one can rewrite (18) as (19) where
This is a new expression that relates the LLPR to the true perturbation . To simplify the above expression, we consider the following factored form of the system matrix:
is defined as
(28) (20)
with (21) Interestingly, the above form of has appeared in different to design a contexts. For example, Huber has used quadratic surrogate for the penalty function with the aim of constructing monotonic optimization algorithms [44, p. 184]. Charbonnier et al. have described certain conditions for edge-pre[25]. Also, note that (20) serving properties in terms of and (21) are similar to (36) in [4], which is an ad hoc expression to predict the contrast recovery when the Huber penalty is used. If the penalty function is quadratic, then is a constant maof the penalty. By (20) and (21), trix, which is the Hessian the second term on the right hand of (19) can be calculated as
where denote ray-dependent factors including detector efficiencies, attenuation, and dead time, and is the geometric projection matrix possibly including space-variant sinogram blurring [1], [4]. Then the Fisher information matrix in (26) can be approximated as
(29) where (30) (31) because is concentrated about its diagonal and , which is interpreted as a measure of certainty, varies slowly with [1]. The approximation in (29) separates into an object-dependent and a precomputable object-independent part . part
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
Suppose the perturbation is localized around pixel . For a will be similarly localreasonable degree of regularization, ized. This localized perturbation assumption allows us to further simplify the approximation for
(32) is concentrated about its diagonal and because slowly [1]. Similarly
varies
(33) Substituting (32) and (33) into (27) yields a simplification of (27) (34) (cf. [1, (34)]). The approximate expression (27) for the LLPR and its simplified version (34) are the key results of this subsection, which can be viewed as a generalization of the approximations for the LLIR in [1]. We now briefly compare the LLPR and LLIR approximations. If we used the first-order Taylor series truncation instead of the Taylor expansion with remainder for in (19), we would obtain the following explicit expression for : (35) where (36) The th column of is the LLIR for pixel [1]. If in (27) is replaced with , then (27) becomes (35) and (36). Thus one can view (35) and (36) as being obtained in (27). The approxiby discarding the dependency of on mation (35) provides an explicit expression for , which is reasonably accurate for quadratic regularization. However, it does not predict the LLPR accurately for nonquadratic regularization particularly in a low SNR case whereas (27) does, as shown in Section IV-B. 2) Efficient Computation Using Local Estimation: One can in (27) by solving the following nonlinear least compute squares problem:
417
However, solving (37) is no less costly than performing noiseby (14). less image reconstructions and computing is localized about pixel , then If the perturbation and will also be concentrated about pixel . Based on this observation, one can downsize the problem (37) using local estimation. First, we define a set of region of support about pixel for some pixels within an such that . Next, we restrict the parameter space to by taking only the pixels in the set as optimization variables for the least squares problem. Equivalently, we use as a weight matrix in (37) where is an indicator function on . Choosing an appropriate depends and the regularization parameters such as on the support of . As increases, the LLPR computed by the local estimation method will become more accurate but the computational cost will increase. This downsizing strategy is somewhat similar to the local backprojection technique used in efficient quadratic penalty design [38] or the small volume support technique used in efficient calculation of resolution and covariance for fully 3-D SPECT [6]. An advantage of computing the LLPR by (37) rather than using (14) through noiseless image reconstructions is that one can use (37) under certain conditions even when only a single noisy dataset is available, as discussed in Section III-C. B. Efficient Monte Carlo Simulation Method for Variance Analysis To estimate the variance of noisy reconstruction , defined in (10), one can use a MC simulation method where a number of noisy datasets are generated and the sample variance is computed from noisy reconstructions for each pixel. However, we present a more efficient MC simulation method here. Neglecting again the nonnegativity constraint, the following first-order optimality conditions are satisfied: (39) (40) where and are defined in in (11), (12), and (13), respectively. As in Section III-A, combining the first-order Taylor around and the zeroth-order series truncation of around with the remainder term Taylor expansion of yields
(41) where the deviation
of
from
is defined as
(37) (42) where
is a diagonal weight matrix and Computing the derivatives as in (24) and (25) and substituting them into (41) leads to
Similarly, one can solve a least squares problem based on (34) by using (38)
(43)
418
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
where
Note that
is approximately white noise, that is, because and
and
. For the factored system model in (28), the right-hand side in (43) can be approximated as
will be . However, the nonlinin in (43) or (46), then earity in nonquadratic regularization manifests itself in nonzero , causing the LLPR to deviate from the LPR particularly for low SNR. If we used the first-order Taylor series truncation for in (40) in place of the Taylor expansion with the remainder term, and the covariwe would obtain an explicit expression for ance of , equal to the covariance of , could be written as
(48) (44) since the means and covariances of both sides above approximately match in view of (28) and (29). Combining (29), (43), and (44) yields
which is equivalent to the covariance expression derived in [2]. For nonquadratic regularization, the above expression becomes inaccurate as the photon count decreases, as shown in Section IV-B. C. Plug-In Method
(45) Again, if our region of interest is localized around pixel , (45) can be locally approximated as (46) The approximate expression (43) for the deviation and its simplified version (46) are the main results of this subsection. in the neighborhood One can estimate the variance of of pixel using (43) or (46) as follows. First, one generates a number of realizations of white noise , and for each reby taking a small alization one solves (43) or (46) for region of support about pixel and using the local estimation method, as described in Section III-A-2. The variance of is estimated by the sample variance of in view of (42). This MC simulation method with local estimation is much less expensive than the ordinary MC simulation method with full image reconstruction, as shown in Section IV-B. The efficient MC simulation method with local estimation can also provide an accurate estimate of the LPR. As the photon count decreases, the accuracy of the LLPR ’s approximation to LPR becomes worse. However, using the MC results, one can predict the LPR accurately even for low SNR. One can rewrite the LPR , defined in (9), as follows:
(47) where . Therefore, can be viewed as a correction term when one predicts the LPR from the LLPR given in (14). Consider a special case where the background image is locally uniform. In this case, since locally and therefore , one can compute by adding , estimated by the sample mean from the MC simulation using (43) or (46), to , estimated by solving (37). This locally uniform background assumption is similar to is not locally uniform, one would [18, Approximation 1]. If from another Monte Carlo simulation have to estimate is linear using instead of in (43) or (46). Note that if
In Sections III-A and III-B, we derived approximate expressions (34) and (46) for the LLPR and variance, respectively, for a factored system model. It is straightforward to use those expressions in a simulation study where the background and perturbation images, and , are known. However, in practice, only a single noisy dataset is available. We address three issues that arise when the methods given in Sections III-A and III-B are to be applied to this practical case. First, one can estimate from the noisy data and substitute in (34) and (46) [1], [2], [4], [5]. The plug-in esthem for varies slowly with the image timate is fairly accurate since . We use the generalized error lookup table and method [45] to estimate in Section IV. Second, the expression in (34) requires , which is unavailable in practice. It is not acceptable to substitute a noisy is sensitive to a noisy estimate reconstruction because of and moreover it is not possible to estimate only the background image in the presence of the signal . So we is locally uniform in the assume the background image region of interest, as is , and we restrict the analysis using the plug-in method to the locally uniform background case. in Based on this assumption, one can replace for some where (34) with and , consequently removing the dependency on the background image. Also, in (46) can be approximated by since under the locally uniform background assumption. is needed in (34); however, it Third, the true perturbation is unknown in practice. But one can investigate the resolution and variance properties for various typical perturbations as in a signal-known-exactly detection task (e.g., [12], [17]). One could also establish a statistical model of the perturbation and study the properties of reconstructed images as in signal-known-statistically tasks (e.g., [14]). D. Prediction of Resolution and Variance for Shift-Invariant Imaging Systems For a localized perturbation , one can use (34) to estimate and use (46) to estimate deviations , from the LLPR
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
which the variance of can be estimated. The expressions (34) and (46) are functions of 1) the true perturbation , 2) the regularization parameters, e.g., and for the Huber penalty, 3) measures of certainty , 4) (the mean reconstruction of the) , background image , and 5) the system response matrix which is data independent. As we add more assumptions, we can further simplify the situation. Assume that the background image , and hence , are locally uniform in the region of inwith in (34) terest. Then one can replace with in (46) as discussed and replace are no longer functions of in Section III-C, and now and the background image. Additionally, assume that the system is and will also be shift-invariant. Then (34) and (46) for shift-invariant, and the perturbation response and its statistical properties are not dependent on the location of the perturbation in the image domain. To summarize, if the background image is locally uniform and the system is shift-invariant, then for a given perturbation and regularization parameters and , one can precompute the LLPR (or LPR) and variance for a range of different ’s. Then the LLPR (or LPR) and variance in practice can be esregardless of the timated by table lookup using an estimate perturbation location and the background image. A similar table lookup strategy has been proposed in [3], although for quadratic regularization only.
419
Fig. 2. Digital phantom used in simulation studies. The background, the left hot disc and the right cold disc have relative emission activities of 2, 3, and 1, respectively. This figure also shows the circular perturbation of diameter 3 pixels in white that we used in the studies presented.
(34) in Section III-A, one can obtain the following approximawhen is used as a regularization function for the LLPR is localized about pixel : tion, assuming the perturbation
(50) where (51) with being defined as the th pair of neighboring is also localized about pixel , the following pixels. Since approximation is valid:
E. Spatially Uniform Perturbation Response
(52)
For nonquadratic regularization, the resolution is a function of the perturbation. When the LLPR , which approximates the LPR well for high SNR, is viewed as a function of is shift-variant even for a perturbation , the system shift-invariant imaging system partly due to the contribution of the background image through the ’s. In this subsection, we make the LLPR system shift-invariant using appropriate penalty weights. This is an extension of the quadratic penalty design for uniform spatial resolution [1], [5] to a nonquadratic penalty case. Note that the nonquadratic penalty is intended to yield nonuniform resolution and consequently what we aim to achieve is not uniform spatial resolution but rather a spatially uniform functional dependence of the LLPR on the perturbation. Here, we assume the imaging system is shift-invariant and the background is locally uniform. Using the implicit expression for the LLPR in (34), we show that the spatially uniform perturbation response of the LLPR can be achieved by using the following certainty-based modified penalty, which is proposed in [1]: (49) are estimated from noisy data, and and are the where same as given in (5). Following a similar approach to deriving
where
is defined in (20). Substituting (52) into (50) yields (53)
and we assume the background image is since locally uniform about pixel . Therefore, using the modified . penalty in (49) removes the dependence of LLPR on It can be seen in (53) that if one considers a system whose input is perturbation and whose output is LLPR , then the system is shift-invariant. Therefore, the modified penalty achieves the spatially uniform perturbation response properties. IV. RESULTS A. Two-Dimensional PET Simulation To validate the methods developed in Section III, 2-D PET emission scans were simulated. The synthetic emission phantom used in this study had relative emission activities of 1, 2, and 3 in the right cold disc, the warm elliptic background, and the left hot disc, respectively, as shown in Fig. 2. Since some of the methods we test are for shift-invariant imaging systems from Sections III-D and III-E, we generated a nearly shift-invariant 2-D PET system matrix using ASPIRE [46] to simplify the validation. The sinograms had 192 radial bins and 160 angles uniformly sampled over 180 . Nonuniform detector efficiencies were considered, and attenuation factors, from a uniform attenuation coefficient 0.095 cm , were included. The background
420
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
contribution factors corresponded to a uniform field of 10% of the total count. Reconstructed images were obtained, unless stated otherwise, using 20 iterations of fast ordered subset separable paraboloidal surrogates (OS-SPS) with 16 subsets [41] followed by 20 iterations of convergent SPS with optimal curvature [47] to ensure convergence, with filtered backprojection (FBP) initialization. The image had 128 128 pixels and a first-order neighborhood system was used. For a perturbation , a uniform circle of diameter 3 pixels was placed in the center of either the background, hot region or cold region, as shown in Fig. 2. For the local estimation described in Sections III, to solve the nonlinear least squares problems in (37), we used the function in the Optimization Toolbox of Matlab (Mathworks, Natick, MA). B. Accuracy of Prediction We evaluate the accuracy of the approximate expressions given in Sections III-A and III-B. The uniform circular activity located in the center of the background is taken as a perturbawith magnitude . Noisy datasets following tion the Poisson distribution were generated for various photon counts (0.25, 0.5, 1, 2, 4, and 8 millions) by scaling in (28) appropriately, in the presence and absence of the perturbation. was estimated by substituting the sample means, The LPR computed from 500 realizations, into the means in (9). The Huber penalty was used for reconstruction and the parameters and where is the total photon were set to count in millions in order to fix the resolution for different counts. The parameter was chosen to be 0.05, which is a reasonably small arbitrary value with respect to the SNR levels tested. Note that for larger values of the Huber penalty will exhibit approximately quadratic behavior. was computed through noiseless reconstrucThe LLPR tions using (14). For a quantitative comparison, we computed the recovery coefficient (RC) [48] and full-width half-maximum (FWHM) of the perturbation responses. The RC is defined as for the LPR , and as for the LLPR . The FWHM was calculated as an angularly averaged radius of a contour line corresponding to half the maximum function value of the perturbation response by using the of Image Reconstruction Toolbox [49]. Fig. 3 shows that the RC and FWHM predicted by the LLPR agrees well with those of the LPR (denoted by “measured”) for high SNR (total count M), and they become more inaccurate as SNR decreases M). The error bars represent the standard errors (total count estimated by a bootstrap method. We also estimate the LLPR based on the LLIR using (35), denoted by “conventional.” The LLIR does not predict the LPR accurately even for high SNR, as shown in Fig. 3. Finally, we estimated the LPR by generating 500 realizations of white noise and implementing the efficient MC simulation method using local estimation with an 11 11 pixel region of support as described in Section III-B. The efficient MC simulation method predicts LPR very accurately even for low SNR, as shown in Fig. 3. We also compared the sample standard deviation at the center pixel estimated from the full MC simulation with the predicted values when the perturbation is present and when absent. As shown in Fig. 4, the efficient MC simulation method using local
Fig. 3. Comparison of measured and predicted recovery coefficients in (a) and FWHMs of LPR in (b) for the circular perturbation in the center of the phantom.
estimation described in Section III-B predicts the measured values accurately, but the conventional method using (48) underestimates the standard deviation, that is, it overestimates the smoothing, particularly for low SNR. To summarize, the LLPR, which can be computed efficiently by the local estimation method, predicts the resolution accurately, particularly for high SNR. The conventional method based on first-order Taylor series truncation is inaccurate in predicting the resolution and standard deviation for low SNR. However, the efficient MC simulation method using local estimation predicts both resolution and standard deviation accurately. Next, we examine the degree of computational saving we obtain in the efficient MC method using local estimation instead of full image reconstruction. We compared computation times for a full image reconstruction using noisy sinogram data and for the local estimation method to solve (46) with an 11 11 pixel region of support using a realization of white noise. The total photon count was 1 M and the Huber penalty was used and . We make a fair comparison of the with
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
421
Fig. 5. Comparison of convergence rates of the (OS)-SPS algorithm for full image reconstruction and the trust region reflective Newton method by the lsqnonlin function of Matlab for local estimation. The solution space for full and that for local estimation was . image reconstruction was
Fig. 4. Comparison of measured and predicted standard deviations of the center pixel (a) when the perturbation is absent and (b) when it is present in the center of the phantom.
convergence rates as follows. First, a solution of each optimization problem is determined. For full image reconstruction, 5000 iterations of the convergent SPS algorithm were run. For and local estimation, the termination criterion parameters, , of the function were reduced until the solution did not change; then several additional iterations of the steepest descent algorithm with a bisection line-search were run to ensure convergence. Next, for each method, the normalized between an iterate and the sodistance lution was computed as a function of computation time. All computation was performed on an AMD Opteron 870 2.0-GHz computer. As shown in Fig. 5, the local estimation is much faster than full image reconstruction as a result of reducing the soluto . tion space from C. Plug-In Method We next examined the effects of using estimated from a single noisy dataset in the plug-in method described in
Section III-C. We compared the LLPR estimated using known and values and that estimated using estimates . The circular perturbation was located in the center of the background, the total photon count was 500 K, and the Huber and was used. The estimates penalty with were obtained by the generalized error lookup table method [45]. We also compared predicted standard deviations. The LLPR was estimated by the local estimation method with an 11 11 region of support and the standard deviation was estimated by the efficient MC method using local estimation. Fig. 6 shows profiles through the LLPR and standard deviation image estimated using known ’s and using noisy estimates ’s. In Fig. 6(a), the error bars are not shown because they are very small. It can be seen that the LLPR and standard deviation predicted using the ’s agree well with those estimated from true values. As we move far from the center, the degree of disagreement increases. This is because we use a finite region of support with local estimation. If we use a larger region of support, we can reduce the inaccuracy with increased computational cost. D. Efficient Prediction for Shift-Invariant Systems As discussed in Section III-D, for a shift-invariant imaging system under the locally uniform background assumption, given perturbation and regularization parameters and , the LPR and standard deviation can be precomputed and tabulated as a function of , regardless of the perturbation location and background value. To investigate this, we used a circular perturbawith . We first precomputed the LPR and tion standard deviation for a range of values (0.004–0.04) by the efficient MC method using local estimation with a 11 11 reand gion of support. The Huber penalty with was used. In order to evaluate the predictability of the precomputed values, we computed the LPR and standard deviation for different perturbation locations (the center of the background, hot disc, and cold disc) and for different photon counts (0.25, 0.5, 1, 2, 4, and 8 millions) by the MC method using full noisy image reconstructions. As shown in Fig. 7, the recovery coefficients and FWHMs from the precomputed LPRs agree well
422
Fig. 6. Horizontal profiles of (a) the LLPR and (b) standard deviation image ^. estimated using true values and using noisy estimates
with the estimated ones obtained from the full MC simulation. The error bars are not shown in the figure because they are very small. The precomputed standard deviations also agree reasonably well with the measured ones, as shown in Fig. 8. E. Spatially Uniform Perturbation Response We evaluated the performance of the certainty-based modified penalty in terms of its ability to produce a spatially uniform perturbation response as described in Section III-E. We comand pared the ordinary Huber penalty with and the modified Huber penalty, based on noisy estimates , given in (49). The total photon count was 1M. We tested for and difcircular perturbations with different magnitudes ferent locations (the center of the background, hot disc, and cold disc). We compared the recovery coefficients of the LPRs for different perturbation locations estimated from MC simulations using full image reconstruction from 500 noisy datasets in each , of the case. Fig. 9 shows plots of the RC, i.e., ; error bars are LPR versus the perturbation magnitude not shown because they are very small. As shown in Fig. 9(a),
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
Fig. 7. (a) Comparison of recovery coefficients predicted by the efficient MC method using local estimation for a range of values and measured ones for different locations of the perturbation and for various photon counts. (b) Same as (a) except that we show the FWHMs of the LPR.
the ordinary Huber penalty yields space-variant perturbation responses because of space-variant maps. However, the modified Huber penalty generates spatially more uniform responses as shown in Fig. 9(b) although some degree of nonuniformity remains. V. DISCUSSION We derived new approximate expressions for the LPR (or LLPR) and variance for nonquadratically regularized image reconstruction. The expressions can be used for predicting resolution and variance accurately and efficiently. Under certain assumptions such as shift-invariant imaging systems and locally uniform background, one can precompute and tabulate the resolution and standard deviation as a function of and regularization parameters (e.g., and ), given a perturbation signal. For shift-variant imaging systems, one could precompute such a lookup table for a subset of all the voxels as in [5], [38]. Although the MC method using local estimation is much more efficient than the conventional MC method
AHN AND LEAHY: ANALYSIS OF RESOLUTION AND NOISE PROPERTIES OF NONQUADRATICALLY REGULARIZED IMAGE RECONSTRUCTION METHODS
Fig. 8. Comparison of standard deviation of the center pixel of a perturbation predicted by the efficient MC method using local estimation for a range of values and measured ones for different perturbation locations and for various photon counts (a) when the perturbation is present and (b) when it is absent. The number of photons follows the same pattern as in Fig. 7.
using full image reconstruction, it is still costly to predict variance (and resolution) particularly for shift-variant imaging systems; further investigation is needed for accelerating the computation. For local estimation, we determined the region of support empirically. It would be possible to determine it, for a nonquadratic penalty , systematically from the correlation , map of a quadratic penalty with which gives an upper bound for the spread of the perturbation responses. We achieved spatially uniform perturbation responses by extending the certainty-based modified penalty proposed in [1]. Although the modified penalty does not yield perfect spatial uniformity, it improves the spatial uniformity of the perturbation re, sponses. If a look-up table for resolution as a function of and were available, then one would be able to choose local regularization parameters ( and ) depending on estimated values, as in [5], to improve uniformity.
423
Fig. 9. Comparison of recovery coefficients of LPR (a) for the Huber penalty and (b) for the modified Huber penalty. This figure shows the RC of the LPR when the circular perturbations are located in the center of the background, hot disc, and cold disc for different perturbation magnitudes.
Future work will include further investigation of other types of penalties and perturbations with respect to image quality, lesion detectability and quantification accuracy. Interestingly, Qi has recently reported from simulation studies that nonquadratic regularization does not work better than quadratic regularization in the context of lesion detection and region-of-interest quantification [50]. Further investigation into the performance of nonquadratic regularization is still needed. ACKNOWLEDGMENT The authors would like to thank Dr. J. Fessler for his insightful comments and constructive suggestions. REFERENCES [1] J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized-likelihood image reconstruction methods: Space-invariant tomographs,” IEEE Trans. Image Process., vol. 5, no. 9, pp. 1346–1358, Sep. 1996.
424
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 27, NO. 3, MARCH 2008
[2] J. A. Fessler, “Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography,” IEEE Trans. Image Process., vol. 5, no. 3, pp. 493–506, Mar. 1996. [3] J. A. Fessler, “Approximate variance images for penalized-likelihood image reconstruction,” in Proc. IEEE Nucl. Sci. Symp. Med. Imag. Conf., 1997, vol. 2, pp. 949–952. [4] J. Qi and R. M. Leahy, “A theoretical study of the contrast recovery and variance of MAP reconstructions for PET data,” IEEE Trans. Med. Imag., vol. 18, no. 4, pp. 293–305, Apr. 1999. [5] J. Qi and R. M. Leahy, “Resolution and noise properties MAP reconstruction for fully 3-D PET,” IEEE Trans. Med. Imag., vol. 19, no. 5, pp. 493–506, May 2000. [6] J. W. Stayman and J. A. Fessler, “Efficient calculation of resolution and covariance for penalized-likelihood reconstruction in fully 3-D SPECT,” IEEE Trans. Med. Imag., vol. 23, no. 12, pp. 1543–1556, Dec. 2004. [7] E. Asma and R. M. Leahy, “Mean and covariance properties of dynamic PET reconstructions from list-mode data,” IEEE Trans. Med. Imag., vol. 25, no. 1, pp. 42–54, Jan. 2006. [8] Y. Zhang-O’Connor and J. A. Fessler, “Fast predictions of variance images for fan-beam transmission tomography with quadratic regularization,” IEEE Trans. Med. Imag., vol. 26, no. 3, pp. 335–346, Mar. 2007. [9] E. Veklerov and J. Llacer, “Stopping rule for the MLE algorithm based on statistical hypothesis testing,” IEEE Trans. Med. Imag., vol. 6, no. 4, pp. 313–319, Dec. 1987. [10] J. A. Fessler and S. D. Booth, “Conjugate-gradient preconditioning methods for shift-variant PET image reconstruction,” IEEE Trans. Image Process., vol. 8, no. 5, pp. 688–99, May 1999. [11] P. Bonetto, J. Qi, and R. M. Leahy, “Covariance approximation for fast and accurate computation of channelized Hotelling observer statistics,” IEEE Trans. Nucl. Sci., vol. 47, no. 4, pp. 1567–1572, Aug. 2000. [12] J. Qi and R. H. Huesman, “Theoretical study of lesion detectability of MAP reconstruction using computer observers,” IEEE Trans. Med. Imag., vol. 20, no. 8, pp. 815–822, Aug. 2001. [13] Y. Xing, I.-T. Hsiao, and G. Gindi, “Rapid calculation of detectability in Bayesian single photon emission computed tomography,” Phys. Med. Biol., vol. 48, no. 22, pp. 3755–3773, Nov. 2003. [14] J. Qi, “Analysis of lesion detectability in Bayesian emission reconstruction with nonstationary object variability,” IEEE Trans. Med. Imag., vol. 23, no. 3, pp. 321–329, Mar. 2004. [15] P. Khurd and G. Gindi, “Fast LROC analysis of Bayesian reconstructed emission tomographic images using model observers,” Phys. Med. Biol., vol. 50, pp. 1519–1532, Apr. 2005. [16] P. Khurd and G. Gindi, “Rapid computation of LROC figures of merit using numerical observers (for SPECT/PET reconstruction),” IEEE Trans. Nucl. Sci., vol. 52, no. 3, pp. 618–626, Jun. 2005. [17] A. Yendiki and J. A. Fessler, “Analysis of observer performance in known-location tasks for tomographic image reconstruction,” IEEE Trans. Med. Imag., vol. 25, no. 1, pp. 28–41, Jan. 2006. [18] J. Qi and R. H. Huesman, “Theoretical study of penalized-likelihood image reconstruction for region of interest quantification,” IEEE Trans. Med. Imag., vol. 25, no. 5, pp. 640–648, May 2006. [19] J. Qi and R. H. Huesman, “Penalized maximum-likelihood image reconstruction for lesion detection,” Phys. Med. Biol., vol. 51, no. 16, pp. 4017–4029, Aug. 2006. [20] S. Geman and D. E. McClure, “Bayesian image analysis: An application to single photon emission tomography,” in Proc. Stat. Comp. Sect. Amer. Stat. Assoc., 1985, pp. 12–18. [21] T. Hebert and R. M. Leahy, “A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors,” IEEE Trans. Med. Imag., vol. 8, no. 2, pp. 194–202, Jun. 1989. [22] P. J. Green, “Bayesian reconstructions from emission tomography data using a modified EM algorithm,” IEEE Trans. Med. Imag., vol. 9, no. 1, pp. 84–93, Mar. 1990. [23] K. Lange, “Convergence of EM image reconstruction algorithms with Gibbs smoothing,” IEEE Trans. Med. Imag., vol. 9, no. 4, pp. 439–446, Dec. 1990. [24] C. Bouman and K. Sauer, “A generalized Gaussian image model for edge-preserving MAP estimation,” IEEE Trans. Image Process., vol. 2, no. 3, pp. 296–310, Jul. 1993. [25] P. Charbonnier, L. Blanc-Féraud, G. Aubert, and M. Barlaud, “Deterministic edge-preserving regularization in computed imaging,” IEEE Trans. Image Process., vol. 6, no. 2, pp. 298–311, Feb. 1997.
[26] J. Nuyts, D. Beque, P. Dupont, and L. Mortelmans, “A concave prior penalizing relative differences for maximum-a-posteriori reconstruction in emission tomography,” IEEE Trans. Nucl. Sci., vol. 49, no. 1, pp. 56–60, Feb. 2002. [27] D. S. Lalush and B. M. W. Tsui, “Simulation evaluation of Gibbs prior distribution for use in maximum a posteriori-EM SPECT reconstructions,” IEEE Trans. Med. Imag., vol. 11, no. 2, pp. 267–275, Jun. 1992. [28] D. S. Lalush and B. M. W. Tsui, “A generalized Gibbs prior for maximum a posteriori reconstruction in SPECT,” Phys. Med. Biol., vol. 38, no. 6, pp. 729–741, Jun. 1993. [29] A. J. Rockmore and A. Macovski, “A maximum likelihood approach to emission image reconstruction from projections,” IEEE Trans. Nucl. Sci., vol. 23, pp. 1428–1432, 1976. [30] M. E. Casey and E. J. Hoffman, “Quantitation in positron emission computed tomography: 7 a technique to reduce noise in accidental coincidence measurements and coincidence efficiency calibration,” J. Comput. Assist. Tomogr., vol. 10, no. 5, pp. 845–850, 1986. [31] J. M. Ollinger, “Model-based scatter correction for fully 3-D PET,” Phys. Med. Biol., vol. 41, no. 1, pp. 153–76, Jan. 1996. [32] A. H. Delaney and Y. Bresler, “Globally convergent edge-preserving regularized reconstruction: An application to limited-angle tomography,” IEEE Trans. Image Process., vol. 7, no. 2, pp. 204–221, Feb. 1998. [33] H. L. Royden, Real Analysis, 3rd ed. Englewood Cliffs, NJ: PrenticeHall, 1988. [34] I. A. Elbakri and J. A. Fessler, “Statistical image reconstruction for polyenergetic X-ray computed tomography,” IEEE Trans. Med. Imag., vol. 21, no. 2, pp. 89–99, Feb. 2002. [35] W. Chlewicki, F. Hermansen, and S. B. Hansen, “Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior,” Phys. Med. Biol., vol. 49, pp. 4717–4730, 2004. [36] Y. Du, B. M. W. Tsui, and E. C. Frey, “Partial volume effect compensation for quantitative brain SPECT imaging,” IEEE Trans. Med. Imag., vol. 24, no. 8, pp. 969–976, Aug. 2005. [37] G. Boening, P. H. Pretorius, and M. A. King, “Study of relative quantification of Tc-99 m with partial volume effect and spillover correction for SPECT oncology imaging,” IEEE Trans. Nucl. Sci., vol. 53, no. 3, pp. 1205–1212, Jun. 2006. [38] J. W. Stayman and J. A. Fessler, “Compensation for nonuniform resolution using penalized-likelihood reconstruction in space-variant imaging systems,” IEEE Trans. Med. Imag., vol. 23, no. 3, pp. 269–284, Mar. 2004. [39] J. W. Stayman and J. A. Fessler, “Regularization for uniform spatial resolution properties in penalized-likelihood image reconstruction,” IEEE Trans. Med. Imag., vol. 19, no. 6, pp. 601–615, June 2000. [40] S. Ahn and R. M. Leahy, “Spatial resolution properties of nonquadratically regularized image reconstruction for PET,” in Proc. IEEE Intl. Symp. Biomed. Imag., 2006, pp. 287–290. [41] S. Ahn and J. A. Fessler, “Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms,” IEEE Trans. Med. Imag., vol. 22, no. 5, pp. 613–626, May 2003. [42] C. R. Rao, Linear Statistical Inference and Its Applications. New York: Wiley, 1973. [43] P. Boxandall and H. Liebeck, Vector Calculus. Oxford, U.K.: Clarendon, 1986. [44] P. J. Huber, Robust Statistics. New York: Wiley, 1981. [45] Q. Li, E. Asma, J. Qi, J. R. Bading, and R. M. Leahy, “Accurate estimation of the Fisher information matrix for the PET image reconstruction problem,” IEEE Trans. Med. Imag., vol. 23, no. 9, pp. 1057–1064, Sep. 2004. [46] J. A. Fessler, ASPIRE 3.0 user’s guide: A sparse iterative reconstruction library Commun. Sign. Process. Lab., Dept. EECS, Univ. Michigan, Ann Arbor, Tech. Rep. 293, Jul. 1995. [47] H. Erdo˘gan and J. A. Fessler, “Monotonic algorithms for transmission tomography,” IEEE Trans. Med. Imag., vol. 18, no. 9, pp. 801–814, Sep. 1999. [48] E. J. Hoffman, S. C. Huang, and M. E. Phelps, “Quantitation in positron emission computed tomography: 1 Effect of object size,” J. Comput. Assist. Tomogr., vol. 3, no. 3, pp. 299–308, 1979. [49] J. A. Fessler, Image reconstruction toolbox [Online]. Available: http:// www.eecs.umich.edu/~fessler [50] J. Qi, “Comparison of lesion detection and quantification in MAP reconstruction with Gaussian and non-Gaussian priors,” Int. J. Biomed. Imag., 2006 [Online]. Available: http://www.hindawi.com/GetArticle. aspx?doi=10.1155/IJBI/2006/87567
