Sufficient Statistics as a Generalization of Binning ... - Semantic Scholar
84
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
Sufficient Statistics as a Generalization of Binning in Spectral X-ray Imaging Adam S. Wang*, Student Member, IEEE, and Norbert J. Pelc
Abstract—It is well known that the energy dependence of X-ray attenuation can be used to characterize materials. Yet, even with energy discriminating photon counting X-ray detectors, it is still unclear how to best form energy dependent measurements for spectral imaging. Common ideas include binning photon counts based on their energies and detectors with both photon counting and energy integrating electronics. These approaches can be generalized to energy weighted measurements, which we prove can form a sufficient statistic for spectral X-ray imaging if the weights used, which we term -weights, are basis attenuation functions that can also be used for material decomposition. To study the performance of these different methods, we evaluate the Cramér-Rao lower bound (CRLB) of material estimates in the presence of quantum noise. We found that the choice of binning and weighting schemes can greatly affect the performance of material decomposition. Even with optimized thresholds, binning condenses information but incurs penalties to decomposition precision and is not robust to changes in the source spectrum or object size, although this can be mitigated by adding more bins or removing photons of certain energies from the spectrum. On the other hand, because -weighted measurements form a sufficient statistic for spectral imaging, the CRLB of the material decomposition estimates is identical to the quantum noise limited performance of a system with complete energy information of all photons. Finally, we show that -weights lead to increased conspicuity over other methods in a simulated calcium contrast experiment. Index Terms—Cramér-Rao lower bound, dual energy, energy weighting, photon counting, spectral imaging, sufficient statistic.
I. INTRODUCTION T is well known that a material’s X-ray attenuation is dependent on the energies of the X-rays and that this dependence can be used for material selective imaging [1], [2]. Therefore, there is information about a measured object contained not only in the total number of photons transmitted through the object but also in the energy of each of these photons. How this information is collected greatly affects the ability of a system to efficiently utilize the energy dependence in what is known as spectral X-ray imaging. Early methods for spectral imaging include the use of dual kVp techniques that
I
Manuscript received May 23, 2010; revised July 15, 2010; accepted July 17, 2010. Date of publication August 03, 2010; date of current version December 30, 2010. This work was supported in part by GE Healthcare and in part by the Lucas Foundation. Asterisk indicates corresponding author. *A. S. Wang is with the Departments of Electrical Engineering and Radiology, Stanford University, Stanford, CA 94305 USA (e-mail: [email protected]). N. J. Pelc is with the Departments of Radiology, Bioengineering, and Electrical Engineering, Stanford University, Stanford, CA 94305 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TMI.2010.2061862
require two different exposures—one at a lower energy than the other—or the use of dual layer detectors. Both methods inherently suffer from spectral overlap, which degrades the spectral imaging performance [3]–[7]. Energy discriminating detectors offer the possibility of eliminating spectral overlap by directly partitioning the transmitted spectrum into separate measurements [8]–[10]. Our project began with work to understand and optimize the partition of the spectrum into bins and led to the observation that weighted measurements, a generalization of binning, can achieve universally optimal performance. Here, we evaluate several methods for forming energy dependent measurements in the context of ideal photon counting detectors with energy discriminating capabilities. The work led us to the discovery of a simple and elegant sufficient statistic for spectral X-ray imaging. II. BACKGROUND For now, we assume that the materials being measured have no K-edges within the detected X-ray spectrum; this will be relaxed later. Spectral imaging in the diagnostic energy range is often referred to as dual-energy imaging, which is based on the observation that X-rays in this energy range primarily interact by two physical mechanisms—photoelectric absorption and Compton scattering—and that in the absence of K-edges the energy dependence of these mechanisms is independent of the material [1]. Thus, materials without K-edges within the detected spectrum behave as if they were some linear combination of material independent basis functions (e.g., photoelectric absorption as a function of energy and Compton scattering as a function of energy). In other words, in the energy range of interest, material 1 will have a linear attenuation coefficient at energy that can be expressed as (1) and where and depend on the material, and are the attenuation functions attributable to the photoelectric and Compton interactions, respectively. Moreover, because all materials without K-edges attenuate X-rays based on the same underlying principles, any material can be expressed as a linear combination of any other two materials. For materials 1, 2, and 3, constants and can be found such that (2) Thus, any object can be described as an amount of material 1 and an amount of material 2. To effectively estimate
0278-0062/$26.00 © 2010 IEEE
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
, the object’s attenuation can be measured with X-rays of different energies. An ideal system based on photon counting detectors would keep track of the energy of each detected photon. Data handling requirements generally preclude this, and instead photons are grouped into energy intervals. For example, a system may measure the number of low and high energy photons, where a preset threshold energy determines whether photons are counted in a low or high energy bin. When photons are recorded in this manner based on thresholds, we say they undergo binning. with maximum energy , For an incident spectrum the expected number of transmitted photons at energy through an object of thickness of material 1 and thickness of material 2 is (3) and where we have simplified the notation of to , and , respectively. However, the observed spectrum will be a random vector , is the observed number of photons of energy at where the detector and is a Poisson random variable with mean that is independent of the observed number of photons at other energies [11]. For notational convenience, we will consider all energies to the nearest integer keV. Performance bounds on basis material decomposition for any system configuration can be derived from the Cramér-Rao lower bound (CRLB) [1], [12]. The CRLB provides a lower bound on the covariance of unbiased estimators of deterministic parameters (e.g., ) and can be computed from the distribution of the measurements used to estimate the parameter. Due to the inherent quantum noise of photon statistics, an unbiased estimate of the amount of each basis material will be noisy and have some covariance matrix . For any unbiased estimate with corresponding covariance matrix , the CRLB provides a lower bound so that in the matrix inequality is possense—understood to mean that the quantity and itive semidefinite. In particular, . Since the lower bound is dependent on the incident spectrum and detector configuration, the CRLB can be used as a metric for jointly optimizing the incident spectrum and bin thresholds to minimize the estimation noise. In our work, we consider a fixed incident spectrum and show that the choice of bin thresholds and, later, weighting functions greatly affects dual energy performance. Furthermore, we use the CRLB to compare different measurement schemes to show that different methods have inherently disparate lower bounds on their material decomposition precision. III. METHODS Using Spektr [13], a Matlab implementation of the X-ray spectral model TASMIP [14], we simulated an X-ray tube with constant tube voltage of 120 kV, 2.5 mm Al total filtration, and 0.2 mAs exposure illuminating an ideal 1 mm photon counting energy discriminating detector at 1 m from the source [Fig. 1(a)]. We take basis materials 1 and 2 to be calcium (density: 1.55 g/cm ) and water (density: 1.00 g/cm ), respectively
85
Fig. 1. Simulated 120 kV spectrum with 2.5 mm Al filtration from Spektr (left) and calcium and water linear attenuation functions used as basis functions for dual energy decomposition (right). (a) Spectrum. (b) Basis attenuation.
[Fig. 1(b)] [15]. We chose these as our basis materials because such a decomposition can have direct applications like bone densitometry. Furthermore, we can easily convert our calcium and water decomposition into any other two basis materials via a linear transformation [2]. IV. BINNING bins with energy thresholds keV, the measured number of counts will be , where for each bin is the observed number of photons of energy and is a random vector of the binned counts. In actual is used to prevent electronic noise from systems, threshold keV to be forming counts, and in this work we assume above the noise floor. Availability of the full detected spectrum is a special case of binning in which we have bins each of width 1 keV and is the most information we could ever acquire from the detector. To compare different binning schemes, we turn to the CRLB, whose formulation for binning is shown in Appendix A. The (A.4) represents the absoCRLB for the full spectrum lute minimum covariance performance of any unbiased estimate of the material decomposition for any information collection method, including binning. The only limitation to our estimation precision is the inherent quantum noise. We compare the CRLB (A.3) of any binning scheme to to illustrate the suboptimality in estimating when binning the detected spectrum. If we normalize by , we get a penalty factor for binning. For instance, if calcium is material 1, represents how much higher the best case , will be variance of the estimated amount of calcium, as a result of binning the full detected spectrum since corresponds to the minimum variance of estimate when binning. Clearly, the performance of binned measurements can be no better than that of the full detected spectrum, which encom. passes all information available at the detector, so itself is a function of the binning threshold The CRLB energies, so naturally we would like to find the threshold(s) that minimize some objective function of the CRLB for any given object and incident spectrum. For this discussion, we only consider threshold energies at integer keV values to keep the set of possible threshold energies finite, and we take as our objective , to minimize the variance of the calcium estimate equivalent to minimizing the penalty factor . Note that the For
86
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
Fig. 3. Optimal thresholds for normal binning and for binning with notch filters . The unoptimized three bin and the resulting penalties for object case is included for comparison with the two bins with a notch filter.
t = (0:5; 20)
Q
Fig. 2. Calcium estimate penalty factor as a function of the threshold energy between two bins for three different objects.
resulting optimal threshold(s) will depend on the object’s material decomposition . A. Two Bins Fig. 2 highlights the shortcomings of binning with a single fixed threshold. It shows the penalty factor for calcium estimation as a function of binning threshold for measurements with the incident spectrum of Fig. 1 and three objects characterized by different amounts of calcium and water. The penalty factor varies drastically as a function of threshold energy and across objects. If the threshold energy is chosen to be 55 keV to , an object with a mateminimize the penalty for rial decomposition of 0.5 cm Ca and 20 cm water, the penalty . This factor is higher than the minimum for object suggests that using two bins is far from ideal (because even the is at least 50% for these objects) and minimum penalty far from robust (because using the wrong threshold incurs an even higher penalty). B. Improving Performance The binning penalty and robustness can be improved by increasing the number of bins, but another interesting possibility is removing photons from the spectrum—either at the source or the detector—to increase spectral separation of the two bins. Here we present an idealized scenario where photons in an energy interval are removed without impacting the remainder of the spectrum. We examine the effect of removing photons from the spectrum on the CRLB and the noise penalty. Removing photons from the source spectrum also reduces radiation dose but that effect is not treated here. If we allow for idealized notch filters that remove all photons within a certain range of energies, these gaps in the spectrum should not fall within a bin—rather, they should straddle the binning thresholds since the goal is to increase spectral separation of bins. Then each bin will effectively have lower and higher energy thresholds ( and ), where keV. The measured counts of bin will now . Note that this scenario is a generalization of be the simple binning case, which can be recovered by removing and . all notch filters so that
To optimize our thresholds in regular binning, we perform an exhaustive search over the set of integer such that keV, finding the set of that minimizes , thus enabling the most precise estimate of given bins and incident spectrum . If notch filters are allowed, the search is over the space of all effective thresholds such . If the exhaustive that search is too taxing, e.g., when the number of bins is large, then branching algorithms can be used. To calculate the noise penalty , comparison was made to the CRLB with the full spectrum (i.e., without photons removed) in all cases. Fig. 3 shows the results for four binning strategies: and 3, with and without all energies, for a fixed object . The effective energy ranges of the bins are shown as bars and notch filters are shown as gaps. The column to the right . We find that removing lists the calcium estimate penalty photons in the middle of the spectrum from either of the two , bins, with the notch filter optimized for the measurement of significantly below that of the normal two bin reduced case. Of course we would expect that retaining these photons in a separate bin would be even better. While this indeed provides a further improvement, the additional benefit is marginal because this is an unoptimized three bin design. The optimized three bin design performs even better and has thresholds that are different than those of the case of two bins with a notch filter. As the number of bins is increased, the penalty factor goes to unity. These observations apply to other objects as well, although each object may have different optimal thresholds. In some hardware implementations (e.g., using comparators) inserting a gap between two bins is essentially as complex as adding an additional bin. In those cases, retaining the photons in an (optimized) additional bin is preferred. However, an important point in the results presented above is that much of the benefit of additional bins can come from providing separation in the spectra contained in each bin. If the notch filter is implemented at the source, the benefit to material decomposition also comes with a benefit in dose reduction without additional complexity at the detector. This strategy of applying notch filters is beneficial when es(and also ) but may not apply for other tasks timating where the total number of photons is more important than spectral separation of the bins, such as measuring the total attenuation or forming monoenergetic images. In those cases, retaining the gap photons in a separate bin is a very significant benefit.
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
87
Fig. 4. Energy weighting functions for three special cases: (a) CIX - PC/EI, (b) binning with all energies, and (c) binning with notch filters.
V. WEIGHTED MEASUREMENTS A further generalization of binning counts is forming weighted measurements. The general concept is that photon counting energy discriminating detectors can be used to form weighted measurements of photon counts, where the weighting are a function of photon energy. The functions , where result is a set of measurements
(4)
The concept of energy weighting or weighting functions is not new [8], [16]–[18], and when at least two measurements with different weighting functions are acquired, dual-energy decomposition can be performed. A special case is a counting and integrating X-ray (CIX) detector that does both photon counting and energy integration, providing measurements with weights and [19]–[21] [Fig. 4(a)]. Another special case is binning—the weighting functions are binary and nonoverlapping [Fig. 4(b)]. When notch filters are considered, the weights are zero for energies filtered out, creating separation between the bins [Fig. 4(c)]. In the most general case, the weighting functions are arbitrary and can take on any real value at every energy so that, unlike binning where photons contribute to no more than one bin measurement, here photons can contribute to all weighted measurements. As with binning, forming weighted measurements reduces the amount of data that has to be stored and processed. Since binning, with or without excluding some energies, are special cases, there must be sets of weighting functions that perform at least as well as binning. We seek to characterize the performance as a function of the choice of weights and ultimately to find the optimal weighting functions. For weighted measurements, the random vector (4) is composed of weighted sums of independent, Poisson distributed to meavalues. Every photon of energy contributes value surement . Except for the special case of binary, nonoverlapping weights, since the weighted measurements are obtained from the same realization of , the elements of vector are no longer independent nor Poisson distributed. In such cases, multivariate Gaussian models are commonly used because of their convenient analytical form, tractability, and ability to accurately capture the first- and second-moments of [12], [22], where the weights determine the mean and covariance (B.1). For a sufficiently high transmitted number of photons, this model works well, but when only a few photons are transmitted through
the object, the weighted sum of discrete Poisson distributions may not be well-approximated by a smooth, continuous distribution. Nonetheless, assuming a jointly Gaussian distribution, the CRLB approach can be used to evaluate estimator performance as a function of incident spectrum and energy weights that should be valid at least at high signal-to-noise ratio (SNR) (Appendix B). For comparison, we also find the exact CRLB . The distribution of a weighted for weighted measurements, sum of Poisson random variables is the convolution of the individual scaled Poisson distributions, but has no known tractable can be form. However, the exact distribution and CRLB computed numerically using empirical characteristic functions (Appendix C). A. Optimal Weights In principle, we can find the optimal weights by brute force, treating the weights as variables over which we minimize some objective function [23], [24]. A more elegant approach, if possible, is to find weights that form a sufficient statistic for the . If it could be fully known detected spectrum shown that a set of weights produces weighted measurements that are sufficient statistics for estimating , then the weighted measurements would contain as much information about material thicknesses as the full spectrum . More formally, let be a random vector whose distribution is a sufficient is parameterized by . Then, a statistic is statistic for if the conditional distribution of given independent of [25]. Intuitively, simple examples of sufficient statistics include counting the number of heads in a sequence of coin flips rather than keeping track of the individual flips to estimate the probability of heads or using the sample mean of a sequence of independent, identically distributed random variables to estimate the mean of a normal or Poisson distribution. To find a sufficient statistic, we use the Factorization Theorem is the joint probability mass [25], which states that if function (pmf) of given , then statistic is sufand ficient for if and only if there exist functions such that, for all and . Under this condition can be used to estimate as well as the original measurements . Thus, for our case, the pmf of the detected spectrum given the material thicknesses needs to be decomposed into the product of one function of the statistic and , and a function of the detected spectrum that does not depend on . We begin by writing the joint pmf of the detected spectrum as the product of the pmf at each energy since they are independent given . We use our knowledge that the photon counts are Poisson
88
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
distributed and assume that the incident spectrum and can be treated as a constant
is known
(5) where
(6)
(7) and (8) Therefore, we have proven that two weighted measurements and as their using the attenuation basis functions weights form a sufficient statistic for dual energy imaging. By definition, this sufficient statistic (i.e., these weighted measurements) should achieve optimal performance for estimating any with any incident spectrum . Our only assumption is that the object’s attenuation can be described by two basis functions. We call these weights -weights. VI. COMPARISON To illustrate that the -weighted measurements form a sufficient statistic, we can compute the CRLB. We expect the performance of the -weighted measurements to be identical to that for all objects. In of the full detected spectrum addition, because the -weighted measurements only require two values derived from the detector signals, for comparison we compute the CRLB for two other schemes that report two measurements: optimal binning with two bins and using a CIX detector that does both photon counting and energy integration. A. Performance Curves cm, and Consider a range of relevant ( cm). For the 120 kVp spectrum [Fig. 1(a)], the minvariance ( , in units of cm ) as a funcimum achievable tion of the amount of water and calcium is shown in Fig. 5(a) as
log contour lines. The uncertainty increases with increasing ob(the ject attenuation. Fig. 5(b) shows the penalty function ) when two bins are used with a increase in the variance of fixed threshold of 57 keV, the threshold that minimizes the maximum (over ) penalty factor, and Fig. 5(c) shows similar data for the CIX detector. Binning with a threshold of 57 keV incurs a penalty by a factor of 3.4 at the lower left corner and by at least a factor of 1.7 throughout. A CIX detector incurs lower penalties, especially at larger but still incurs a penalty of . at least 22% and up to a factor of two in Var to for -weights when using the The ratio of Gaussian model is shown in Fig. 5(d) for . This ratio is near 1 at low object thickness, as expected, but drops below 1 for thicker objects where the number of photons is low. However, we know that any unbiased estimator cannot truly have better performance than an estimator having knowledge of the full detected spectrum. Therefore, the multivariate Gaussian model may not be accurate for weighted measurements when the number of photons detected is low. For example, the expected number of photons detected for the thickest object is only 39. Poisson distributions with low mean counts can only take on a set of discrete values and are positively skewed, unlike Gaussian distributions. This conclusion is supported by the fact that the exact CRLB of the -weighted measurements from the empirical characteristic function (Appendix C) is well within 0.1% of the CRLB of the full detected spectrum [Fig. 5(e)]. The accuracy of this method simply depends on the discretization of the Fourier space variable and becomes a numerical issue. These results show that -weighted measurements achieve as knowledge of the the same performance in estimating full spectrum. Although not shown here, the and for -weighted measurements also match those of the full spectrum, while those of binning and the CIX exhibit counterparts under the similar penalty ratios to their different measurement schemes. It is possible to have objects that decompose into a negative amount of calcium or water (but not both). While only quadrant of the calcium/water decomposition space I is shown in Fig. 5, our results extend smoothly into the valid regions of quadrants II and IV, and none of our methods preclude negative amounts of either basis material. B. Wedge Phantom We also simulated the dual-energy performance of the three different binning or weighting schemes on a phantom designed to compare the calcium detectability of the methods. Consider a projection image with 0.2 mAs exposure of a water wedge phantom varying from 0 to 40 cm thickness in the horizontal direction with square calcium contrast elements ranging in area from 1 mm to 1 cm overlaid on top [Fig. 6(a)]. These elements have a thickness such that the predicted (full area) SNR of each calcium element from the full detected spectrum is 4 (at the threshold of detectability). Because the necessary amount of calcium to maintain a constant SNR increases as the water gets thicker and as the elements becomes smaller, the calcium images are displayed on a spatially varying grayscale window in the horizontal direction, where is specified as a function of horizontal position [Fig. 6(b) and (c)].
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
89
Q
Fig. 5. CRLB performance over a wide range of material thicknesses for (a) the full detected spectrum (cm ) and the multiplicative penalty factor for (b) two bins with a 57 keV threshold (no notch), (c) counting and integrating, (d) optimal -weights with a Gaussian model, and (e) optimal -weights with the exact = 57 keV. (c) numerically computed CRLB. Note the different contour levels in the subfigures. (a) Var(^ ): Full spectrum. (b) Var(^ ) penalty: 2 bins, Var(^ ) penalty: CIX. (d) Var(^ ) penalty: -weights, Gaussian. (e) Var(^ ) penalty: -weights, Exact.
t
t
t
The numerical phantom was designed so that all contrast elements project onto an integer number of detectors and are perfectly aligned with the detector grid so as to avoid any partial volume effects. Furthermore, while it is a practical concern, we leave aside the issue of spatial resolution, count rates, and energy resolution by assuming an idealized photon counting detector with no cross-talk, no count rate losses, and perfect energy resolution. Transmission of the X-ray spectrum and detection with energy resolution of 1 keV and perfect detection efficiency was simulated. Quantum noise was incorporated by independently sampling the appropriate Poisson distribution at every 1 keV energy. The simulated noisy detected spectra were then subjected to each of the measurement schemes so that all binned and weighted measurements were obtained from the same realization of detected spectra. Maximum-likelihood estimation (MLE) is a commonly used method for estimating given noisy measurements [9], [26], [27] and was a convenient choice given the extensive use of the log-likelihood functions (A.2) and (B.1) in deriving the CRLB. We used Matlab’s Optimization Toolbox (v3.1.1) to find the maximum-likelihood solution
(9) Although more sophisticated methods may exist for estimating , these are beyond the focus of this paper. For the CIX and weighted measurements, we used the Gaussian model. We compare the ML decomposition of measurements taken with: 1) two
t
t
bins with a 57 keV threshold (no notch); 2) a CIX detector (i.e., ); 3) -weights. The resulting ML decompositions for the calcium component are shown in Fig. 6. Detection of the targets is challenging in all the images since the phantom was designed to make this so, but the -weights image is noticeably superior. The images illustrate that using -weights to form two measurements allows for increased conspicuity of all contrast elements as compared with a two bin or CIX approach. The estimates are shown with . As the spatially varying grayscale window expected from the CRLB comparisons (Fig. 5), the noise in the calcium image is highest when using two bins and lowest when using -weights. The empirical SNR as measured with 10 000 realizations of this experiment is consistent with the theoretical SNR using the CRLB (Fig. 5), although a bias in the estimates becomes evident at the thick end of the wedge. As predicted by Fig. 5(c), the advantage of -weights over CIX is largest in the lower left corner of the wedge and drops off to the upper right. VII. DISCUSSION AND CONCLUSION A simple explanation for the improvement of binning performance from including a notch between bins is that the sensitivity of the ML estimate to detected photons depends on the photon energy. Fig. 7 plots the expected transmitted spectrum and the change in the material decomposition estimates using the full detected spectrum when one additional photon of energy is detected beyond the expected number of counts for the object . Photons near 64 keV provide very little information to the calcium estimate, whereas additional photons at lower energies would induce the MLE to predict less
90
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
0
Fig. 6. Calcium thickness estimation task for overhead projection with different weighting schemes, displayed on a spatially varying grayscale range [ (x) + x)], where (a) is an illustration of the water wedge with calcium contrast elements of different sizes and thicknesses overlaid, (b) is the true calcium thickness, (c) shows the grayscale window as a function of horizontal direction x, and the estimated calcium thickness is from (d) two bins with a 57 keV threshold (no notch), (e) counting and integrating, and (f) optimal -weights. (a) Wedge phantom. (b) True t . (c) Spatially varying . (d) t^ : 2 bins, = 57 keV. (e) t^ : CIX. (f) t^ : -weights.
(
Fig. 7. The expected transmitted spectrum for object t = (0:5; 20) and the sensitivity of estimates t^ and t^ to the number of photons at energy j when using the full detected spectrum. The shaded region is the range of photon energies that should be removed when estimating t with two bins.
calcium and more water and vice-versa for additional photons at higher energies because the ratio of calcium to water attenuation is much higher at lower energies. While photons in the energy notch contain some information, once binned, photons of different energies within bin lose their individual energy identities and cannot be processed any differently than other photons in the same bin that have much more information. They therefore weaken the overall material discrimination ability. Material decomposition is not only insensitive to photons in the shaded energy range (Fig. 7), but these photons also move the average energy of each bin closer to each other and reduce the spectral separation between the bins. It has been shown for
dual kVp and dual layer methods that spectral separation is desired [3]–[7]. Introducing a notch increases the spectral separation. While this is at the expense of fewer photons per bin, it nonetheless improves the precision of the calcium (and water) estimate. In analogous prior work on dual energy imaging with dual layer detectors, it was found that adding a metal filter between the front and back layers improves material separation performance. In that case and as in our work, discarding photons of intermediate energy improves the precision of material decomposition even though it decreases the number of photons that contribute to the measurements [5]–[7]. Ideally, X-rays of the discarded energies would be removed from the incident spectrum to reduce radiation dose, and the net effect is equivalent to binning with a notch between bins. If a notch filter is realizable (either at the source or detector), material decomposition precision can be improved with the same data storage, transmission, and estimation requirements as without one. In practice such selectivity may not be possible to achieve, but the large fraction of omitted photons (Fig. 7) highlights the important role that source filtration and optimization can play in dose efficiency. It is important to note that the implementation of two bins with a notch at the detector using comparators requires the same number of comparators as using three bins. Without a doubt, forming a third bin from the discarded photons instead and utilizing it enables better material decomposition estimation. Although adding more bins and notch filter(s) improves the calcium estimation task performance, the optimal thresholds are a function of object size and the incident spectrum. No fixed binning scheme is universally optimal. Therefore, binning cannot
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
Fig. 8. Two sets of equivalent weights where the weights on the left are simply the linear attenuation functions of calcium, water and iodine and the weights on the right are linearly independent combinations of the attenuation functions. (a) Weights ; , and (b) Linear combinations.
form a sufficient statistic for the full detected spectrum, although for a large number of bins, the difference in performance may be negligible. Using -weights, the detector performance is optimized for all spectra and objects, and thus we have decoupled the joint optimization of the spectrum and detector weights into simply an optimization of the spectrum. The -weighted measurements could perhaps be formed at the detector using analog processing with two parallel circuits and [i.e., an energy dependent that map pulse energy to gain of and ] followed by integrators that and , respectively, or digitally by determining represent the energy of a detected photon from its pulse height and adding and to two different accumulators that represent and , respectively. The values of and can be embedded in hardware or stored in a look-up table. Incorporating the -weights directly into the detectors may require significant hardware redesign. Alternatively, the -weights can be applied to binned counts (assuming there are multiple bins) to reduce the dimensionality of the measurements to only two, although this effectively imposes weights that are a piecewise constant approximation to the attenuation curves and would have to be further investigated. Because any linearly independent comand , the attenuation due to the photoelectric binations of effect and Compton scattering, respectively, can be used as and , there are infinitely many selections for and available. Should the object possibly contain an additional material that deviates from the dual basis material with attenuation assumption, such as a material with a K-edge within the detected spectrum, the attenuation of this material can be used as an additional basis function. The sufficient statistic proof (5) can be extended to show that only three weighted measureand . To illustrate the ments are needed, with weights previous two points, the selected weights could be and [Fig. 8(a), where material 3 is iodine, density: 4.93 g/cm ], or they could be a more complicated linear combination, such as normalized and flattened weights, where the third set of weights has values near 0 for energies greater than 60 keV [Fig. 8(b)]. The performance would be identical since any set of three weighted measurements is entirely recoverable from any other set of three weighted measurements through a simple linear transformation, even in the presence of quantum noise. Using the fact that -weighted measurements form a sufficient statistic, we have shown that a multivariate Gaussian
91
model for the distribution of weighted measurements may not be accurate at low photon counts (the CRLB is lower than possible). This suggests that using a Gaussian model for the estimation of the thicknesses can yield biased estimates. The loss of accuracy in a Gaussian model is not limited to just -weights and extends to CIX detectors and binning because the underlying problem is that the measurements are a weighted sum of Poisson random variables instead of the commonly assumed Gaussian distribution. Although bias correction methods and efficient estimators are beyond the scope of this paper, we expect them to be important issues in accurate and SNR-efficient dual energy decomposition. One solution may be to select a more appropriate likelihood function, such as a Gaussian approximation when the number of photons is large and switching to a discrete distribution at lower counts using the empirical characteristic function. Future work will also focus on other practical considerations of forming sufficient statistics in photon counting detectors, such as incorporating spectral response functions [9], [12] and forming weighted sums of binned counts [8], [10], [17], [18]. In conclusion, the choice of binning and weighting schemes can greatly affect the performance of dual-energy material decomposition. Even with optimal thresholds, binning incurs penalties to decomposition precision and is not robust to changes in the source spectrum or object size, although this can be mitigated by adding more bins. On the other hand, two -weighted measurements in theory form a sufficient statistic for recording the detected photons and in our simulated experiment do well when implemented with a multivariate Gaussian approximation to the true probability mass function. The -weighted measurements outperform binned or CIX measurements in our contrast wedge phantom simulation. Estimating two unknown parameters is at the core of dual-energy theory, and we have proven that only two -weighted measurements are sufficient to accomplish this without any loss of information.
APPENDIX A BINNING CRLB To find the CRLB, we compute the Fisher information matrix, whose elements are given by
(A.1)
is the log-likelihood function of acquiring meawhere surements for object and the expected value of the argument is over , a random variable. Then, if is an unbiased estimator of , the CRLB provides a lower bound on the covariance matrix of by stating that , where . For binned measurements, each bin is the sum of independent Poisson random variables. Therefore, each bin count is itself a Poisson random variable, parameterized by for the case when all energies are used and when
92
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
notch filters are used. For measurements of an object characterized by thickness vector , the log-likelihood function (ignoring the constant term) is (A.2) As shown in [12], the resulting CRLB for binning,
The characteristic function of a distribution is the Fourier transform of its probability density function. For a random variable having a Poisson distribution with mean , the characteristic function is
, is
(C.1) (C.2)
(A.3) . One interpretation
where
APPENDIX C EMPIRICAL CHARACTERISTIC FUNCTIONS
is that each expected measurement is weighted by of that is representative of over energies to . In the some special case of binning when we use the full detected spectrum , the CRLB for the full spectrum is
. The characteristic function has where in this section the property that if the random variable is multiplied by a constant , then random variable will have a characteristic . Furthermore, the characteristic function function of a sum of independent random variables is equivalent to the product of their characteristic functions [29]. That is, random variable
(C.3) In the case of two weighted measurements, we can also consider the characteristic function of the multivariate random variable . For
(C.4)
(A.4) where
has a characteristic function
Note that because the characteristic function of a random variable is the Fourier transform of its distribution
.
(C.5)
APPENDIX B GAUSSIAN CRLB taken with weighting For measurements functions , the expected value of mea. In addition, , the covarisurements will be . ance matrix of measurements will be Ignoring the constant term, the log-likelihood function for a Gaussian model of weighted measurements is
This can be numerically computed by discretizing and evalat these points, giving what is known as the empiruating ical characteristic function [30], [31]. Then, taking the 2-D inverse discrete Fourier transform of the empirical characteristic . The resofunction gives us the empirical distribution of lution of depends on the sampled range of , while the range of depends on the spacing of the samples of [32]. Expanding this method, we can also compute the Fisher matrix exactly
(B.1) The elements of the Fisher information matrix are derived in [12], [28] and given by (C.6)
(B.2) Therefore, the CRLB for weighted measurements when modeled by a Gaussian distribution is .
where
is the support of . Then the exact CRLB comes from . We have shown that can be found using (C.5), and the other terms in the integral can be numerically computed as well. These terms all contain partial derivative(s), but because the variables and are the Fourier transform duals
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
of one another, the partial derivative with respect to can move inside the Fourier transform. For example,
(C.7) The other partial derivatives can be similarly computed numerically for any . Then, because a discretized form of each term in the integrand of (C.6) can be found, the Riemann integral can be found by summing the expression over all samples . to arrive at
REFERENCES [1] R. E. Alvarez and A. Macovski, “Energy-selective reconstructions in X-ray computerized tomography,” Phys. Med. Biol., vol. 21, no. 5, pp. 733–744, 1976. [2] L. A. Lehmann, R. E. Alvarez, A. Macovski, W. R. Brody, N. J. Pelc, S. J. Riederer, and A. L. Hall, “Generalized image combinations in dual KVP digital radiography,” Med. Phys., vol. 8, no. 5, pp. 659–667, 1981. [3] F. Kelcz, P. M. Joseph, and S. K. Hilal, “Noise considerations in dual energy CT scanning,” Med. Phys., vol. 6, no. 5, pp. 418–425, 1979. [4] A. N. Primak, J. C. Ramirez Giraldo, X. Liu, L. Yu, and C. H. McCollough, “Improved dual-energy material discrimination for dual-source CT by means of additional spectral filtration,” Med. Phys., vol. 36, no. 4, pp. 1359–1369, 2009. [5] G. T. Barnes, R. A. Sones, M. M. Tesic, D. R. Morgan, and J. N. Sanders, “Detector for dual-energy digital radiography,” Radiology, vol. 156, no. 2, pp. 537–540, 1985. [6] H. N. Cardinal and A. Fenster, “Theoretical optimization of a split septaless xenon ionization detector for dual-energy chest radiography,” Med. Phys., vol. 15, no. 2, pp. 167–180, 1988. [7] G. M. Stevens and N. J. Pelc, “Depth-segmented detector for X-ray absorptiometry,” Med. Phys., vol. 27, no. 5, pp. 1174–1184, 2000. [8] J. Karg, D. Niederlöhner, J. Giersch, and G. Anton, “Using the Medipix2 detector for energy weighting,” Nucl. Instrum. Meth. A, vol. 546, pp. 306–311, 2005. [9] J. P. Schlomka, E. Roessl, R. Dorscheid, S. Dill, G. Martens, T. Istel, C. Bäumer, C. Herrmann, R. Steadman, G. Zeitler, A. Livne, and R. Proksa, “Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography,” Phys. Med. Biol., vol. 53, no. 15, pp. 4031–4047, 2008. [10] P. M. Shikhaliev, “Tilted angle CZT detector for photon counting/energy weighting X-ray and CT imaging,” Phys. Med. Biol., vol. 51, no. 17, pp. 4267–4287, 2006. [11] A. Macovski, Medical Imaging Systems. Englewood Cliffs, NJ: Prentice-Hall, 1983. [12] E. Roessl and C. Herrmann, “Cramér-Rao lower bound of basis image noise in multiple-energy X-ray imaging,” Phys. Med. Biol., vol. 54, no. 5, pp. 1307–1318, 2009.
93
[13] J. H. Siewerdsen, A. M. Waese, D. J. Moseley, S. Richard, and D. A. Jaffray, “Spektr: A computation tool for X-ray spectral analysis and imaging system optimization,” Med. Phys., vol. 31, no. 11, pp. 3057–3067, Nov. 2004. [14] J. M. Boone and J. A. Seibert, “An accurate method for computer-generating tungsten anode X-ray spectra from 30 to 140 kV,” Med. Phys., vol. 24, no. 11, pp. 1661–1670, Nov. 1997. [15] M. J. Berger, J. H. Hubbell, S. M. Seltzer, J. Chang, J. S. Coursey, R. Sukumar, and D. S. Zucker, XCOM: Photon Cross Section Database (Version 1.3) National Institute of Standards and Technology. Gaithersburg, MD [Online]. Available: http://physics.nist.gov/xcom, Sep. 2009 [16] M. J. Tapiovaara and R. F. Wagner, “SNR and DQE analysis of broad spectrum X-ray imaging,” Phys. Med. Biol., vol. 30, no. 6, pp. 519–529, 1985. [17] R. N. Cahn, B. Cederstrom, M. Danielsson, A. Hall, M. Lundqvist, and D. Nygren, “Detective quantum efficiency dependence on X-ray energy weighting in mammography,” Med. Phys., vol. 26, no. 12, pp. 2680–2683, 1999. [18] J. Giersch, D. Niederlöhner, and G. Anton, “The influence of energy weighting on X-ray imaging quality,” Nucl. Instrum. Meth. A, vol. 531, pp. 68–74, 2004. [19] E. Roessl, A. Ziegler, and R. Proksa, “On the influence of noise correlations in measurement data on basis image noise in dual-energy like X-ray imaging,” Med. Phys., vol. 34, no. 3, pp. 959–966, 2007. [20] E. Kraft, P. Fischer, M. Karagounis, M. Koch, H. Kruger, I. Peric, N. Wermes, C. Herrmann, A. Nascetti, M. Overdick, and W. Ruetten, “Counting and integrating readout for direct conversion X-ray imaging: Concept, realization and first prototype measurements,” IEEE Trans. Nucl. Sci., vol. 54, no. 2, pp. 383–390, Apr. 2007. [21] H. Krüger, J. Fink, E. Kraft, N. Wermes, P. Fischer, I. Peric, C. Herrman, M. Overdick, and W. Rütten, “CIX—A detector for spectral enhanced X-ray imaging in simultaneous counting and integrating,” in Proc. SPIE, 2008, vol. 6913, p. 69130P. [22] M. Firsching, A. P. Butler, N. Scott, N. G. Anderson, T. Michel, and G. Anton, “Contrast agent recognition in small animal CT using the Medipix2 detector,” Nucl. Instrum. Meth. A, vol. 607, pp. 179–182, 2009. [23] A. S. Wang and N. J. Pelc, “Optimal energy thresholds and weights for separating materials using photon counting X-ray detectors with energy discriminating capabilities,” in Proc. SPIE, 2009, vol. 7258, p. 725872. [24] D. Niederlöhner, J. Karg, J. Giersch, and G. Anton, “The energy weighting technique: Measurements and simulations,” Nucl. Instrum. Meth. A, vol. 546, pp. 37–41, 2005. [25] G. Casella and R. L. Berger, “Principles of data reduction,” in Statistical Inference, 2nd ed. Belmont, CA: Duxbury Press, 2001. [26] J. A. Fessler, I. Elbakri, P. Sukovic, and N. H. Clinthorne, “Maximum-likelihood dual-energy tomographic reconstruction,” in Proc. SPIE, 2002, vol. 4684, p. 468405. [27] J. Xu, E. C. Frey, K. Taguchi, and B. M. W. Tsui, “A poisson likelihood iterative reconstruction algorithm for material decomposition in CT,” in Proc. SPIE, 2007, vol. 6510, p. 65101Z. [28] S. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Upper Saddle River, NJ: Prentice Hall, 1993. [29] A. Papoulis and S. U. Pillai, “Functions of one random variable,” in Probability, Random Variables and Stochastic Processes, 4th ed. Boston, MA: McGraw-Hill, 2001. [30] A. Feuerverger and R. A. Mureika, “The empirical characteristic function and its applications,” Ann. Stat., vol. 5, no. 1, pp. 88–97, 1977. [31] A. Feuerverger and P. McDunnough, “On some Fourier methods for inference,” J. Am. Stat. Assoc., vol. 76, no. 374, pp. 379–387, 1981. [32] A. V. Oppenheim, R. W. Schafer, and J. R. Buck, Discrete-Time Signal Processing, 2nd ed. Upper Saddle River, NJ: Prentice Hall, 1999.
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
Sufficient Statistics as a Generalization of Binning in Spectral X-ray Imaging Adam S. Wang*, Student Member, IEEE, and Norbert J. Pelc
Abstract—It is well known that the energy dependence of X-ray attenuation can be used to characterize materials. Yet, even with energy discriminating photon counting X-ray detectors, it is still unclear how to best form energy dependent measurements for spectral imaging. Common ideas include binning photon counts based on their energies and detectors with both photon counting and energy integrating electronics. These approaches can be generalized to energy weighted measurements, which we prove can form a sufficient statistic for spectral X-ray imaging if the weights used, which we term -weights, are basis attenuation functions that can also be used for material decomposition. To study the performance of these different methods, we evaluate the Cramér-Rao lower bound (CRLB) of material estimates in the presence of quantum noise. We found that the choice of binning and weighting schemes can greatly affect the performance of material decomposition. Even with optimized thresholds, binning condenses information but incurs penalties to decomposition precision and is not robust to changes in the source spectrum or object size, although this can be mitigated by adding more bins or removing photons of certain energies from the spectrum. On the other hand, because -weighted measurements form a sufficient statistic for spectral imaging, the CRLB of the material decomposition estimates is identical to the quantum noise limited performance of a system with complete energy information of all photons. Finally, we show that -weights lead to increased conspicuity over other methods in a simulated calcium contrast experiment. Index Terms—Cramér-Rao lower bound, dual energy, energy weighting, photon counting, spectral imaging, sufficient statistic.
I. INTRODUCTION T is well known that a material’s X-ray attenuation is dependent on the energies of the X-rays and that this dependence can be used for material selective imaging [1], [2]. Therefore, there is information about a measured object contained not only in the total number of photons transmitted through the object but also in the energy of each of these photons. How this information is collected greatly affects the ability of a system to efficiently utilize the energy dependence in what is known as spectral X-ray imaging. Early methods for spectral imaging include the use of dual kVp techniques that
I
Manuscript received May 23, 2010; revised July 15, 2010; accepted July 17, 2010. Date of publication August 03, 2010; date of current version December 30, 2010. This work was supported in part by GE Healthcare and in part by the Lucas Foundation. Asterisk indicates corresponding author. *A. S. Wang is with the Departments of Electrical Engineering and Radiology, Stanford University, Stanford, CA 94305 USA (e-mail: [email protected]). N. J. Pelc is with the Departments of Radiology, Bioengineering, and Electrical Engineering, Stanford University, Stanford, CA 94305 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TMI.2010.2061862
require two different exposures—one at a lower energy than the other—or the use of dual layer detectors. Both methods inherently suffer from spectral overlap, which degrades the spectral imaging performance [3]–[7]. Energy discriminating detectors offer the possibility of eliminating spectral overlap by directly partitioning the transmitted spectrum into separate measurements [8]–[10]. Our project began with work to understand and optimize the partition of the spectrum into bins and led to the observation that weighted measurements, a generalization of binning, can achieve universally optimal performance. Here, we evaluate several methods for forming energy dependent measurements in the context of ideal photon counting detectors with energy discriminating capabilities. The work led us to the discovery of a simple and elegant sufficient statistic for spectral X-ray imaging. II. BACKGROUND For now, we assume that the materials being measured have no K-edges within the detected X-ray spectrum; this will be relaxed later. Spectral imaging in the diagnostic energy range is often referred to as dual-energy imaging, which is based on the observation that X-rays in this energy range primarily interact by two physical mechanisms—photoelectric absorption and Compton scattering—and that in the absence of K-edges the energy dependence of these mechanisms is independent of the material [1]. Thus, materials without K-edges within the detected spectrum behave as if they were some linear combination of material independent basis functions (e.g., photoelectric absorption as a function of energy and Compton scattering as a function of energy). In other words, in the energy range of interest, material 1 will have a linear attenuation coefficient at energy that can be expressed as (1) and where and depend on the material, and are the attenuation functions attributable to the photoelectric and Compton interactions, respectively. Moreover, because all materials without K-edges attenuate X-rays based on the same underlying principles, any material can be expressed as a linear combination of any other two materials. For materials 1, 2, and 3, constants and can be found such that (2) Thus, any object can be described as an amount of material 1 and an amount of material 2. To effectively estimate
0278-0062/$26.00 © 2010 IEEE
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
, the object’s attenuation can be measured with X-rays of different energies. An ideal system based on photon counting detectors would keep track of the energy of each detected photon. Data handling requirements generally preclude this, and instead photons are grouped into energy intervals. For example, a system may measure the number of low and high energy photons, where a preset threshold energy determines whether photons are counted in a low or high energy bin. When photons are recorded in this manner based on thresholds, we say they undergo binning. with maximum energy , For an incident spectrum the expected number of transmitted photons at energy through an object of thickness of material 1 and thickness of material 2 is (3) and where we have simplified the notation of to , and , respectively. However, the observed spectrum will be a random vector , is the observed number of photons of energy at where the detector and is a Poisson random variable with mean that is independent of the observed number of photons at other energies [11]. For notational convenience, we will consider all energies to the nearest integer keV. Performance bounds on basis material decomposition for any system configuration can be derived from the Cramér-Rao lower bound (CRLB) [1], [12]. The CRLB provides a lower bound on the covariance of unbiased estimators of deterministic parameters (e.g., ) and can be computed from the distribution of the measurements used to estimate the parameter. Due to the inherent quantum noise of photon statistics, an unbiased estimate of the amount of each basis material will be noisy and have some covariance matrix . For any unbiased estimate with corresponding covariance matrix , the CRLB provides a lower bound so that in the matrix inequality is possense—understood to mean that the quantity and itive semidefinite. In particular, . Since the lower bound is dependent on the incident spectrum and detector configuration, the CRLB can be used as a metric for jointly optimizing the incident spectrum and bin thresholds to minimize the estimation noise. In our work, we consider a fixed incident spectrum and show that the choice of bin thresholds and, later, weighting functions greatly affects dual energy performance. Furthermore, we use the CRLB to compare different measurement schemes to show that different methods have inherently disparate lower bounds on their material decomposition precision. III. METHODS Using Spektr [13], a Matlab implementation of the X-ray spectral model TASMIP [14], we simulated an X-ray tube with constant tube voltage of 120 kV, 2.5 mm Al total filtration, and 0.2 mAs exposure illuminating an ideal 1 mm photon counting energy discriminating detector at 1 m from the source [Fig. 1(a)]. We take basis materials 1 and 2 to be calcium (density: 1.55 g/cm ) and water (density: 1.00 g/cm ), respectively
85
Fig. 1. Simulated 120 kV spectrum with 2.5 mm Al filtration from Spektr (left) and calcium and water linear attenuation functions used as basis functions for dual energy decomposition (right). (a) Spectrum. (b) Basis attenuation.
[Fig. 1(b)] [15]. We chose these as our basis materials because such a decomposition can have direct applications like bone densitometry. Furthermore, we can easily convert our calcium and water decomposition into any other two basis materials via a linear transformation [2]. IV. BINNING bins with energy thresholds keV, the measured number of counts will be , where for each bin is the observed number of photons of energy and is a random vector of the binned counts. In actual is used to prevent electronic noise from systems, threshold keV to be forming counts, and in this work we assume above the noise floor. Availability of the full detected spectrum is a special case of binning in which we have bins each of width 1 keV and is the most information we could ever acquire from the detector. To compare different binning schemes, we turn to the CRLB, whose formulation for binning is shown in Appendix A. The (A.4) represents the absoCRLB for the full spectrum lute minimum covariance performance of any unbiased estimate of the material decomposition for any information collection method, including binning. The only limitation to our estimation precision is the inherent quantum noise. We compare the CRLB (A.3) of any binning scheme to to illustrate the suboptimality in estimating when binning the detected spectrum. If we normalize by , we get a penalty factor for binning. For instance, if calcium is material 1, represents how much higher the best case , will be variance of the estimated amount of calcium, as a result of binning the full detected spectrum since corresponds to the minimum variance of estimate when binning. Clearly, the performance of binned measurements can be no better than that of the full detected spectrum, which encom. passes all information available at the detector, so itself is a function of the binning threshold The CRLB energies, so naturally we would like to find the threshold(s) that minimize some objective function of the CRLB for any given object and incident spectrum. For this discussion, we only consider threshold energies at integer keV values to keep the set of possible threshold energies finite, and we take as our objective , to minimize the variance of the calcium estimate equivalent to minimizing the penalty factor . Note that the For
86
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
Fig. 3. Optimal thresholds for normal binning and for binning with notch filters . The unoptimized three bin and the resulting penalties for object case is included for comparison with the two bins with a notch filter.
t = (0:5; 20)
Q
Fig. 2. Calcium estimate penalty factor as a function of the threshold energy between two bins for three different objects.
resulting optimal threshold(s) will depend on the object’s material decomposition . A. Two Bins Fig. 2 highlights the shortcomings of binning with a single fixed threshold. It shows the penalty factor for calcium estimation as a function of binning threshold for measurements with the incident spectrum of Fig. 1 and three objects characterized by different amounts of calcium and water. The penalty factor varies drastically as a function of threshold energy and across objects. If the threshold energy is chosen to be 55 keV to , an object with a mateminimize the penalty for rial decomposition of 0.5 cm Ca and 20 cm water, the penalty . This factor is higher than the minimum for object suggests that using two bins is far from ideal (because even the is at least 50% for these objects) and minimum penalty far from robust (because using the wrong threshold incurs an even higher penalty). B. Improving Performance The binning penalty and robustness can be improved by increasing the number of bins, but another interesting possibility is removing photons from the spectrum—either at the source or the detector—to increase spectral separation of the two bins. Here we present an idealized scenario where photons in an energy interval are removed without impacting the remainder of the spectrum. We examine the effect of removing photons from the spectrum on the CRLB and the noise penalty. Removing photons from the source spectrum also reduces radiation dose but that effect is not treated here. If we allow for idealized notch filters that remove all photons within a certain range of energies, these gaps in the spectrum should not fall within a bin—rather, they should straddle the binning thresholds since the goal is to increase spectral separation of bins. Then each bin will effectively have lower and higher energy thresholds ( and ), where keV. The measured counts of bin will now . Note that this scenario is a generalization of be the simple binning case, which can be recovered by removing and . all notch filters so that
To optimize our thresholds in regular binning, we perform an exhaustive search over the set of integer such that keV, finding the set of that minimizes , thus enabling the most precise estimate of given bins and incident spectrum . If notch filters are allowed, the search is over the space of all effective thresholds such . If the exhaustive that search is too taxing, e.g., when the number of bins is large, then branching algorithms can be used. To calculate the noise penalty , comparison was made to the CRLB with the full spectrum (i.e., without photons removed) in all cases. Fig. 3 shows the results for four binning strategies: and 3, with and without all energies, for a fixed object . The effective energy ranges of the bins are shown as bars and notch filters are shown as gaps. The column to the right . We find that removing lists the calcium estimate penalty photons in the middle of the spectrum from either of the two , bins, with the notch filter optimized for the measurement of significantly below that of the normal two bin reduced case. Of course we would expect that retaining these photons in a separate bin would be even better. While this indeed provides a further improvement, the additional benefit is marginal because this is an unoptimized three bin design. The optimized three bin design performs even better and has thresholds that are different than those of the case of two bins with a notch filter. As the number of bins is increased, the penalty factor goes to unity. These observations apply to other objects as well, although each object may have different optimal thresholds. In some hardware implementations (e.g., using comparators) inserting a gap between two bins is essentially as complex as adding an additional bin. In those cases, retaining the photons in an (optimized) additional bin is preferred. However, an important point in the results presented above is that much of the benefit of additional bins can come from providing separation in the spectra contained in each bin. If the notch filter is implemented at the source, the benefit to material decomposition also comes with a benefit in dose reduction without additional complexity at the detector. This strategy of applying notch filters is beneficial when es(and also ) but may not apply for other tasks timating where the total number of photons is more important than spectral separation of the bins, such as measuring the total attenuation or forming monoenergetic images. In those cases, retaining the gap photons in a separate bin is a very significant benefit.
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
87
Fig. 4. Energy weighting functions for three special cases: (a) CIX - PC/EI, (b) binning with all energies, and (c) binning with notch filters.
V. WEIGHTED MEASUREMENTS A further generalization of binning counts is forming weighted measurements. The general concept is that photon counting energy discriminating detectors can be used to form weighted measurements of photon counts, where the weighting are a function of photon energy. The functions , where result is a set of measurements
(4)
The concept of energy weighting or weighting functions is not new [8], [16]–[18], and when at least two measurements with different weighting functions are acquired, dual-energy decomposition can be performed. A special case is a counting and integrating X-ray (CIX) detector that does both photon counting and energy integration, providing measurements with weights and [19]–[21] [Fig. 4(a)]. Another special case is binning—the weighting functions are binary and nonoverlapping [Fig. 4(b)]. When notch filters are considered, the weights are zero for energies filtered out, creating separation between the bins [Fig. 4(c)]. In the most general case, the weighting functions are arbitrary and can take on any real value at every energy so that, unlike binning where photons contribute to no more than one bin measurement, here photons can contribute to all weighted measurements. As with binning, forming weighted measurements reduces the amount of data that has to be stored and processed. Since binning, with or without excluding some energies, are special cases, there must be sets of weighting functions that perform at least as well as binning. We seek to characterize the performance as a function of the choice of weights and ultimately to find the optimal weighting functions. For weighted measurements, the random vector (4) is composed of weighted sums of independent, Poisson distributed to meavalues. Every photon of energy contributes value surement . Except for the special case of binary, nonoverlapping weights, since the weighted measurements are obtained from the same realization of , the elements of vector are no longer independent nor Poisson distributed. In such cases, multivariate Gaussian models are commonly used because of their convenient analytical form, tractability, and ability to accurately capture the first- and second-moments of [12], [22], where the weights determine the mean and covariance (B.1). For a sufficiently high transmitted number of photons, this model works well, but when only a few photons are transmitted through
the object, the weighted sum of discrete Poisson distributions may not be well-approximated by a smooth, continuous distribution. Nonetheless, assuming a jointly Gaussian distribution, the CRLB approach can be used to evaluate estimator performance as a function of incident spectrum and energy weights that should be valid at least at high signal-to-noise ratio (SNR) (Appendix B). For comparison, we also find the exact CRLB . The distribution of a weighted for weighted measurements, sum of Poisson random variables is the convolution of the individual scaled Poisson distributions, but has no known tractable can be form. However, the exact distribution and CRLB computed numerically using empirical characteristic functions (Appendix C). A. Optimal Weights In principle, we can find the optimal weights by brute force, treating the weights as variables over which we minimize some objective function [23], [24]. A more elegant approach, if possible, is to find weights that form a sufficient statistic for the . If it could be fully known detected spectrum shown that a set of weights produces weighted measurements that are sufficient statistics for estimating , then the weighted measurements would contain as much information about material thicknesses as the full spectrum . More formally, let be a random vector whose distribution is a sufficient is parameterized by . Then, a statistic is statistic for if the conditional distribution of given independent of [25]. Intuitively, simple examples of sufficient statistics include counting the number of heads in a sequence of coin flips rather than keeping track of the individual flips to estimate the probability of heads or using the sample mean of a sequence of independent, identically distributed random variables to estimate the mean of a normal or Poisson distribution. To find a sufficient statistic, we use the Factorization Theorem is the joint probability mass [25], which states that if function (pmf) of given , then statistic is sufand ficient for if and only if there exist functions such that, for all and . Under this condition can be used to estimate as well as the original measurements . Thus, for our case, the pmf of the detected spectrum given the material thicknesses needs to be decomposed into the product of one function of the statistic and , and a function of the detected spectrum that does not depend on . We begin by writing the joint pmf of the detected spectrum as the product of the pmf at each energy since they are independent given . We use our knowledge that the photon counts are Poisson
88
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
distributed and assume that the incident spectrum and can be treated as a constant
is known
(5) where
(6)
(7) and (8) Therefore, we have proven that two weighted measurements and as their using the attenuation basis functions weights form a sufficient statistic for dual energy imaging. By definition, this sufficient statistic (i.e., these weighted measurements) should achieve optimal performance for estimating any with any incident spectrum . Our only assumption is that the object’s attenuation can be described by two basis functions. We call these weights -weights. VI. COMPARISON To illustrate that the -weighted measurements form a sufficient statistic, we can compute the CRLB. We expect the performance of the -weighted measurements to be identical to that for all objects. In of the full detected spectrum addition, because the -weighted measurements only require two values derived from the detector signals, for comparison we compute the CRLB for two other schemes that report two measurements: optimal binning with two bins and using a CIX detector that does both photon counting and energy integration. A. Performance Curves cm, and Consider a range of relevant ( cm). For the 120 kVp spectrum [Fig. 1(a)], the minvariance ( , in units of cm ) as a funcimum achievable tion of the amount of water and calcium is shown in Fig. 5(a) as
log contour lines. The uncertainty increases with increasing ob(the ject attenuation. Fig. 5(b) shows the penalty function ) when two bins are used with a increase in the variance of fixed threshold of 57 keV, the threshold that minimizes the maximum (over ) penalty factor, and Fig. 5(c) shows similar data for the CIX detector. Binning with a threshold of 57 keV incurs a penalty by a factor of 3.4 at the lower left corner and by at least a factor of 1.7 throughout. A CIX detector incurs lower penalties, especially at larger but still incurs a penalty of . at least 22% and up to a factor of two in Var to for -weights when using the The ratio of Gaussian model is shown in Fig. 5(d) for . This ratio is near 1 at low object thickness, as expected, but drops below 1 for thicker objects where the number of photons is low. However, we know that any unbiased estimator cannot truly have better performance than an estimator having knowledge of the full detected spectrum. Therefore, the multivariate Gaussian model may not be accurate for weighted measurements when the number of photons detected is low. For example, the expected number of photons detected for the thickest object is only 39. Poisson distributions with low mean counts can only take on a set of discrete values and are positively skewed, unlike Gaussian distributions. This conclusion is supported by the fact that the exact CRLB of the -weighted measurements from the empirical characteristic function (Appendix C) is well within 0.1% of the CRLB of the full detected spectrum [Fig. 5(e)]. The accuracy of this method simply depends on the discretization of the Fourier space variable and becomes a numerical issue. These results show that -weighted measurements achieve as knowledge of the the same performance in estimating full spectrum. Although not shown here, the and for -weighted measurements also match those of the full spectrum, while those of binning and the CIX exhibit counterparts under the similar penalty ratios to their different measurement schemes. It is possible to have objects that decompose into a negative amount of calcium or water (but not both). While only quadrant of the calcium/water decomposition space I is shown in Fig. 5, our results extend smoothly into the valid regions of quadrants II and IV, and none of our methods preclude negative amounts of either basis material. B. Wedge Phantom We also simulated the dual-energy performance of the three different binning or weighting schemes on a phantom designed to compare the calcium detectability of the methods. Consider a projection image with 0.2 mAs exposure of a water wedge phantom varying from 0 to 40 cm thickness in the horizontal direction with square calcium contrast elements ranging in area from 1 mm to 1 cm overlaid on top [Fig. 6(a)]. These elements have a thickness such that the predicted (full area) SNR of each calcium element from the full detected spectrum is 4 (at the threshold of detectability). Because the necessary amount of calcium to maintain a constant SNR increases as the water gets thicker and as the elements becomes smaller, the calcium images are displayed on a spatially varying grayscale window in the horizontal direction, where is specified as a function of horizontal position [Fig. 6(b) and (c)].
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
89
Q
Fig. 5. CRLB performance over a wide range of material thicknesses for (a) the full detected spectrum (cm ) and the multiplicative penalty factor for (b) two bins with a 57 keV threshold (no notch), (c) counting and integrating, (d) optimal -weights with a Gaussian model, and (e) optimal -weights with the exact = 57 keV. (c) numerically computed CRLB. Note the different contour levels in the subfigures. (a) Var(^ ): Full spectrum. (b) Var(^ ) penalty: 2 bins, Var(^ ) penalty: CIX. (d) Var(^ ) penalty: -weights, Gaussian. (e) Var(^ ) penalty: -weights, Exact.
t
t
t
The numerical phantom was designed so that all contrast elements project onto an integer number of detectors and are perfectly aligned with the detector grid so as to avoid any partial volume effects. Furthermore, while it is a practical concern, we leave aside the issue of spatial resolution, count rates, and energy resolution by assuming an idealized photon counting detector with no cross-talk, no count rate losses, and perfect energy resolution. Transmission of the X-ray spectrum and detection with energy resolution of 1 keV and perfect detection efficiency was simulated. Quantum noise was incorporated by independently sampling the appropriate Poisson distribution at every 1 keV energy. The simulated noisy detected spectra were then subjected to each of the measurement schemes so that all binned and weighted measurements were obtained from the same realization of detected spectra. Maximum-likelihood estimation (MLE) is a commonly used method for estimating given noisy measurements [9], [26], [27] and was a convenient choice given the extensive use of the log-likelihood functions (A.2) and (B.1) in deriving the CRLB. We used Matlab’s Optimization Toolbox (v3.1.1) to find the maximum-likelihood solution
(9) Although more sophisticated methods may exist for estimating , these are beyond the focus of this paper. For the CIX and weighted measurements, we used the Gaussian model. We compare the ML decomposition of measurements taken with: 1) two
t
t
bins with a 57 keV threshold (no notch); 2) a CIX detector (i.e., ); 3) -weights. The resulting ML decompositions for the calcium component are shown in Fig. 6. Detection of the targets is challenging in all the images since the phantom was designed to make this so, but the -weights image is noticeably superior. The images illustrate that using -weights to form two measurements allows for increased conspicuity of all contrast elements as compared with a two bin or CIX approach. The estimates are shown with . As the spatially varying grayscale window expected from the CRLB comparisons (Fig. 5), the noise in the calcium image is highest when using two bins and lowest when using -weights. The empirical SNR as measured with 10 000 realizations of this experiment is consistent with the theoretical SNR using the CRLB (Fig. 5), although a bias in the estimates becomes evident at the thick end of the wedge. As predicted by Fig. 5(c), the advantage of -weights over CIX is largest in the lower left corner of the wedge and drops off to the upper right. VII. DISCUSSION AND CONCLUSION A simple explanation for the improvement of binning performance from including a notch between bins is that the sensitivity of the ML estimate to detected photons depends on the photon energy. Fig. 7 plots the expected transmitted spectrum and the change in the material decomposition estimates using the full detected spectrum when one additional photon of energy is detected beyond the expected number of counts for the object . Photons near 64 keV provide very little information to the calcium estimate, whereas additional photons at lower energies would induce the MLE to predict less
90
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
0
Fig. 6. Calcium thickness estimation task for overhead projection with different weighting schemes, displayed on a spatially varying grayscale range [ (x) + x)], where (a) is an illustration of the water wedge with calcium contrast elements of different sizes and thicknesses overlaid, (b) is the true calcium thickness, (c) shows the grayscale window as a function of horizontal direction x, and the estimated calcium thickness is from (d) two bins with a 57 keV threshold (no notch), (e) counting and integrating, and (f) optimal -weights. (a) Wedge phantom. (b) True t . (c) Spatially varying . (d) t^ : 2 bins, = 57 keV. (e) t^ : CIX. (f) t^ : -weights.
(
Fig. 7. The expected transmitted spectrum for object t = (0:5; 20) and the sensitivity of estimates t^ and t^ to the number of photons at energy j when using the full detected spectrum. The shaded region is the range of photon energies that should be removed when estimating t with two bins.
calcium and more water and vice-versa for additional photons at higher energies because the ratio of calcium to water attenuation is much higher at lower energies. While photons in the energy notch contain some information, once binned, photons of different energies within bin lose their individual energy identities and cannot be processed any differently than other photons in the same bin that have much more information. They therefore weaken the overall material discrimination ability. Material decomposition is not only insensitive to photons in the shaded energy range (Fig. 7), but these photons also move the average energy of each bin closer to each other and reduce the spectral separation between the bins. It has been shown for
dual kVp and dual layer methods that spectral separation is desired [3]–[7]. Introducing a notch increases the spectral separation. While this is at the expense of fewer photons per bin, it nonetheless improves the precision of the calcium (and water) estimate. In analogous prior work on dual energy imaging with dual layer detectors, it was found that adding a metal filter between the front and back layers improves material separation performance. In that case and as in our work, discarding photons of intermediate energy improves the precision of material decomposition even though it decreases the number of photons that contribute to the measurements [5]–[7]. Ideally, X-rays of the discarded energies would be removed from the incident spectrum to reduce radiation dose, and the net effect is equivalent to binning with a notch between bins. If a notch filter is realizable (either at the source or detector), material decomposition precision can be improved with the same data storage, transmission, and estimation requirements as without one. In practice such selectivity may not be possible to achieve, but the large fraction of omitted photons (Fig. 7) highlights the important role that source filtration and optimization can play in dose efficiency. It is important to note that the implementation of two bins with a notch at the detector using comparators requires the same number of comparators as using three bins. Without a doubt, forming a third bin from the discarded photons instead and utilizing it enables better material decomposition estimation. Although adding more bins and notch filter(s) improves the calcium estimation task performance, the optimal thresholds are a function of object size and the incident spectrum. No fixed binning scheme is universally optimal. Therefore, binning cannot
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
Fig. 8. Two sets of equivalent weights where the weights on the left are simply the linear attenuation functions of calcium, water and iodine and the weights on the right are linearly independent combinations of the attenuation functions. (a) Weights ; , and (b) Linear combinations.
form a sufficient statistic for the full detected spectrum, although for a large number of bins, the difference in performance may be negligible. Using -weights, the detector performance is optimized for all spectra and objects, and thus we have decoupled the joint optimization of the spectrum and detector weights into simply an optimization of the spectrum. The -weighted measurements could perhaps be formed at the detector using analog processing with two parallel circuits and [i.e., an energy dependent that map pulse energy to gain of and ] followed by integrators that and , respectively, or digitally by determining represent the energy of a detected photon from its pulse height and adding and to two different accumulators that represent and , respectively. The values of and can be embedded in hardware or stored in a look-up table. Incorporating the -weights directly into the detectors may require significant hardware redesign. Alternatively, the -weights can be applied to binned counts (assuming there are multiple bins) to reduce the dimensionality of the measurements to only two, although this effectively imposes weights that are a piecewise constant approximation to the attenuation curves and would have to be further investigated. Because any linearly independent comand , the attenuation due to the photoelectric binations of effect and Compton scattering, respectively, can be used as and , there are infinitely many selections for and available. Should the object possibly contain an additional material that deviates from the dual basis material with attenuation assumption, such as a material with a K-edge within the detected spectrum, the attenuation of this material can be used as an additional basis function. The sufficient statistic proof (5) can be extended to show that only three weighted measureand . To illustrate the ments are needed, with weights previous two points, the selected weights could be and [Fig. 8(a), where material 3 is iodine, density: 4.93 g/cm ], or they could be a more complicated linear combination, such as normalized and flattened weights, where the third set of weights has values near 0 for energies greater than 60 keV [Fig. 8(b)]. The performance would be identical since any set of three weighted measurements is entirely recoverable from any other set of three weighted measurements through a simple linear transformation, even in the presence of quantum noise. Using the fact that -weighted measurements form a sufficient statistic, we have shown that a multivariate Gaussian
91
model for the distribution of weighted measurements may not be accurate at low photon counts (the CRLB is lower than possible). This suggests that using a Gaussian model for the estimation of the thicknesses can yield biased estimates. The loss of accuracy in a Gaussian model is not limited to just -weights and extends to CIX detectors and binning because the underlying problem is that the measurements are a weighted sum of Poisson random variables instead of the commonly assumed Gaussian distribution. Although bias correction methods and efficient estimators are beyond the scope of this paper, we expect them to be important issues in accurate and SNR-efficient dual energy decomposition. One solution may be to select a more appropriate likelihood function, such as a Gaussian approximation when the number of photons is large and switching to a discrete distribution at lower counts using the empirical characteristic function. Future work will also focus on other practical considerations of forming sufficient statistics in photon counting detectors, such as incorporating spectral response functions [9], [12] and forming weighted sums of binned counts [8], [10], [17], [18]. In conclusion, the choice of binning and weighting schemes can greatly affect the performance of dual-energy material decomposition. Even with optimal thresholds, binning incurs penalties to decomposition precision and is not robust to changes in the source spectrum or object size, although this can be mitigated by adding more bins. On the other hand, two -weighted measurements in theory form a sufficient statistic for recording the detected photons and in our simulated experiment do well when implemented with a multivariate Gaussian approximation to the true probability mass function. The -weighted measurements outperform binned or CIX measurements in our contrast wedge phantom simulation. Estimating two unknown parameters is at the core of dual-energy theory, and we have proven that only two -weighted measurements are sufficient to accomplish this without any loss of information.
APPENDIX A BINNING CRLB To find the CRLB, we compute the Fisher information matrix, whose elements are given by
(A.1)
is the log-likelihood function of acquiring meawhere surements for object and the expected value of the argument is over , a random variable. Then, if is an unbiased estimator of , the CRLB provides a lower bound on the covariance matrix of by stating that , where . For binned measurements, each bin is the sum of independent Poisson random variables. Therefore, each bin count is itself a Poisson random variable, parameterized by for the case when all energies are used and when
92
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 1, JANUARY 2011
notch filters are used. For measurements of an object characterized by thickness vector , the log-likelihood function (ignoring the constant term) is (A.2) As shown in [12], the resulting CRLB for binning,
The characteristic function of a distribution is the Fourier transform of its probability density function. For a random variable having a Poisson distribution with mean , the characteristic function is
, is
(C.1) (C.2)
(A.3) . One interpretation
where
APPENDIX C EMPIRICAL CHARACTERISTIC FUNCTIONS
is that each expected measurement is weighted by of that is representative of over energies to . In the some special case of binning when we use the full detected spectrum , the CRLB for the full spectrum is
. The characteristic function has where in this section the property that if the random variable is multiplied by a constant , then random variable will have a characteristic . Furthermore, the characteristic function function of a sum of independent random variables is equivalent to the product of their characteristic functions [29]. That is, random variable
(C.3) In the case of two weighted measurements, we can also consider the characteristic function of the multivariate random variable . For
(C.4)
(A.4) where
has a characteristic function
Note that because the characteristic function of a random variable is the Fourier transform of its distribution
.
(C.5)
APPENDIX B GAUSSIAN CRLB taken with weighting For measurements functions , the expected value of mea. In addition, , the covarisurements will be . ance matrix of measurements will be Ignoring the constant term, the log-likelihood function for a Gaussian model of weighted measurements is
This can be numerically computed by discretizing and evalat these points, giving what is known as the empiruating ical characteristic function [30], [31]. Then, taking the 2-D inverse discrete Fourier transform of the empirical characteristic . The resofunction gives us the empirical distribution of lution of depends on the sampled range of , while the range of depends on the spacing of the samples of [32]. Expanding this method, we can also compute the Fisher matrix exactly
(B.1) The elements of the Fisher information matrix are derived in [12], [28] and given by (C.6)
(B.2) Therefore, the CRLB for weighted measurements when modeled by a Gaussian distribution is .
where
is the support of . Then the exact CRLB comes from . We have shown that can be found using (C.5), and the other terms in the integral can be numerically computed as well. These terms all contain partial derivative(s), but because the variables and are the Fourier transform duals
WANG AND PELC: SUFFICIENT STATISTICS AS A GENERALIZATION OF BINNING IN SPECTRAL X-RAY IMAGING
of one another, the partial derivative with respect to can move inside the Fourier transform. For example,
(C.7) The other partial derivatives can be similarly computed numerically for any . Then, because a discretized form of each term in the integrand of (C.6) can be found, the Riemann integral can be found by summing the expression over all samples . to arrive at
REFERENCES [1] R. E. Alvarez and A. Macovski, “Energy-selective reconstructions in X-ray computerized tomography,” Phys. Med. Biol., vol. 21, no. 5, pp. 733–744, 1976. [2] L. A. Lehmann, R. E. Alvarez, A. Macovski, W. R. Brody, N. J. Pelc, S. J. Riederer, and A. L. Hall, “Generalized image combinations in dual KVP digital radiography,” Med. Phys., vol. 8, no. 5, pp. 659–667, 1981. [3] F. Kelcz, P. M. Joseph, and S. K. Hilal, “Noise considerations in dual energy CT scanning,” Med. Phys., vol. 6, no. 5, pp. 418–425, 1979. [4] A. N. Primak, J. C. Ramirez Giraldo, X. Liu, L. Yu, and C. H. McCollough, “Improved dual-energy material discrimination for dual-source CT by means of additional spectral filtration,” Med. Phys., vol. 36, no. 4, pp. 1359–1369, 2009. [5] G. T. Barnes, R. A. Sones, M. M. Tesic, D. R. Morgan, and J. N. Sanders, “Detector for dual-energy digital radiography,” Radiology, vol. 156, no. 2, pp. 537–540, 1985. [6] H. N. Cardinal and A. Fenster, “Theoretical optimization of a split septaless xenon ionization detector for dual-energy chest radiography,” Med. Phys., vol. 15, no. 2, pp. 167–180, 1988. [7] G. M. Stevens and N. J. Pelc, “Depth-segmented detector for X-ray absorptiometry,” Med. Phys., vol. 27, no. 5, pp. 1174–1184, 2000. [8] J. Karg, D. Niederlöhner, J. Giersch, and G. Anton, “Using the Medipix2 detector for energy weighting,” Nucl. Instrum. Meth. A, vol. 546, pp. 306–311, 2005. [9] J. P. Schlomka, E. Roessl, R. Dorscheid, S. Dill, G. Martens, T. Istel, C. Bäumer, C. Herrmann, R. Steadman, G. Zeitler, A. Livne, and R. Proksa, “Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography,” Phys. Med. Biol., vol. 53, no. 15, pp. 4031–4047, 2008. [10] P. M. Shikhaliev, “Tilted angle CZT detector for photon counting/energy weighting X-ray and CT imaging,” Phys. Med. Biol., vol. 51, no. 17, pp. 4267–4287, 2006. [11] A. Macovski, Medical Imaging Systems. Englewood Cliffs, NJ: Prentice-Hall, 1983. [12] E. Roessl and C. Herrmann, “Cramér-Rao lower bound of basis image noise in multiple-energy X-ray imaging,” Phys. Med. Biol., vol. 54, no. 5, pp. 1307–1318, 2009.
93
[13] J. H. Siewerdsen, A. M. Waese, D. J. Moseley, S. Richard, and D. A. Jaffray, “Spektr: A computation tool for X-ray spectral analysis and imaging system optimization,” Med. Phys., vol. 31, no. 11, pp. 3057–3067, Nov. 2004. [14] J. M. Boone and J. A. Seibert, “An accurate method for computer-generating tungsten anode X-ray spectra from 30 to 140 kV,” Med. Phys., vol. 24, no. 11, pp. 1661–1670, Nov. 1997. [15] M. J. Berger, J. H. Hubbell, S. M. Seltzer, J. Chang, J. S. Coursey, R. Sukumar, and D. S. Zucker, XCOM: Photon Cross Section Database (Version 1.3) National Institute of Standards and Technology. Gaithersburg, MD [Online]. Available: http://physics.nist.gov/xcom, Sep. 2009 [16] M. J. Tapiovaara and R. F. Wagner, “SNR and DQE analysis of broad spectrum X-ray imaging,” Phys. Med. Biol., vol. 30, no. 6, pp. 519–529, 1985. [17] R. N. Cahn, B. Cederstrom, M. Danielsson, A. Hall, M. Lundqvist, and D. Nygren, “Detective quantum efficiency dependence on X-ray energy weighting in mammography,” Med. Phys., vol. 26, no. 12, pp. 2680–2683, 1999. [18] J. Giersch, D. Niederlöhner, and G. Anton, “The influence of energy weighting on X-ray imaging quality,” Nucl. Instrum. Meth. A, vol. 531, pp. 68–74, 2004. [19] E. Roessl, A. Ziegler, and R. Proksa, “On the influence of noise correlations in measurement data on basis image noise in dual-energy like X-ray imaging,” Med. Phys., vol. 34, no. 3, pp. 959–966, 2007. [20] E. Kraft, P. Fischer, M. Karagounis, M. Koch, H. Kruger, I. Peric, N. Wermes, C. Herrmann, A. Nascetti, M. Overdick, and W. Ruetten, “Counting and integrating readout for direct conversion X-ray imaging: Concept, realization and first prototype measurements,” IEEE Trans. Nucl. Sci., vol. 54, no. 2, pp. 383–390, Apr. 2007. [21] H. Krüger, J. Fink, E. Kraft, N. Wermes, P. Fischer, I. Peric, C. Herrman, M. Overdick, and W. Rütten, “CIX—A detector for spectral enhanced X-ray imaging in simultaneous counting and integrating,” in Proc. SPIE, 2008, vol. 6913, p. 69130P. [22] M. Firsching, A. P. Butler, N. Scott, N. G. Anderson, T. Michel, and G. Anton, “Contrast agent recognition in small animal CT using the Medipix2 detector,” Nucl. Instrum. Meth. A, vol. 607, pp. 179–182, 2009. [23] A. S. Wang and N. J. Pelc, “Optimal energy thresholds and weights for separating materials using photon counting X-ray detectors with energy discriminating capabilities,” in Proc. SPIE, 2009, vol. 7258, p. 725872. [24] D. Niederlöhner, J. Karg, J. Giersch, and G. Anton, “The energy weighting technique: Measurements and simulations,” Nucl. Instrum. Meth. A, vol. 546, pp. 37–41, 2005. [25] G. Casella and R. L. Berger, “Principles of data reduction,” in Statistical Inference, 2nd ed. Belmont, CA: Duxbury Press, 2001. [26] J. A. Fessler, I. Elbakri, P. Sukovic, and N. H. Clinthorne, “Maximum-likelihood dual-energy tomographic reconstruction,” in Proc. SPIE, 2002, vol. 4684, p. 468405. [27] J. Xu, E. C. Frey, K. Taguchi, and B. M. W. Tsui, “A poisson likelihood iterative reconstruction algorithm for material decomposition in CT,” in Proc. SPIE, 2007, vol. 6510, p. 65101Z. [28] S. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory. Upper Saddle River, NJ: Prentice Hall, 1993. [29] A. Papoulis and S. U. Pillai, “Functions of one random variable,” in Probability, Random Variables and Stochastic Processes, 4th ed. Boston, MA: McGraw-Hill, 2001. [30] A. Feuerverger and R. A. Mureika, “The empirical characteristic function and its applications,” Ann. Stat., vol. 5, no. 1, pp. 88–97, 1977. [31] A. Feuerverger and P. McDunnough, “On some Fourier methods for inference,” J. Am. Stat. Assoc., vol. 76, no. 374, pp. 379–387, 1981. [32] A. V. Oppenheim, R. W. Schafer, and J. R. Buck, Discrete-Time Signal Processing, 2nd ed. Upper Saddle River, NJ: Prentice Hall, 1999.
