A framework for performance characterization of energy-resolving photon-counting detectors
aa r X i v : . [ phy s i c s . m e d - ph ] D ec A framework for performance characterization of energy-resolving photon-countingdetectors
Mats Persson a) Departments of Bioengineering and Radiology, Stan-ford University, Stanford 94305, California, USA(Corresponding author. E-mail: [email protected])
Paurakh L. Rajbhandary
Departments of Electrical Engineering and Radiology, Stanford University,Stanford 94305, California, USA
Norbert J. Pelc
Departments of Bioengineering, Radiology and Electrical Engineering,Stanford University, Stanford 94305, California, USA (Dated: 18 December 2019) 1 urpose : Photon-counting energy resolving detectors are subject to intense researchinterest, and there is a need for a general framework for performance assessment ofthese detectors. The commonly used linear-systems theory framework, which mea-sures detector performance in terms of noise-equivalent quanta (NEQ) and detectivequantum efficiency (DQE) is widely used for characterizing conventional X-ray de-tectors but does not take energy-resolving capabilities into account. The purpose ofthis work is to extend this framework to encompass energy-resolving photon-countingdetectors and elucidate how the imperfect energy response and other imperfectionsin real-world detectors affect imaging performance, both for feature detection and formaterial quantification tasks.
Method : We generalize NEQ and DQE to matrix-valued quantities as functionsof spatial frequency, and show how these matrices can be calculated from simpleMonte Carlo simulations. To demonstrate how the new metrics can be interpreted,we compute them for simplified models of fluorescence and Compton scatter in aphoton-counting detector and for a Monte Carlo model of a CdTe detector with0 . × . pixels. Results : Our results show that the ideal-linear-observer performance for any detec-tion or material quantification task can be calculated from the proposed generalizedNEQ and DQE metrics. We also demonstrate that the proposed NEQ metric isclosely related to a generalized version of the Cram´er-Rao lower bound commonlyused for assessing material quantification performance. Off-diagonal elements in theNEQ and DQE matrices are shown to be related to loss of energy information dueto imperfect energy resolution. The Monte Carlo model of the CdTe detector pre-dicts a zero-frequency dose efficiency relative to an ideal detector of 0.86 and 0.65 fordetecting water and bone, respectively. When the task instead is to quantify thesematerials, the corresponding values are 0.34 for water and 0.26 for bone.
Conclusions : We have developed a framework for assessing the performance ofphoton-counting energy-resolving detectors and shown that the matrix-valued NEQand DQE metrics contain sufficient information for calculating the dose efficiencyfor both detection and quantification tasks, the task having any spatial and energydependence. This framework will be beneficial for the development and optimizationof photon-counting X-ray detectors. 2ey words: X-ray imaging, photon counting, spectral imaging, linear-systems theory,detective quantum efficiency, noise-equivalent quanta, Cram´er-Rao lower bound.3 . INTRODUCTION
The development of energy-resolving photon-counting detectors for medical X-ray imag-ing is an active research field.
Systems with such detectors are currently commerciallyavailable for mammography , and prototypes have been demonstrated for computed to-mography (CT) . An energy-resolving photon-counting detector uses several electronicthreshold levels to separate the registered counts into energy bins depending on the depositedenergy from each event. This enables spectroscopic imaging, i.e. using the information con-tained in the energy distribution of the incoming spectrum to obtain more informationabout the imaged object. The energy information can be used for improving the the signal-difference-to-noise ratio (SDNR) by optimal weighting , for removing beam-hardeningartifacts and for generating material-selective images. The ongoing development of improved energy-resolving photon-counting X-ray detectorsraises the question of how detector performance should be measured. Using a relevant metricis important in order to optimize and compare detector designs comprising different detectormaterials and pixel sizes or with different properties of the readout electronics such as theenergy resolution and anti-coincidence logic.The performance of conventional, non-energy-resolving X-ray detectors, such as energy-integrating detectors and photon-counting detectors without energy-resolving capabilities,is commonly measured in terms of the detective quantum efficiency (DQE). This metric,which is ideally expressed as a function of spatial frequency, measures how well the detectoris able to use the incoming photon statistics for each spatial frequency compared to an idealdetector. A closely related metric is the noise-equivalent quanta (NEQ), which measuresthe number of photons that an ideal detector would require in order to obtain the samesignal-difference-to-noise ratio as the studied detector. An overview of the theory of DQEand NEQ for ordinary (non-spectral) detectors can be found in Ref. 16.Several studies have been published on how to model the frequency-dependent detectivequantum efficiency of photon-counting detectors, but these studies do not take energyinformation into account. At the same time, the impact of detector nonidealities on spec-tral imaging tasks has been studied by several authors . However, these publicationsonly study the zero-spatial-frequency performance and do not investigate the frequency-dependent spectral performance. Shikhaliev et al. studied the effect of characteristic X-rays4n the spectrum shape and on the spatial resolution. Richard and Siewerdsen used linear-systems theory to study the performance of a dual-energy system. Also, the detectabilityfor specific tasks with both spectral and spatial dependence has been studied by Fredenberget al. , Yveborg et al. and Chen et al. A common way to combine information from several spectral channels is to make aweighted sum of the energy-selective images, with weights chosen to be optimal for theconsidered imaging task.
A good measure of the performance of an energy-resolvingdetector should therefore take optimal weighting into account. Taking the spatial-frequencydependence of the signal and noise into account when calculating optimal weights can im-prove detectability in an energy-weighted image, either by calculating an optimized weightfor each image or by letting the weights themselves be functions of spatial frequency .This is equivalent to applying different spatial filters to the different energy bin images be-fore forming the weighted sum, an approach studied in the dual-energy case by Richard andSiewerdsen. Another way to utilize spectral X-ray measurements is to perform basis material de-composition, which builds on the observation that the linear attenuation coefficient of anysubstance in the human body can be approximated well as a linear combination of a smallnumber of basis functions, typically two in the absence and three in the presence of a high-atomic-number contrast agent with a K-edge in the diagnostic energy range. Basis materialdecomposition then amounts to using the energy-resolved X-ray data to estimate the amountof each of these basis materials along each projection line, thereby completely characterizingthe energy-dependent attenuation of the object.
This is useful in particular for CT, toremove beam-hardening artifacts and characterize object composition. Assessment of thepotential performance of X-ray detectors for basis material decomposition is typically madeusing the Cram´er-Rao lower bound (CRLB), which gives a lower bound for the variance ina basis material image estimated from the measured data. Most previous studies of detectorperformance utilizing the CRLB have studied only the zero-frequency (large-area) perfor-mance, although an extension to higher spatial frequencies has been published recently. The purpose of this work is to extend and unify the above-mentioned methods of detectorcharacterization into a general framework for comparing the performance of energy-resolvingimaging detectors, encompassing both energy-weighting and material-decomposition perfor-mance. To this end, we derive an expression for the detectability achieved by an ideal linear5bserver, taking into account the full spatial-frequency dependence of the available energyinformation. We also propose generalized NEQ and DQE metrics for measuring energy-resolving detector performance. Furthermore, we study the performance attainable in mate-rial quantification tasks by deriving an expression for the frequency-dependent CRLB, andshow that the CRLB is closely related to the proposed generalized NEQ metric. As with anyNEQ and DQE formalism, we assume a linear system, although the relevance to situationswhere nonlinear imaging is used is discussed. The outline of the article is as follows: In Sec.II, we derive the expressions for detectability, NEQ, DQE and CRLB and show how theyare related to each other. To show how the proposed metrics can be used in practice, weshow in Sec. III how this framework can be applied to two simple models of detectors withnonideal energy response, and to a more realistic simulation model of a photon-countingdetector. In Sec. IV we present the resulting NEQ and DQE metrics for each of these threemodels, and we discuss the implications of these results in Sec. V and VI.
II. THEORYII.A. Matrix-valued NEQ and DQE
In the following, we use parentheses to denote that a quantity is a function of a continuousvariable and subscripts when the variable is discrete. Most of our notation is similar toRef. 16 and 18. Consider a situation where a distribution of X-ray photons is incident ona detector. Let q ( r , E ) be a random process describing the number of incident photonsper area and energy as a function of position r = ( x, y ) T and energy E . This photondistribution is measured by an energy-resolving detector and registered as counts in energybins d s n ,k where n = ( n x , n y ) T is the discrete coordinate of a detector pixel in the ( x, y ) planeand k = 1 , . . . , N b is the index of the energy bin. To avoid boundary effects, we assumethat the detector and the photon distribution are of infinite extent. Both q and d s arerandom processes, with expected values q and d s , respectively. It is also useful to definea presampling detected signal d k ( r ), whose value at every point r is the number of countsthat would be registered in each energy bin k by a fictitious pixel centered at r . Therefore, d s n ,k = d k ( r n ) where r n = ( n x ∆ x , n y ∆ y ) T is the center of pixel n and ∆ x and ∆ y are thepixel center-to-center distances in the x and y directions, respectively.6he continuous Fourier transform of q ( r , E ) with respect to r is a random process givenby Q ( u , E ) = F [ q ( r , E )] ( u , E ) = R R q ( r , E )e − πi u · r d r , and D k ( u ) is defined analogously asthe transform of d k ( r ). The 2D discrete-space Fourier transform of d s n ,k with respect to n isgiven by D sk ( u ) = F DS (cid:2) d s n ,k (cid:3) ( u ) = P ∞ n = −∞ d s n ,k e − πi u · r n . We now assume that the system islinear and shift-invariant, which allows us to introduce the system point-spread function forthe k th energy bin h k ( ∆r , E ) that relates the expected values of the the incoming photondistribution and the registered count distribution: d k ( r ) = ∆ x ∆ y Z R Z E max h k ( r − r ′ , E ) q ( r ′ , E )d E d r ′ , (1)where E max is the maximum incident energy. Here, the pixel cell area ∆ x ∆ y has beenbroken out of the definition of h k ( ∆r , E ), meaning that h k can be interpreted as a pre-sampling distribution of registered quanta per unit pixel area in the k th energy bin whenthe input signal q ( r , E ) is a beam incident on a single point of the detector. h k ( ∆r , E )thus contains information about both the quantum detection efficiency and the spatial andenergy resolution of the detector. Denoting the Fourier transform of h k with respect to ∆ x and ∆ y by H k , D k ( u ) = ∆ x ∆ y Z E max H k ( u , E ) Q ( u , E )d E. (2)In analogy with conventional X-ray detectors, it is possible to define a presampling mod-ulation transfer function [MTF pre ] k ( u , E ) = | H k ( u , E ) /H k ( , E ) | for each incident energyand for each energy bin.To calculate the performance of the detector for a given imaging task, we also needto know the covariance between different measurements d s n ,k . Assuming that the noise iswide-sense stationary, this is given by the cross-covariance matrix K s with elements K s ∆n ,k,k ′ = Cov( d s n ,k , d s n + ∆n ,k ′ ) . (3)The cross-spectral density of d s n ,k can then be calculated as W sk,k ′ ( u ) = ∞ X ∆ n = −∞ K s ∆n ,k,k ′ e − πi ( u · ∆r n ) , (4)where k and k ′ are indices of two energy bins and ∆r n = r n + ∆n − r n . For k = k ′ , W sk,k ′ ( u )is the noise power spectrum (NPS) in energy bin k , whereas W sk,k ′ ( u ) for k = k ′ describesthe frequency dependence of correlations between different energy bins.7ow assume that the studied detector is used to discriminate between two cases: a object-absent (background) case with incident photon distribution q ( r , E ) and an object-presentcase where the incident photon distribution is q obj ( r , E ) = q ( r , E ) + ∆ q ( r , E ), i.e. a signal-known-exactly/background-known-exactly (SKE/BKE) task. The difference in the expectedpresampled signal between the two cases is ∆ d n ,k with Fourier transform∆ D k ( u ) = ∆ x ∆ y Z E max H k ( u , E )∆ Q ( u , E )d E, (5)where ∆ Q ( u , E ) is the expected signal difference in the Fourier domain. The expecteddifference in the sampled signal is ∆ d s n ,k with discrete-space Fourier transform ∆ D sk ( u ) = x ∆ y P ∞ m u ,m v = −∞ ∆ D k ( u − m u ∆ x , v − m v ∆ y ).Our measure of the detector performance for this task will be the squared signal-difference-to-noise ratio (“detectability”) d ′ of the optimal linear observer, which tellsus how well this model observer can discriminate between the two cases (object present orabsent) given the available measurements in all pixels and all energy bins. Viewing d s n ,k and∆ d s n ,k for all k and n as elements of vectors d s and ∆d s (see Sec. 13.2.12 of Ref. 37), d ′ can be expressed as d ′ = (cid:16) ∆d s (cid:17) T Cov( d s ) − ∆d s . (6)We will now transform (6) to the Fourier domain and express this formula in a way thatseparates parameters specific to the detector from parameters specific to the discriminationtask. To accomplish this, we need to assume that signal aliasing is negligible, so that∆ D sk ( u ) ≈ x ∆ y ∆ D k ( u ). This is true if the pixel size is sufficiently small compared toother resolution-limiting factors, such as the detector point-spread function or the focal spotsize, or if the object is inherently band limited, or if an oversampling scheme is used i.e.several measurements are performed for different detector positions relative to the object.To obtain an expression that is easily comparable to non-pixelated systems, we also expressthe noise correlations in terms of the cross-spectral density W d + ( u ) = x ∆ y W s ( u ) of thesampled pulse train signal d + k ( r ) = P ∞ n = −∞ d s n ,k δ ( r − r n ). These manipulations, which arefound in Appendix A, yield d ′ = Z Nyq X k X k ′ Z E max Z E max H k ( u , E ) ∗ ∆ Q ( u , E ) ∗ · (cid:2) W − d + (cid:3) ∗ k,k ′ ( u ) H k ( u , E ′ )∆ Q ( u , E ′ )d E d E ′ d u . (7)8here ∗ denotes complex conjugate and Nyq denotes the Nyquist region n u : | u | < x , | v | < y o .Note that both E and E ′ are incident energies.To rewrite (7) in a form that is more similar to the corresponding expression for a con-ventional detector, we introduce the relative signal difference ∆ S ( u , E ) = ∆ Q ( u , E ) /q ( E )and define a frequency-dependent matrix NEQ ( u ) with elements given byNEQ( u , E, E ′ ) = X k X k ′ q ( E ) H k ( u , E ) ∗ (cid:2) W − d + (cid:3) ∗ k,k ′ ( u ) H k ′ ( u , E ′ ) q ( E ′ ) . (8)This gives d ′ = Z Nyq Z E max Z E max ∆ S ( u , E ) ∗ NEQ( u , E, E ′ )∆ S ( u , E ′ )d E d E ′ d u , (9)which is analogous to the formula for conventional X-ray detectors (Ref. 16, Eq. 2.150, fornon-pixelated detectors): d ′ = Z | ∆ S ( u ) | NEQ( u )d u , (10)where (Ref. 16, Eqs. 2.193 and 2.209)NEQ = q G MTF ( u )NPS dig ( u ) = q G MTF ( u )∆ x ∆ y NPS d + ( u ) = q H ( u )NPS d + ( u ) . (11)In this equation, G is the large-area gain (average registered counts per pixel divided bythe expected input quanta q ), MTF pre is the presampling modulation transfer function and | H ( u ) | = G MTF pre / (∆ x ∆ y ) is a transfer function defined in analogy with H k ( u , E ). NEQ ( u ) is thus a natural generalization of the noise-equivalent quanta used to describeconventional detectors. It is a frequency-dependent matrix with continuous indices E and E ′ encoding all the information about the detector that is relevant for evaluating the detectionperformance. It depends on the background case noise characteristics, and thus on theincident spectrum in that setting, but not on the discrimination task. The discriminationtask function ∆ S determines how the different frequency components of NEQ are weightedtogether to give the detectability. The
NEQ matrix therefore defines a quadratic form in∆ S ( u ) giving the contribution to d ′ at each spatial frequency.To attain the performance limit given by d ′ , the model observer needs to take data fromall the energy bins into account and give an optimal weight to each energy bin at each spatial9requency, which is difficult for a human observer. For any given task, however, frequency-dependent optimal weighting of the energy bin images can be used to form a single image, forwhich an ideal linear observer attains the same performance for that task. This observationallows the following interpretation: d ′ is the maximum achievable detectability in a singleimage that is formed as a weighted sum of the bin images, where the weights are frequency-dependent and optimized for the given detection task. Note that the optimization of theweights ensures that this detectability is always greater than (or equal to) the detectabilityin an image formed by simply adding the bin images together with equal weight (i.e. animage formed from all counts above the lowest threshold).It is important to keep in mind that both E and E ′ represent incident energies. It mayseem surprising at first that NEQ( u , E, E ′ ) depends on two incident energy indices insteadof one, but this is necessitated by the fact that a detector with imperfect energy resolutionmay confuse incident photons of different incident energies with each other (e.g., a photonfor which the full energy is recorded and a higher energy photon for which some of the energyescapes). Any description of the detector performance must therefore be able to describehow well the detector can distinguish between incident photons with any combination ofenergies E and E ′ . The total detectability of a feature is therefore the sum of detectabilitycontributions from all pairs of energies ( E, E ′ ).The above description applies to a continuous-to-discrete imaging system, which is therelevant type of system for multibin photon-counting detectors. Note, however, that the anal-ysis, apart from the sampling step, can be adapted to a continuous-to-continuous imagingsystem by replacing n x , n y and k by continuous variables. An ideal energy-resolving imagingsystem is a continuous-to-continuous imaging system that registers the correct position andenergy for each incoming photon. Its transfer function is therefore H ( u , E, ε ) = δ ( E − ε ),where ε denotes the registered energy, and its noise cross-spectrum is white, i.e. W q ( u , ε, ε ′ )is constant with respect to u . For an ideal photon-counting detector, the input signal isPoisson distributed and the registered counts at different energies are uncorrelated, so that W q ( u , ε, ε ′ ) = q ( ε ) δ ( ε − ε ′ ) where q ( ε ) is the expected number of incident photons per energyand area (Ref. 16, Sec. 2.6.2.3). For a pixelated but otherwise ideal detector, W d + ( u , ε, ε ′ )tends to W q ( u , ε, ε ′ ) as ∆ x and ∆ y tend to 0, which gives NEQ ideal ( u , E, E ′ ) = q ( E ) δ ( E − E ′ )for an ideal detector.In analogy with the conventional definition of DQE, DQE( u ) = NEQ( u ) /q , we can now10efine a DQE matrix by normalizing
NEQ by the ideal-detector performance. Therefore, letDQE( u , E, E ′ ) = q ( E ) − / NEQ( u , E, E ′ ) q ( E ′ ) − / , (12)or equivalently,DQE( u , E, E ′ ) = q ( E ) / "X k X k ′ H k ( u , E ) ∗ (cid:2) W − d + (cid:3) ∗ k,k ′ ( u ) H k ′ ( u , E ′ ) q ( E ′ ) / (13)This definition gives DQE ideal ( u , E, E ′ ) = δ ( E − E ′ ) for the ideal detector. The fact thatno detector performs better than an ideal one is now formulated as NEQ ≤ NEQ ideal where ≤ should be interpreted in the matrix inequality sense, i.e. A ≤ B if B − A is positivesemidefinite. Proving this inequality mathematically would require analyzing how the noisecross-spectral density is related to the signal transfer function, e.g. using cascaded systemsanalysis, which is beyond the scope of this work. We therefore view it as a physicallymotivated requirement that any realistic detector model must satisfy.As a simple example of a nonideal detector, consider a detector that has imperfect quan-tum detection efficiency η ( E ) < u , E, E ′ ) = η ( E ) q ( E ) δ ( E − E ′ ) and DQE( u , E, E ′ ) = η ( E ) δ ( E − E ′ ). Our definition of DQE can thus be seen as a generalization of the energy-dependent detection efficiency η ( E ).The diagonal elements of the DQE matrix can be interpreted as the dose efficiency perenergy for detecting a spectrum change at one single energy, whereas the off-diagonal termsreflect the degree to which different detected energies interfere with each other. Thisinterference may be constructive or destructive, depending on whether the off-diagonalterms are positive or negative. To understand why the
DQE matrix has dimensions ofinverse energy unlike the conventional DQE which is dimensionless, note that a task func-tion which equals ∆ S ( u ) = ∆ S in a small interval ∆ E around a single energy E gives d ′ ( u ) ≈ q ( E )DQE( u , E , E )∆ E ∆ S , which can be compared to d ′ ( u ) ≈ q ( E )∆ E ∆ S for the ideal detector. DQE task = d ′ ( u ) /d ′ ( u ) ≈ DQE( u , E , E )∆ E is thus proportionalto ∆ E if DQE( u , E , E ) is nonsingular. This is the case for binning detectors since thesehave difficulty detecting a change in a narrow energy band measured in a background ofnoise contributions from a wide energy interval.11 I.B. NEQ and DQE in material basis
The
NEQ matrix contains information about the performance of the detector for anyspectral discrimination task. In practice, this matrix would be cumbersome to report, e.g.in a publication, especially if the energy scale is discretized into a large number of steps. Forexample, if the energy variable is discretized in 1 keV-steps from 0 to 140 keV, one needsto supply a 140 ×
140 matrix as a function of spatial frequency in order to characterizethe detector performance. However, the
NEQ matrix can be expressed more compactly byexploiting the basis material concept i.e. the approximation that all substances likely to bepresent in the field of view have attenuation coefficients µ ( E ) in a low-dimensional space,spanned by a small number N m of basis functions f l ( E ): µ ( E ) = P N m l =1 a l f l ( E ). As longas one is only interested in the detector performance for differentiating materials whoseattenuation coefficients lie within this space, and as long as the signal difference is smallenough that a linearized description may be used, one only needs to know the restrictionof the
NEQ matrix to this subspace to be able to predict the performance for the relevanttasks.For a homogeneously illuminated object and in the absence of detected scatter from theobject, we can assume that the photon distribution incident on the detector is given by q ( r , E ) = Φ( E )exp − N m X l =1 A l ( r ) f l ( E ) ! , (14)where Φ( E ) is the number of photons per area and energy incident on the object and A l ( r )are the basis coefficients integrated along the X-ray beam path. We will study the task ofdifferentiating between two situations, A bg l ( r ) and A obj l ( r ) = A bg l ( r ) + ∆ A ( r ), both smallperturbations of a homogeneous baseline A l ( r ) = A l . In the small-signal approximation, q obj ( r , E ) ≈ q ( E ) + P N m l =1 ∂q ( E ) ∂A l (cid:0) A obj ( r ) − A l (cid:1) , and similarly for q bg ( r , E ). Here, q ( E ) is theexpected photon flux at the detector, given by (14) for A l ( r ) = A l , l = 1 , . . . , N m and thederivative ∂q ( E ) ∂A l is evaluated for A l = A l . The difference in expected photon flux betweenthe two cases is approximately ∆ q ( r , E ) = P N m l =1 ∂q ( E ) ∂A l ∆ A l ( r ). Letting ˜ A ( u ) denote thecontinuous Fourier transform of A ,∆ Q ( u , E ) = N m X l =1 ∂q ( E ) ∂A l ∆ ˜ A l ( u ) (15)12nd∆ S ( u , E ) = N m X l =1 L B l ( E ) 1 q tot ∂q tot ∂A l ∆ ˜ A l ( u ) = N m X l =1 L B l ( E )∆ S B l ( u ) , (16)where we have defined L B l ( E ) = q tot q ( E ) ∂q ( E ) ∂A l / ∂q tot ∂A l with q tot = R E max q ( E )d E . L B l ( E ) can beregarded as an element of a transformation matrix from the basis of individual energies tothe basis B of differential spectrum changes corresponding to differential path length changesof basis materials { f , f , . . . , f N m } . We have also defined the signal difference vector in thebasis B , ∆S B ( u ), by ∆ S B l ( u ) = q tot ∂q tot ∂A l ∆ ˜ A l ( u ). Eq. 9 can now be expressed as d ′ = 1 q tot2 Z Nyq X l X l ′ ∆ ˜ A l ( u ) ∗ ∂q tot ∂A l NEQ B l,l ′ ( u ) ∂q tot ∂A l ′ ∆ ˜ A l ′ ( u )d u (17)= Z Nyq ∆S B ( u ) † NEQ B ( u ) ∆S B ( u )d u , where the matrix elements of the NEQ matrix in basis B are given byNEQ B l,l ′ ( u ) = Z E max Z E max L B l ( E )NEQ( u , E, E ′ ) L B l ′ ( E ′ )d E d E ′ , (18)or, in a form that is easier to compute:NEQ B l,l ′ ( u ) = X k X k ′ H B k,l ( u ) ∗ (cid:2) W − d + (cid:3) ∗ k,k ′ ( u ) H B k ′ ,l ′ ( u ) . (19)Here, H B k,l ( u ) = R E max E =0 q tot (cid:16) ∂q ( E ) ∂A l / ∂q tot ∂A l (cid:17) H k ( u , E )d E is the transfer function, relating themagnitude of a small relative modulation of the input photon fluence corresponding to achange in A l at spatial frequency u to the resulting change in registered counts per area inenergy bin k at that frequency.Like the NEQ for non-energy resolving detectors, NEQ B has units of area − . Each diagonalcomponent NEQ B ll of the latter specifies the number of quanta per area that an ideal detectorneeds to measure to achieve the same detectability as the studied detector, when the taskis to detect the addition of a small amount of basis material l . The DQE matrix in basis B is defined in analogy with (12):DQE B l,l ′ ( u ) = NEQ B l,l ′ ( u ) q NEQ B , ideal l,l · NEQ B , ideal l ′ l ′ , (20)13here our notation reflects that NEQ B , ideal is independent of spatial frequency. The DQE defined in this way is a dose-normalized quantity that encodes the performance of the de-tector for all different tasks. For any specific detection task, the performance relative to anideal detector is given by the task-specific DQE, which is a scalar-valued quantity:DQE task ( u ) = d ′ ( u ) d ′ ( u ) = R E max R E max ∆ S ( u , E ) ∗ NEQ( u , E, E ′ )∆ S ( u , E ′ )d E d E ′ R E max | ∆ S ( u , E ) | q ( E )d E (21)= ∆S B ( u ) † NEQ B ( u ) ∆S B ( u ) ∆S B ( u ) † NEQ B , ideal ( u ) ∆S B ( u ) , where we have used the notation d ′ ( u ) for the d ′ contribution from frequency u . Eq.21 shows that DQE task for detecting a small difference at frequency u in one of the basismaterials used in the decomposition, basis material l , is simply given by the diagonal elementDQE B l,l ( u ) of DQE B . The off-diagonal terms of DQE B l,l , on the other hand, specify the degreeof constructive or destructive interference when the task involves detecting a difference in twoor more basis functions. A positive off-diagonal term corresponds to a positive interferenceeffect on detectability when both basis material path lengths are increased simultaneously.DQE task for detecting differences due to a material with another spectral response can beobtained by expressing the DQE matrix in any basis which includes the linear attenuation ofthat feature as a basis function. When transforming to another basis of differential spectrumchanges, NEQ transforms according to the normal rules for coordinate changes of quadraticforms, but there is no simple transformation rule for
DQE since its components are obtainedas ratios of components of two matrices both of which are basis-dependent.In particular, we note that Eq. 21 can be used to calculate the maximum and minimumDQE task for any detection task where the linear attenuation coefficient of the feature todetect is a linear combination of the basis functions in B . The change of variables ∆S B ′ ( u ) = (cid:2) NEQ B , ideal ( u ) (cid:3) / ∆S B ( u ) transforms (21) intoDQE task ( u ) = h ∆S B ′ ( u ) † NEQ B , ideal ( u ) †− / NEQ B ( u ) NEQ B , ideal ( u ) − / ∆S B ′ ( u ) i / (cid:13)(cid:13)(cid:13) ∆S B ′ (cid:13)(cid:13)(cid:13) , (22)with maximum and minimum values given by the eigenvalues of NEQ B , ideal ( u ) †− / NEQ B ( u ) NEQ B , ideal ( u ) − / .14 I.C. Material decomposition
We now leave the study of detection-task performance for a moment and turn to anotherquestion: how well can the material composition of an object be quantified using measure-ments with the detector? One could try to estimate A l ( r ) directly, but a more fruitfulapproach is to estimate its Fourier transform ˜ A l ( u ). A limit on estimation performance isgiven by the Cram´er-Rao lower bound (CRLB), which gives a lower limit of the covariancematrix of the estimated parameter for any unbiased estimator and has been shown to be auseful tool for analyzing estimation accuracy for spectral CT tasks . Although the CRLBonly provides a lower bound, in practice the variance of the maximum likelihood estimatorusually agrees well with the CRLB for material decomposition tasks, as has been shown insimulations . The authors are not aware of a tractable form for the CRLB for correlatedPoisson variables in general, but as long as there are at least a few tens of photons in everyenergy bin and measurement, one can approximate the joint probability distribution of thebin counts as multivariate Gaussian. This approximation is valid both with and withoutpileup taken into account, but in the pileup-free case both mean and variance of each indi-vidual measurement will be equal to d s n ,k . The components of the Fisher information matrix F θ for a real-valued vector parameter θ from a measured multivariate Gaussian randomvariable with mean µ and covariance matrix Σ is given by: [ F θ ] ij = (cid:18) ∂ µ ∂ θ i (cid:19) T Σ − ∂ µ ∂ θ j + Tr (cid:18) Σ − ∂ Σ ∂ θ i Σ − ∂ Σ ∂ θ j (cid:19) . (23)The CRLB now states that Cov( ˆ θ ) ≥ F θ − for any unbiased estimator ˆ θ of θ .In our case, µ and Σ are the mean and covariance of the measured counts vector d s , whilethe parameter θ = (cid:16) Re ˜ A , Im ˜ A (cid:17) T to be estimated contains the real and imaginary partsof the vector ˜ A obtained by concatenating the vectors ˜ A ( u ) for all spatial frequencies in theintersection of the right half-plane { u : u > u = 0 , v > } and the Nyquist region. Since˜ A is the Fourier transform of a real-valued function and hence conjugate-symmetric in u ,we only estimate it on half the Nyquist region in order to avoid getting a singular covariancematrix. (Readers who are skeptical about our use of the CRLB derived in a finite-dimensionalsetting for estimating a function in an infinite-dimensional space can think of the function˜ A ( u ) as being approximated by a large but finite number of delta functions at discrete15ample points in the u plane. ˜ A then becomes a vector of delta function coefficients.)In the absence of pileup, we note that both d s and Cov( d s ), and therefore also the firstterm in (23), are proportional to the pre-patient photon flux, whereas the second term isindependent of this quantity. As long as the photon fluence is high enough, we can thereforeapproximate the second term with zero. We then obtain the Fisher information matrix F θ for θ as F θ ≈ ∂ d s ∂ θ ! T Cov( d s ) − ∂ d s ∂ θ . (24)As shown in Appendix B, this leads to the following expression for the CRLB for thecomponents of ˜ A ( u ):Cov (cid:16) ˆ˜ A l ( u ) , ˆ˜ A l ′ ( u ′ ) (cid:17) ≥ q tot2 (cid:2) NEQ B ( u ) − (cid:3) l,l ′ ∂q tot ∂A l ∂q tot ∂A l ′ δ ( u − u ′ ) . (25)Here, the complex covariance is defined as Cov( z, w ) = E [( z − E( z ))( w − E( w )) ∗ ] where Edenotes expected value. Eq. 25 specifies the achievable variance and covariance in the basisimages, and thereby answers our question about the material quantification performance. Itis a natural extension of the commonly used zero-frequency CRLB to frequency-dependentestimation tasks, which allows us to describe the full noise correlation structure in thematerial-decomposed images, rather than just the pixel-wise variance.In analogy with the definition (21) of DQE task for a detection task we can define theDQE task for quantification of basis material l asDQE task ( u ) = Var( ˆ˜ A ideal l )Var( ˆ˜ A l ( u )) , (26)where Var (cid:16) ˆ˜ A ideal l (cid:17) is the CRLB of the variance of ˜ A l when measured with an ideal detector.Eq. 25 shows that NEQ B ( u ), which describes performance for detection tasks, and thefrequency-dependent CRLB, which describes the detector performance for material quan-tification tasks, are very closely related. We will now conclude our theory presentation witha demonstration of another close connection between the detection and material quantifi-cation frameworks. To this end, we note that it is possible to make a frequency-dependentoptimally weighted sum of the decomposed basis images, just as this can be done with16he original energy bin images, and we then investigate what detection performance can beachieved in the weighted image if the weights are chosen optimally. The lower bound forthe variance in the basis images given by the CRLB translates to an upper bound of theachievable detectability in any such weighted basis image. This upper bound is given by theoptimal-linear-observer detectability and can be expressed as (see Appendix C): d ′ ≤ Z R ∆S B ( u ) † NEQ B ( u ) ∆S B ( u )d u . (27)By comparing (27) with (17), we see that the same detectability is obtained by com-bining the different basis images as from combining the original energy bin images, as longenough statistics is available that the estimator variance will be close to the CRLB. This isa generalization of the analogous result for zero frequency derived by Alvarez. We there-fore conclude that the process of basis material decomposition preserves the informationcontained in the image, as long as the variance of the estimator is close to attaining theCRLB.We emphasize that while this optimal weighting of basis images is theoretically interest-ing, in many situations there are other ways of using the basis images that are preferable.For example, the estimated basis images can be viewed directly as they are, e.g. as contrastagent distribution maps or virtual non-contrast (VNC) images. Although each individualbasis image has lower detectability for specific tasks compared to a frequency-dependentweighted sum of all basis images, the individual images may be more useful in practicesince they provide information about material composition and may provide higher lesionconspicuity. Another way to use the basis images is to generate virtual monoenergetic im-ages (VMI) as p mono ( E ) = P N m l =1 A l ( r ) f l ( E ). This is useful because the image valuesare easy to interpret. Note, however, that the basis images in this case are combined usingfrequency-independent weight functions, meaning that there is no guarantee that there isan energy for which the resulting monoenergetic image has optimal detectability. III. MATERIALS AND METHODS
To illustrate the matrix-valued NEQ and DQE measures, we will study their form for twosimple models (one modeling K-escape and one modeling Compton scatter in the detector),and also for a more realistic simulation model of a CdTe detector. Each of these three17odels is described in one subsection (III.A-III.C) below, and the results for each model arepresented and discussed in corresponding subsections of Sec. IV and V.
III.A. Fluorescence model
As a simplified model of fluorescence K-escape, we assume that an incident photon ofenergy E will deposit either its entire energy, with probability 1 − p F ( E ), or an energy E − E F with probability p F ( E ). E F is thus the energy of the fluorescence photon, and the probabilityof fluorescence escape p F ( E ) is a function of E with p F ( E ) = 0 for E ≤ E F . For simplicity,we neglect detection of secondary (fluorescence) photons in other detector elements, e.g.all fluorescence photons escape or the detector is assumed to have coincidence logic thateliminates secondary events but is not able to reconstruct the original photon energies. Wealso assume that the detector has infinitely many bins, infinitely small pixels and perfectenergy resolution, i.e. the fluorescence escape is the only degrading factor. The transferfunction is then H ( u , E, ε ) = (1 − p F ( E )) δ ( E − ε ) + p F ( E ) δ ( E − E F − ε ) and an incidentspectrum q ( E ) results in a measured spectrum d ( ε ) = (1 − p ( ε )) q ( ε ) + p ( ε + E F ) q ( ε + E F ). Since each photon is only registered once, the measurements at different energiesare independent random variables, and the cross-spectral density is given by W ( u , ε, ε ′ ) = d ( ε ) δ ( ε − ε ′ ). Inserting these results into (8) gives, after some algebra,NEQ( u , E, E ′ ) = (cid:20) (1 − p F ( E )) d ( E ) + p F ( E ) d ( E − E F ) (cid:21) q ( E ) δ ( E − E ′ )+(1 − p F ( E )) q ( E ) p F ( E ′ ) q ( E ′ ) d ( E ) δ ( E ′ − ( E + E F ))+ p F ( E ) q ( E )(1 − p F ( E ′ )) q ( E ′ ) d ( E ′ ) δ ( E ′ − ( E − E F )) . (28)The DQE can then be obtained from (12).To plot the DQE matrix for a simple example, we discretized energy in steps of 1 keV andassumed E F = 25 keV and p F ( E ) = 0 . E ) for E > E F . For simplicity wealso assumed a rectangular incident spectrum: q ( E ) = q = q tot / ( E − E ) for E < E < E and 0 otherwise. In order to study the effect of spectral overlap, we studied two differentincident spectra: one that is nonzero between E = 40 keV and E = 60 keV, which givesnonoverlapping deposited spectra for the K-escape peak and photopeak events, and one that18s nonzero between E = 40 keV and E = 120 keV, which produces overlap between thephotopeak and K-escape spectra. We also used (21) to calculate DQE task for different taskfunctions ∆ S ( E ) which were, for simplicity, taken to be piecewise constant, equal to either0 or ± ∆ S , as described below. III.B. Scatter model
In order to demonstrate how the proposed framework treats interactions where the energyinformation is lost, we also studied a simple model for Compton scatter in the detector. Inthis model, we assume that an incident photon is either photoabsorbed with probability1 − p S , depositing all its energy, or Compton scattered in the detector with probability p S = 0 .
5. In the latter case it deposits an energy E S which we for simplicity assume tobe fixed at 10 keV independent of the incident energy. We also assume that all scatteredphotons are stopped by blocking lamellae or escape from the detector, so that each photonis counted only once. Like in the fluorescence model, we assumed that the incident spectrumis constant, equal to q = q tot / ( E − E ) between E = 40 keV and E = 120 keV and 0otherwise. Once again, the photopeak and distorted spectra do not overlap. This gives thedetected spectrum as d ( ε ) = p S q tot δ ( E − E s ) + (1 − p S ) q ( E ) and, after some calculations,DQE( E, E ′ ) = p S E − E + (1 − p S ) δ ( E − E ′ ). III.C. CdTe simulation
To study a more realistic case, we used Monte Carlo simulation (pyPENELOPE , basedon PENELOPE ) to model a CdTe detector. This simulation model is similar to the modelstudied previously in Ref. 36. A special version of pyPENELOPE was compiled in orderto include secondary photons. For each 1 keV-step from 1 to 120 keV, a pencil beam ofphotons incident on a 3 mm thick slab of CdTe was simulated, and the location and energyof each photon interaction were recorded.We simulated charge sharing according to the uniform spherical charge cloud model ofTaguchi , where the diameter d of the charge cloud is related to the deposited energy E as d = d ( E/E ref ) / , and we used d = 30 µ m at E ref = 70 keV. By dividing the detectorinto 0 . × . pixels and calculating the charge cloud volume fraction located within19he borders of each pixel, we obtained the total charge collected in each pixel. The chargecontributions from all interactions stemming from each individual photon were then summed,and the total number of registered photons in each of five energy bins (25-40, 41-54, 55-64,65-77 and 78-120 keV) were recorded. We repeated this for a grid of 16 ×
16 positions ofthe pencil beam in one quadrant of the detector cell, extending from the center to the pixelborder in the positive x and y directions. Exploiting mirror symmetry in two dimensions,we obtained the point-spread function h k ( ∆r , E ) on a grid with 30 ×
30 sample points foreach detector cell in a grid of 3 × x and y coordinates.To obtain the autocovariance function, we used the same pencil beam simulation butwith a randomized position of each incident photon to simulate homogeneous illuminationof the center pixel. For each incident energy E , we recorded the number of incident photons N n ,k, n ′ ,k ′ ( E ) leading to a registered count in both energy bin k of pixel n and energy bin k ′ ofpixel n ′ = n + ∆n . The covariance of the number of registered photons N n ,k ( E ) and N n ′ ,k ′ ( E )in these two pixel-bin pairs is given by Cov ( N n ,k ( E ) , N n ′ ,k ′ ( E )) = N n ,k, n ′ ,k ′ ( E ), since thephotons registered in either of N n ,k ( E ) and N n ′ ,k ′ ( E ), but in not both, are independentbetween these two measurements. By summing over n in the grid of 3 × K s ∆n ,k,k ′ ( E ) for each incident energy E . K s ∆n ,k,k ′ ( E ) was calculatedfor − ≤ n x , n y ≤
1, i.e. only correlations between nearest neighbors was taken into account.We symmetrized the autocovariance function using mirror symmetry in x and y and withrespect to interchanging the x and y coordinates.We used discrete Fourier transformation to calculate the transfer function H k ( u , E ) andenergy-dependent cross-spectral density W sk,k ′ ( u , E ). We then calculated the noise cross-spectral density of the bin counts W sk,k ′ ( u ) by combining contributions from different energiesusing a 120 kVp tungsten anode X-ray spectrum with a prepatient photon fluence of 4 · mm − and 12 ◦ anode angle, filtered through 2.5 mm Al (prepatient filtration) and100 mm water. The spectrum was obtained from Spektr 3.0 and the linear attenuationcoefficients were obtained from NIST . Electronic noise and pileup were not included in thesimulation. The NEQ and
DQE matrices were calculated in the basis of monoenergies usingEqs. 8 and 12 and in the basis B = { Water , Bone } using Eqs. 19 and 20. By calculating theeigenvalues of NEQ B , ideal ( u ) †− / NEQ B ( u ) NEQ B , ideal ( u ) − / , we obtained the the maximum20nd minimum DQE task for any detection task. We also calculated the task-specific DQEsfor quantifying the amount of cortical bone and water in a two-basis decomposition usingEqs. 25 and 26. IV. RESULTS
In the below sections IV.A-IV.C, we present the results of the calculations of NEQ andDQE for each of the three models described in Sec. III.A-III.C: the fluorescence and scattertoy models and the CdTe model.
IV.A. Fluorescence model
The results for the fluorescence model are shown in Fig. 1, with Fig. 1(a-c) for thenonoverlapping spectrum case and Fig. 1(d-i) for the overlapping spectrum case. In Fig.1(c,f-i), different task functions are shown together with the diagonal of the DQE matrixand the corresponding DQE task values. The
DQE matrix was discretized by replacing thedelta function with a square of width 1 keV. For the overlapping spectrum case (Fig. 1e),the delta function coefficients in the DQE matrix are 0 . − .
84 on the diagonal and 0 . IV.B. Scatter model
The results for the Compton scatter model are shown in Fig. 2. Fig. 2(a) shows theincident and deposited spectra. The
DQE matrix was discretized by replacing the deltafunction with a square of width 1 keV. Fig. 2(c-d) show task functions for (c) densityimaging and (d) spectral imaging, together with the diagonal of the DQE matrix and theresulting DQE task . IV.C. CdTe simulation
Fig. 3 shows the simulated point-spread functions for the CdTe model for monochromaticbeams of 40, 70 and 100 keV. The corresponding transfer functions H k ( u , E ) are shownin Fig. 4 for u ranging from u = 0 mm − to three times the Nyquist frequency. The21atrix elements of the cross-spectral density are shown as a function of frequency in Fig.5, for u from 0 mm − to the Nyquist frequency. Since the autocorrelation function ismirror-symmetric in this model, the cross-spectral density is also symmetric, i.e. extends bymirroring in zero and in the Nyquist frequency.The zero-frequency NEQ and
DQE matrices in the basis of monoenergies are displayedin Fig. 6 together with plots of their respective diagonals. Fig. 7(a) shows the componentsof the
NEQ B matrix for the CdTe model and for an ideal detector, as a function of spatialfrequency for basis materials water and bone. The unattenuated photon fluence of 4 · mm − before the object gives 4 . · mm − after the object. The corresponding DQE B is shown for the same basis functions in Fig. 7(b), whereas Fig. 7(c) shows the largest andsmallest DQE task for any detection task. Finally, Fig. 7(d) shows DQE task for quantifyingwater and bone in a two-material decomposition, calculated with Eq. 26. V. DISCUSSION
In this section, we comment on the results obtained for the fluorescence and scatter toymodels and for the CdTe simulation model (Sec. V.A-V.C). Then, in Sec. V.D we discussthe validity and limitations of the proposed theoretical framework.
V.A. Fluorescence model
In the simple fluorescence model, the NEQ matrix (28) and the corresponding DQEmatrix are independent of spatial frequency, nonzero along the diagonal for the energiescontained in the measured spectrum and, if there is spectral overlap, along two off-diagonallines (Fig. 1(b,e)). To interpret this result, note that according to Eq. 9 the detectabilityof a feature with differential signal ∆ S ( u , E ) at frequency ( u ) is obtained as a sum ofcontributions ∆ S ( u , E ) ∗ NEQ( u , E, E ′ )∆ S ( u , E ) from all pairs of energies E, E ′ . The firstterm of (28), which contains the diagonal elements of the matrix, shows that a relative signaldifference at energy E contributes to d ′ with two terms. The first term is the contributionfrom the non-fluorescence (photopeak) interactions at energy E and the second term is thecontribution from the fluorescence interactions with original energy E and registered energy E − E F . The two off-diagonal contributions in (28) are nonzero only for E ′ = E + E F and22 ′ = E − E F , respectively. These terms contain the contribution from overlap between thephotopeak and the K-escape peak. The denominator for each of these contributions is thequantum noise variance, equal to the number of registered events at the deposited energy.When the incident spectrum is narrow enough that the photopeak and K-escape spectrado not overlap (Fig. 1(a-c)), there are no off-diagonal terms and the DQE matrix is a deltafunction along the diagonal with a coefficient of 1. Even though the K-escape photons areregistered with the wrong energy, there is no risk that these will be confused with photopeakphotons, and since there is no ambiguity about the true energy, the K-escape photons donot cause any performance degradation at all and the DQE task for any imaging task is 1.When the spectrum is broad enough to give overlap between the photopeak and K-escapeevents, on the other hand, off-diagonal elements appear in the
DQE matrix and the valuesof the diagonal elements decrease (Fig. 1(d-i)). The diagonal elements of the DQE matrixindicate how much ∆ S ( u , E ) at each energy E alone contribute to the total d ′ , regardlessof what other energies are present in ∆ S ( u , E ). In this case, the diagonal DQE elementsare highest near the maximum energy, where the recorded spectrum does not include anyK-escape component, and near the minimum energy, where the K-escape events from thesephotons are registered below the minimum incident energy and do not have to compete withthe noise from photopeak events.The off-diagonal terms, on the other hand, specify the contribution to d ′ caused by thepresence of a signal at both energy E and E ′ . This contribution can be either positive ornegative depending on the task. Since the off-diagonal terms are positive in this model,the detectability increases if ∆ S has the same sign at both energies and decreases if ∆ S has opposite signs at the two energies. A case when ∆ S has the same sign at E and E ′ can correspond e.g. to a density imaging task, where the total number of registeredphotons is more important than the energy distribution. In this case the simultaneouspresence of photons at both energies in ∆ S is actually beneficial, since the d ′ contributionis proportional to the square of the differential signal, and increasing spectral overlap inthe signal (while keeping the background noise fixed) helps concentrate more signal at asingle measurement energy. Note, however, that d ′ can never exceed the value that wouldhave been obtained if all photons were registered correctly on the first place. On the otherhand, a case where ∆ S has opposite signs at the two energies is a material discriminationtask, e.g. K-edge imaging. Here, the simultaneous presence of energies E and E ′ decreases23he detectability since spectral overlap degrades the ability of the detector to measure truespectral differences.In the model studied here, the off-diagonal terms are nonzero only for E − E ′ = E F .The impact of the fluorescence on detectability is therefore dependent on whether the taskfunction ∆ S ( u , E ) involves energies with spectral overlap or not. This is demonstrated inFig. 1(f-i). When the task is to detect a signal difference in a narrow energy band, theDQE task is simply given by the diagonal value of the DQE matrix for that energy range (Fig.1(f)). This is less than unity because the registered change in photon flux must be detectedin the presence of a larger noise contribution from co-registered K-escape and photopeakevents. On the other hand, when the task is to detect a signal difference across a broadrange of energies, the off-diagonal elements cause constructive interference and DQE task = 1,as for an ideal detector (Fig. 1(g)). For a density imaging task, the off-diagonal elementsthus improve DQE task .For a spectral imaging task, i.e. when trying to detect an attenuation difference betweentwo energy bands, the off-diagonal elements have a different effect (Fig. 1(h-i)). For aspectral task where the spectral bands do not have any overlap due to fluorescence, DQE task is given by the average over the two energy bands of the delta function coefficient along thediagonal of
DQE , i.e. the off-diagonal elements do not affect the result. If the bands arepositioned so that K-escape events from one of them overlaps with photopeak events fromthe other, the off-diagonal elements give negative interference and decrease DQE task . In thiscase, the off-diagonal elements are detrimental, since cross-talk from one energy band to theother makes it more difficult to detect a spectral difference.
V.B. Scatter model
In the Compton scatter model (Fig. 2) half the incident spectrum is registered at itscorrect energy, while half is registered as a peak at 10 keV (Fig. 2(a)) This is reflectedin the form of the DQE matrix (Fig. 2(b)), which is the sum of two contributions. Thephotoelectric interactions give a diagonal with delta functions whose coefficient is 0.5, i.e.half the DQE of an ideal detector. The Compton interactions provide information that aphoton was detected but do not contain any energy information. Their contribution is evenlydistributed as a weak background over the entire square formed by the spectral supports of24 ( E ) and q ( E ′ ). The total area of the diagonal line and the background is the same sinceeach of them are generated from 50% of the photons.For a density imaging task where ∆ S is constant as a function of energy, d ′ is given bythe integral over nonzero part of the DQE matrix, so that the two contributions add and giveDQE task = 1. Since only the total number of registered photons matters in this case, thefact that some of them are registered at the wrong energy does not degrade the performance.For a spectral imaging task on the other hand, the detectability is obtained as an integralof the NEQ over four squares in the (
E, E ′ ) diagram, two near the diagonal with positivesign and two away from the diagonal with negative sign. The uniform background term inthe NEQ matrix therefore gives a net contribution of 0, so that DQE task = 0 .
5, reflectingonly the photoelectric contribution. The Compton events, which are all registered with thesame energy, thus contain as much information as photoelectric events for density imagingtasks but no information at all for a pure spectral task. A real-world imaging task willtypically be an intermediate between these two types, meaning that the Compton eventswill contribute some information but less than the photoelectric events. Also note thatthis model is simplified in the sense that, in a real detector, the deposited energy from aCompton interaction is random with a distribution function that depends on the incidentenergy and the scattering angle, meaning that the Compton interactions do contain someenergy information.
V.C. CdTe simulation
As shown in Fig. 3, the point-spread function h k ( ∆r , E ) of the CdTe simulation modeldepends on the incident energy. When a monoenergetic beam impinges on the interior of apixel, the majority of counts are registered in the energy bin that includes the true energy:bin 1, 4 and 5 for 30, 70 and 100 keV, respectively. For a 30 keV beam, nothing is registeredin the other bins, since pileup is not included in the model and the energies of differentphotons thus cannot be added together. For higher incident energies, some photons aremisregistered in lower energy bins due to charge sharing and fluorescence. For example,a 100 keV beam incident precisely on the border between two pixels will be registered inenergy bin 2 (41-54 keV) since only about half the photon energy, 50 keV, is deposited in thestudied pixel. If the same beam hits inside the pixel but close to the border, a fraction (up25o 20%) of the photons are registered in bin 4 (65-77 keV) since they only lose a minor partof their energy to charge sharing or fluorescence. Finally, just outside the pixel border, themajority of the events are registered in energy bin 1 (25-40 keV) since most of the energyof these photons is deposited in the neighboring pixel.The shape of the point-spread function is reflected in the transfer function H k ( u , E ) (Fig.4). In each case, the transfer function for the energy bin encompassing the incident energyfalls off from a zero-frequency value close to 0.8, signifying that about 80% of the totalnumber of photons are registered in that energy bin, and crosses 0 near twice the Nyquistfrequency, reflecting the nearly rectangular point-spread function. The other energy binsdetect smaller fractions of the photons and fall off more rapidly, indicating that the blurringfrom the inter-pixel cross-talk suppresses high spatial frequencies.The low-pass character of the cross-talk can also be seen from the cross-spectral densityplotted in Fig. 5. As seen in these figures, the diagonal elements of W d + are nearly constantfunctions of spatial frequency, which indicates that the noise in these is uncorrelated betweenpixels. The exception is the lowest energy bin, which has a low-frequency component in itsnoise stemming from the fact that some photons are counted in this energy bin in morethan one pixel. The cross-elements of W d + exhibit strong low-frequency correlations sincephotons are typically registered in one energy bin in the pixel of incidence and in anotherenergy bin in a neighboring pixel. The strongest such correlation is found between energybins 1 and 2. This shows that the correlation structure in the measured image requires ajoint description in terms of both spatial and energy correlations.The zero-frequency NEQ and
DQE matrices in the basis of monoenergies (Fig. 6) exhibita block-like structure with the block borders corresponding to the energy thresholds. Theseplots show that the largest cross-terms, and thereby the strongest coupling between differentenergies, is found within each energy bin, i.e. in the diagonal blocks. Note that the relativeintensity of the blocks resembles the three diagonal bands caused by the loss of a fixed amountof energy in the simplified fluorescence model (Fig. 1(e)). Figs 6(b,d) also show that boththe NEQ and DQE are largest for the energies corresponding to the characteristic X-raypeaks in the spectrum. All energies falling within a certain energy bin are detected againstthe same amount of background noise, and since the detectability is a quadratic functionof the signal difference, a certain relative signal difference ∆ S is gives a larger detectabilitywhen the spectral density is higher. It is therefore easier to detect an attenuation difference26ight at a characteristic X-ray energy of the source than at a neighboring energy. This alsoexplains why the matrix-valued DQE is dependent on the X-ray spectrum shape.The NEQ B and DQE B are plotted as functions of spatial frequency in Fig. 7(a-b). Notethat the NEQ values for the ideal detector are slightly higher than the transmitted photonfluence of 4 . · mm − , in particular for bone, reflecting the higher detectability of an idealenergy-resolving detector compared to an ideal photon-counting detector with no energydiscrimination. The task-specific DQE values for detecting water and bone are obtaineddirectly from Fig. 7(b) as the diagonal elements of the
DQE B matrix. The plots showthat the performance for both detection of bone and soft tissue suffer relatively little fromthe degraded energy resolution caused by charge sharing and fluorescence, since the zero-frequency DQE task is 0.86 for water and 0.65 for bone. The sinc-like decrease towards higherfrequencies is caused by the pixel aperture. The cross-term is positive, meaning that thepresence of both soft tissue and bone in a detection task will cause constructive interference,i.e. increase detectability.The maximum DQE task for any imaging task is 0.94, while the minimum is 0.23 (Fig.7(c)). Both water and bone detection are therefore tasks for which the performance of thestudied detector comes relatively close to the maximum. For material quantification, onthe other hand, the DQE task is severely degraded by the imperfect energy resolution of thedetector (0.34 for water and 0.26 for bone) and comes close to its lowest possible value forany detection task. (Fig. 7(d)). Thus, the mis-registration of photon energies has a limitedeffect on detecting a small change in the amount of water or bone, since this also affects thetotal number of photons, but a large effect on the ability to determine what material causeda loss of flux, which requires good energy resolution.In the model studied here, we have neglected electronic noise for the purpose of makingthe effects of the spectrum shape and nonideal energy response stand out clearly. In areal detector, electronic noise will blur the detected spectrum, so that the transition of thepoint-spread function at the pixel border becomes smoother compared to Fig. 3 and thecharacteristic X-rays are less prominent in the NEQ matrix compared to Fig. 6.27 .D. Applicability of the framework
Since a completely realistic detector model would be mathematically intractable, the pro-posed framework builds on a number of idealized assumptions. For the detectability formulain basis format (17) to be valid, the feature to be detected must be weakly attenuating, i.e.either small or low-contrast. The detector response is assumed to be linear and shift invari-ant, and the noise is assumed to be stationary. These assumptions are not true in generalin CT, since the photon flux can vary greatly between different projection lines. This makesthe noise non-stationary and can give pileup in some projections, meaning that the responsewill be nonlinear and position-dependent. However, our model is valid in regions of the im-age where the photon flux is slowly varying, so that the detector response is approximatelylinear. The primary usefulness of this framework therefore lies in its ability to aid under-standing and optimization of detector performance for model tasks where a feature is to bedetected against a slowly varying background. On the other hand, imaging of phantoms,in either simulation or experiment, will be needed in order to measure the performancefor more complex tasks. Also note that the assumption that aliasing is negligible meansthat one must use caution when using the framework to predict detection performance forhigh-frequency tasks.There are also effects that are not included in the work presented here but could beincluded in the future. One example is the impact of scatter from the object that is notincluded here but can be modeled as additive noise. Also, our present model does not includepile-up and electronic readout noise, but these effects can be taken into account by futureextensions of the framework.The ideal-linear-observer performance derived here may not always mimic the perfor-mance of a human observer. However, if the measured count numbers are large enough tobe well described by Gaussian statistics, the optimal-linear-observer performance equals theperformance of the ideal observer and can therefore be expected to give an approximateupper limit to the detection performance achievable with advanced image processing.To further investigate the meaning of the NEQ and DQE matrices when a reconstruc-tion or image processing algorithm is applied to the raw data, we start by studying thecase where the post-processing algorithm is a linear, shift-invariant transformation. In thiscase, the Fourier transformed signal will be D postproc k ( u ) = P l U kl ( u ) D l ( u ) where U kl ( u )28re the components of a frequency-dependent transformation matrix. The combination ofthe detector and the algorithm can be seen as a composite detection system with a trans-fer function H sys k ( u , E ) = P l U kl ( u ) H k ( u , E ) and cross-spectral density W sys d + ( u ) fulfilling (cid:2) W sys d + (cid:3) kk ′ ( − u ) = P ll ′ U kl ( u ) [ W d + ] ll ′ ( − u ) U ∗ k ′ l ′ ( u ). (The latter identity follows from Eq.A2 and the corresponding formula for (cid:2) W sys d + (cid:3) kk ′ ( − u ) with d s replaced by F − U ( u ) F DS d s .)Now, Eqs. 8, 12 and 25 show that this leaves the frequency-dependent NEQ and DQE ma-trices as well as the frequency-dependent CRLB unchanged. These performance metrics arethus invariant under all linear shift-invariant processing algorithms and are therefore wellsuited for analyzing imaging performance. Also note that the transfer function H sys k ( u , E )and cross-spectral density W sys d + ( u ) of the composite system are useful quantities in theirown right for studying the effect of linear postprocessing algorithms, e.g. spectral distortioncorrection or deblurring algorithms.In contrast to linear algorithms, nonlinear image processing algorithms, such as manyiterative reconstruction and denoising algorithms, cannot be readily described by thepresent framework, although they may in some cases be linearized and thereby analyzed inan approximative sense. The NEQ and DQE metrics should thus be regarded as measures ofhow much information is available in the raw output data from the detector, i.e. they provideinformation about the amount of information available as input to the algorithm. Underthe reasonable assumption that providing more input information to most such algorithmswill also lead to better output images, the proposed metrics can be expected to be usefulpredictors of image quality also when nonlinear processing is applied. However, furtherinvestigations will be necessary in order to establish the validity of this assumption. VI. CONCLUSION
In this work, we have shown how the framework of linear-systems theory can be extendedto describe energy-resolving detectors. We have demonstrated a natural way of generalizingthe NEQ and DQE metrics to matrix-valued quantities containing information about boththe spatial resolution, detection efficiency and energy resolution of the detector. Further-more, we have demonstrated how basis material decomposition can be used to express thesematrices in a compact form and that they are closely related to a generalized version ofthe Cram´er-Rao lower bound which describes the detector performance for material decom-29osition. We have thus merged two approaches for detector performance assessment, thelinear-systems framework for describing detection task performance and the CRLB approachfor describing material decomposition performance.Although photon-counting detectors have been studied in the examples presented here,the proposed framework can also be extended to other spectral imaging systems, such asdual-layer detectors. In future work, we will also use the present framework to study morerealistic detector systems, where pileup and electronic noise are taken into account. Anothertopic for future work is to develop methods for measuring the matrix-valued NEQ and DQEexperimentally, for detector modules and for complete imaging systems. This new frame-work for detector performance characterization will help the development of photon-countingX-ray imaging systems by facilitating comparisons between different detector designs andelucidating the trade-off between different parameters, such as spatial resolution, energyresolution and dose efficiency.
ACKNOWLEDGMENTS
This study was supported by NIH Grant U01 EB01714003. The authors would like tothank Moa Yveborg and Jesse Tanguay for helpful discussions.
DISCLOSURE OF CONFLICTS OF INTEREST
Mats Persson is stockholder of and consultant for Prismatic Sensors AB. Norbert J. Pelcis consultant for Prismatic Sensors AB.
Appendix A: Fourier-domain covariance matrix
In this section we derive the Fourier-domain expression for the optimal-linear-observerdetectability. The matrix elements of the discrete-space Fourier operator and its inverse is( F DS ) n ( u ) = e − πi u · r n and (cid:0) F − (cid:1) n ( u ) = ∆ x ∆ y e πi u · r n . Eq. 6 then gives d ′ = (cid:16) F DS ∆d s (cid:17) † (cid:16) F † DS (cid:17) − Cov( d s ) − F − F DS ∆d s = (cid:0) ∆D s (cid:1) † (cid:16) F DS Cov( d s ) F † DS (cid:17) − ∆D s , (A1)30here † denotes conjugate-transpose. Since d s is wide-sense stationary, h F DS Cov( d s ) F † DS i k,k ′ ( u , u ′ ) = ∞ X n = −∞ ∞ X n ′ = −∞ Cov( d s n ,k , d s n ′ ,k ′ )e πi ( − u · r n + u ′ · r n ′ ) = ∞ X n = −∞ ∞ X ∆ n = −∞ K s ∆n ,k,k ′ e πi (( u ′ − u ) · r n + u ′ · ∆r n = ∞ X n = −∞ e πi ( u ′ − u ) · r n ∞ X ∆n = −∞ K s ∆n ,k,k ′ e πi u ′ · ∆r n = ∞ X n = −∞ e πi ( u ′ − u ) · r n W sk,k ′ ( − u ′ ) = 1∆ x ∆ y δ ( u − u ′ ) W sk,k ′ ( − u ′ ) . (A2)within the Nyquist region Nyq = n u : | u | < x , | v | < y o . Here, ∆n = n ′ − n and ∆r n = r n ′ − r n . Substituting into (A1) and exploiting the conjugate-symmetry of W sk,k ′ ( u )gives d ′ = (cid:0) ∆D s (cid:1) † (cid:16) F DS Cov( d s ) F † DS (cid:17) − ∆D s = ∆ x ∆ y Z Nyq X k X k ′ ∆ D sk ( u ) ∗ (cid:2) ( W s ) − (cid:3) ∗ k,k ′ ( u )∆ D sk ′ ( u )d u . (A3)where ∗ denotes complex conjugate.We will now express the noise correlations in terms of W d + ( u ), the cross-spectral den-sity of the sampled pulse train signal d + k ( r ) = P ∞ n = −∞ d s n ,k δ ( r − r n ). By observing that K s ∆n ,k,k ′ = ( K d ) k,k ′ ( ∆r n ), where K d is the autocovariance of the presampling signal d ,and using ( K d + ) k,k ′ ( ∆r ) = x ∆ y ( K d ) k,k ′ ( ∆r ) P n δ ( ∆r − ∆r n ) (Ref. 16, eq. 2.108), weobtain W sk,k ′ ( u ) = ∆ x ∆ y ( W d + ) k,k ′ ( u ). Assuming that aliasing is negligible, ∆ D sk ( u ) ≈ x ∆ y ∆ D k ( u ), and using Eq. 5 then allows us to express (A3) as d ′ = Z Nyq X k X k ′ Z E max Z E max H k ( u , E ) ∗ ∆ Q ( u , E ) ∗ · (cid:2) W − d + (cid:3) ∗ k,k ′ ( u ) H k ( u , E ′ )∆ Q ( u , E ′ )d E d E ′ d u . (A4)which is (7). Appendix B: Derivation of the CRLB
In this section we derive the form of the CRLB for ˜ A ( u ) in the Fourier domain. Toobtain the CRLB for the material decomposition ˜ A we introduce a change of variables by31efining ˜ A = (cid:16) ˜ A , ˜ A ∗ (cid:17) T = T θ . Here T is the block matrix I i II − i I , where I is the identitymatrix on the vector space on which Re ˜ A , Im ˜ A are defined. Using Eq. 24 and the identity T − = T † we get the CRLB for ˜ A asCov (cid:16) ˆ˜ A (cid:17) − = (cid:0) T Cov ( θ ) T † (cid:1) − = 14 T (cid:2) Cov ( θ ) − (cid:3) T † (B1) ≤ T ∂ d s ∂ θ ! T Cov( d s ) − ∂ d s ∂ θ T † = ∂ d s ∂ ˜ A , ∂ d s ∂ ˜ A ∗ ! † Cov( d s ) − ∂ d s ∂ ˜ A , ∂ d s ∂ ˜ A ∗ ! = F ˜ A , where F ˜ A is the Fisher matrix for ˜ A . Here the derivatives with respect to A and A ∗ aredefined as Wirtinger derivatives, which allows differentiating both real and complex functionswith respect to a complex variable (See Ref. 51, ch. A2). This gives F ˜ A = F ˜ A ˜ F ˜ A ˜ F ∗ ˜ A F ∗ ˜ A with F ˜ A = (cid:16) ∂ d s ∂ ˜ A (cid:17) † Cov( d s ) − (cid:16) ∂ d s ∂ ˜ A (cid:17) and ˜ F ˜ A = (cid:16) ∂ d s ∂ ˜ A (cid:17) † Cov( d s ) − (cid:16) ∂ d s ∂ ˜ A (cid:17) ∗ Since Eq. 15 gives ∂Q ( u ′ ,E ) ∂ ˜ A l ( u ) = ∂q ( E ) ∂A l δ ( u ′ − u ) in the linear approximation, we can calculatethe Jacobian of d s using the chain rule of differentiation: ∂d s n ,k ∂ ˜ A l ( u ) = Z R Z E max ∂d s n ,k ∂D k ( u ′ ) ∂D k ( u ′ ) ∂Q ( u ′ , E ) ∂Q ( u ′ , E ) ∂ ˜ A l ( u ) d E d u ′ (B2)= ∆ x ∆ y Z E max e πi u · r n H k ( u , E ) ∂q ( E ) ∂A l d E = ∆ x ∆ y q tot ∂q tot ∂A l H B k,l ( u ) e πi u · r n . where the delta function has allowed us to eliminate the integral over u ′ . Note that theconjugate symmetry of H B k,l ( u ) implies that ∂d s n ,k ∂ ˜ A l ( u ) ∗ = (cid:16) ∂d s n ,k ∂ ˜ A l ( u ) (cid:17) ∗ = ∆ x ∆ y q tot ∂q tot ∂A l H B k,l ( − u ) e − πi u · r n .Inserted into (24) this gives the elements of F ˜ A as( F ˜ A ) l,l ′ ( u , u ′ ) = ∞ X n = −∞ ∞ X n ′ = −∞ X k X k ′ ∂d s n ,k ∂ ˜ A l ( u ) ! ∗ (cid:2) Cov( d s ) − (cid:3) n , n ′ ,k,k ′ ∂d s n ′ ,k ′ ∂ ˜ A l ′ ( u ′ )(B3)= ∆ x ∆ y q tot2 ∂q tot ∂A l ∂q tot ∂A l ′ X k X k ′ H B k,l ( u ) ∗ ∞ X n = −∞ ∞ X n ′ = −∞ e − πi u · r n (cid:2) Cov( d s ) − (cid:3) n , n ′ ,k,k ′ e πi u ′ · r n ′ H B k ′ ,l ′ ( u ′ )= ∆ x ∆ y q tot2 ∂q tot ∂A l ∂q tot ∂A l ′ h H B† F DS Cov( d s ) − F † DS H B i l,l ′ ( u , u ′ )= 1 q tot2 ∂q tot ∂A l ∂q tot ∂A l ′ (cid:20) H B† (cid:16) F DS Cov( d s ) F † DS (cid:17) − H B (cid:21) l,l ′ ( u , u ′ ) , where H B is the transfer matrix matrix with elements H B k,l ( u ). Similarly, the values of ˜ F ˜ A (cid:16) ˜ F ˜ A (cid:17) l,l ′ ( u , u ′ ) = 1 q tot2 ∂q tot ∂A l ∂q tot ∂A l ′ (cid:20) H B† (cid:16) F DS Cov( d s ) F † DS (cid:17) − H B (cid:21) l,l ′ ( u , − u ′ ) . (B4)Using (A2), we obtain the Fisher matrix elements as( F ˜ A ) l,l ′ ( u , u ′ ) = 1 q tot2 ∂q tot ∂A l ∂q tot ∂A l ′ X k X k ′ H B k,l ( u ) ∗ δ ( u − u ′ ) (cid:2) W − d + (cid:3) k,k ′ ( − u ) H B k ′ ,l ′ ( u ) (B5)= 1 q tot2 ∂q tot ∂A l ∂q tot ∂A l ′ δ ( u − u ′ )NEQ B l,l ′ ( u ) . A similar calculation shows that (cid:16) ˜ F ˜ A (cid:17) l,l ′ ( u , u ′ ) is proportional to δ ( u + u ′ ). Since we re-strict ourselves to estimating the Fourier coefficient in the right half-plane { u : u > u = 0 , v > } , (cid:16) ˜ F ˜ A (cid:17) l,l ′ ( u , u ′ ) = , which givesCov (cid:16) ˆ˜ A (cid:17) − ≤ F ˜ A = F ˜ A
00 F ∗ ˜ A , (B6)with components of F ˜ A given by (B5).The case u = requires special treatment since the fact that Im ˜ A ( ) = makes F ˜ A singular if u = is included. However, (B5) is sufficient for our purposes since Cov (cid:16) ˆ˜ A l ( u ) (cid:17) can be extrapolated to u = as long as the NEQ is a smooth function of u near the origin.Since different spatial frequencies are independent, the Fisher information matrix F ˜ A isblock diagonal and the block corresponding to spatial frequency u is given by a rescaledversion of NEQ B ( u ), where B is the set of basis functions in which ˜ A is expressed. F ˜ A canthus be inverted separately for each u . This allows the CRLB (B6) to be expressed asCov (cid:16) ˆ˜ A l ( u ) , ˆ˜ A l ′ ( u ′ ) (cid:17) ≥ q tot2 (cid:2) NEQ B ( u ) − (cid:3) l,l ′ ∂q tot ∂A l ∂q tot ∂A l ′ δ ( u − u ′ ) . (B7)which is (25). Appendix C: Detectability in a combination of basis images
In this section we derive the upper limit of detectability in a frequency-dependentweighted linear combination of basis images. Combining the CRLB for θ , Cov( ˆ θ ) ≥ F θ − ,33ith the formula for the ideal linear observer (Sec. 13.2.12 of Ref. 37), d ′ ≤ ∆ θ T F θ ∆ θ = 14 ∆ θ T T † TF θ T † T ∆ θ = ∆ ˜ A † F ˜ A ∆ ˜ A = (C1) ∆ ˜ A † F ˜ A ∆ ˜ A + ∆ ˜ A T F ∗ ˜ A ∆ ˜ A ∗ . The block-diagonality of the Fisher information matrix allows us to treat different spatialfrequencies independently, and therefore the first term in this sum is an integral over thehalf-plane u >
0. Since ∆ ˜ A l ( u ) and NEQ B ( u ) are both conjugate-symmetric, including thesecond term is equivalent to extending the integral to the entire u plane: d ′ ≤ q tot2 Z R X l X l ′ ∂q tot ∂A l ∆ ˜ A l ( u ) † (cid:2) NEQ B ( u ) (cid:3) l,l ′ ∆ ˜ A l ′ ( u ) ∂q tot ∂A l ′ d u (C2)= Z R ∆S B ( u ) † NEQ B ( u ) ∆S B ( u )d u . which is (27). REFERENCES a) Author to whom correspondence should be addressed. Electronic mail:[email protected] K. Taguchi and J. S. Iwanczyk, “Vision 20/20: Single photon counting x-ray detectors inmedical imaging,” Med. Phys. , 100901 (2013). P. M. Shikhaliev, T. Xu, and S. Molloi, “Photon counting computed tomography: Conceptand initial results,” Medical Physics , 427–436 (2005). H. Bornefalk and M. Danielsson, “Photon-counting spectral computed tomography usingsilicon strip detectors: a feasibility study,” Phys. Med. Biol. , 1999–2022 (2010). P. M. Shikhaliev, “Energy-resolved computed tomography: first experimental results,”Phys. Med. Biol. , 5595–5613 (2008). J. Giersch, D. Niederl¨ohner, and G. Anton, “The influence of energy weighting on x-rayimaging quality,” Nucl. Instrum. Meth. A , 68 – 74 (2004), proceedings of the 5thInternational Workshop on Radiation Imaging Detectors. J. Berglund, H. Johansson, M. Lundqvist, B. Cederstr¨om, and E. Fredenberg, “En-ergy weighting improves dose efficiency in clinical practice: implementation on a spectralphoton-counting mammography system,” Journal of Medical Imaging , 031003 (2014).34 J. P. Schlomka, E. Roessl, R. Dorscheid, S. Dill, G. Martens, T. Istel, C. Bumer, C. Her-rmann, R. Steadman, G. Zeitler, A. Livne, and R. Proksa, “Experimental feasibilityof multi-energy photon-counting K-edge imaging in pre-clinical computed tomography,”Phys. Med. Biol. , 4031–4047 (2008). M. Persson, B. Huber, S. Karlsson, X. Liu, H. Chen, C. Xu, M. Yveborg, H. Bornefalk,and M. Danielsson, “Energy-resolved CT imaging with a photon-counting silicon-stripdetector,” Physics in Medicine and Biology , 6709 (2014). Z. Yu, S. Leng, S. M. Jorgensen, Z. Li, R. Gutjahr, B. Chen, A. F. Halaweish, S. Kap-pler, L. Yu, E. L. Ritman, and C. H. McCollough, “Evaluation of conventional imagingperformance in a research whole-body CT system with a photon-counting detector array,”Phys. Med. Biol. , 1572 (2016). J. P. Ronaldson, R. Zainon, N. J. A. Scott, S. P. Gieseg, A. P. Butler, P. H. Butler, andN. G. Anderson, “Toward quantifying the composition of soft tissues by spectral CT withMedipix3,” Med. Phys. , 6847–6857 (2012). R. N. Cahn, B. Cederstr¨om, M. Danielsson, A. Hall, M. Lundqvist, and D. Nygren,“Detective quantum efficiency dependence on x-ray energy weighting in mammography,”Medical Physics , 2680–2683 (1999). T. G. Schmidt, “Optimal “image-based” weighting for energy-resolved CT,” Med. Phys. , 3018–3027 (2009). R. E. Alvarez and A. Macovski, “Energy-selective reconstructions in x-ray computerisedtomography,” Phys. Med. Biol. , 733–744 (1976). E. Roessl and R. Proksa, “K-edge imaging in x-ray computed tomography using multi-binphoton counting detectors,” Phys. Med. Biol. , 4679 (2007). P. Sharp, ed.,
ICRU Report 54: Medical Imaging - the Assessment of Image Quality (Inter-national Commision on Radiation Units and Measurements, Bethesda, MD, USA, 1996). I. A. Cunningham, “Applied linear-systems theory,” (SPIE Press, Bellingham, Washing-ton, USA, 2000) Chap. 2, pp. 79–159. R. J. Acciavatti and A. D. A. Maidment, “A comparative analysis of OTF,NPS, and DQE in energy integrating and photon counting digital x-ray detectors,”Medical Physics , 6480–6495 (2010). J. Tanguay, S. Yun, H. K. Kim, and I. A. Cunningham, “The detective quan-tum efficiency of photon-counting x-ray detectors using cascaded-systems analyses,”35edical Physics , 041913 (2013), 041913. J. Xu, W. Zbijewski, G. Gang, J. W. Stayman, K. Taguchi, M. Lundqvist, E. Fredenberg,J. A. Carrino, and J. H. Siewerdsen, “Cascaded systems analysis of photon countingdetectors,” Medical Physics , 101907 (2014), 101907. K. Stierstorfer, “Modeling the frequency-dependent detective quantum efficiency ofphoton-counting x-ray detectors,” Medical Physics , In press. H. Chen, C. Xu, M. Persson, and M. Danielsson, “Optimization of beam quality forphoton-counting spectral computed tomography in head imaging: simulation study,”Journal of Medical Imaging , 043504 (2015). J. Cammin, S. Kappler, T. Weidinger, and K. Taguchi, “Evaluation of mod-els of spectral distortions in photon-counting detectors for computed tomography,”Journal of Medical Imaging , 023503 (2016). K. Taguchi, M. Zhang, E. C. Frey, X. Wang, J. S. Iwanczyk, E. Nygard, N. E.Hartsough, B. M. W. Tsui, and W. C. Barber, “Modeling the performance of aphoton counting x-ray detector for CT: Energy response and pulse pileup effects,”Medical Physics , 1089–1102 (2011). J. Tanguay, S. Yun, H. K. Kim, and I. A. Cunningham, “Detective quantum efficiency ofphoton-counting x-ray detectors,” Medical Physics , 491–509 (2015). S. Faby, J. Maier, S. Sawall, D. Simons, H.-P. Schlemmer, M. Lell, and M. Kachelrieß,“An efficient computational approach to model statistical correlations in photon countingx-ray detectors,” Medical Physics , 3945–3960 (2016). K. Taguchi, C. Polster, O. Lee, K. Stierstorfer, and S. Kappler, “Spatio-energetic crosstalk in photon counting detectors: Detector model and correlated poisson data generator,”Medical Physics , 6386–6404 (2016). P. L. Rajbhandary, S. S. Hsieh, and N. J. Pelc, “Effect of spatio-energycorrelation in PCD due to charge sharing, scatter, and secondary photons,”Proc. SPIE , 101320V–101320V–8 (2017). P. M. Shikhaliev, S. G. Fritz, and J. W. Chapman, “Photon counting multi-energy x-ray imaging: Effect of the characteristic x rays on detector performance,”Medical Physics , 5107–5119 (2009). S. Richard and J. H. Siewerdsen, “Cascaded systems analysis of noise reduction algorithmsin dual-energy imaging,” Medical Physics , 586–601 (2008).36 E. Fredenberg, M. Hemmendorff, B. Cederstr¨om, M. ˚Aslund, and M. Daniels-son, “Contrast-enhanced spectral mammography with a photon-counting detector,”Medical Physics , 2017–2029 (2010). M. Yveborg, M. Danielsson, and H. Bornefalk, “Theoretical comparison of a dual energysystem and photon counting silicon detector used for material quantification in spectralCT,” IEEE Transactions on Medical Imaging , 796–806 (2015). H. Chen, M. Danielsson, and C. Xu, “Size-dependent scanning parameters (kVp andmAs) for photon-counting spectral CT system in pediatric imaging: simulation study,”Physics in Medicine and Biology , 4105 (2016). M. J. Tapiovaara and R. Wagner, “SNR and DQE analysis of broad spectrum x-ray imag-ing,” Physics in Medicine and Biology , 519 (1985). H. Bornefalk, “Task-based weights for photon counting spectral x-ray imaging,”Medical Physics , 6065–6073 (2011). M. Yveborg, M. Persson, and H. Bornefalk, “Optimalfrequency-based weighting for spectral x-ray projection imaging,”IEEE Transactions on Medical Imaging , 779–787 (2015). P. L. Rajbhandary, M. Persson, and N. J. Pelc, “Frequency dependent dqe of photoncounting detector with spectral degradation and cross-talk,” Proc.SPIE , 10573 –10573 – 10 (2018). H. H. Barrett and K. J. Myers,
Foundations of Image Science (John Wiley & Sons, Inc.,Hoboken, New Jersey, 2004). E. Roessl and C. Herrmann, “Cram´er-Rao lower bound of basis image noise in multiple-energy x-ray imaging,” Physics in Medicine & Biology , 1307 (2009). E. Roessl, B. Brendel, K. J. Engel, J. P. Schlomka, A. Thran, and R. Proksa,“Sensitivity of photon-counting based k-edge imaging in x-ray computed tomography,”IEEE Transactions on Medical Imaging , 1678–1690 (2011). S. M. Kay,
Fundamentals of Statistical Signal Processing: Estimation Theory (PrenticeHall PTR, Upper Saddle River, New Jersey, 1993). R. E. Alvarez, “Near optimal energy selective x-ray imaging system performance withsimple detectors,” Medical Physics , 822–841 (2010). R. Alvarez and E. Seppi, “A comparison of noise and dosein conventional and energy selective computed tomography,”37EEE Transactions on Nuclear Science , 2853–2856 (1979). 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,”Medical Physics , 659–667 (1981). F. S. Philippe T. Pinard, Hendrix Demers and R. Gauvin, “pyPENELOPE software pack-age,” pypenelope.sourceforge.net/. A. M. Hernandez and J. M. Boone, “Tungsten anode spectral model using in-terpolating cubic splines: Unfiltered x-ray spectra from 20 kV to 640 kV,”Medical Physics , 042101 (2014). J. Punnoose, J. Xu, A. Sisniega, W. Zbijewski, and J. H. Siewerdsen, “Technicalnote: spektr 3.0-A computational tool for x-ray spectrum modeling and analysis,”Medical Physics , 4711–4717 (2016). M. J. Berger, J. H. Hubbell, S. M. Seltzer, J. Chang, J. S. Coursey, R. Suku-mar, D. S. Zucker, and K. Olsen, “XCOM: Photon Cross Section Database,”http://physics.nist.gov/xcom. National Institute of Standards and Technology, Gaithers-burg, MD (2005). J.-B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A three-dimensional statisticalapproach to improved image quality for multislice helical CT,” Medical Physics , 4526–4544. K. M. Brown, S. Zabi´c, and G. Shechter, “Impact of spectral separation in dual-energyCT with anti-correlated statistical reconstruction,” Presented at The thirteenth Inter-national Meeting on Fully Three-Dimensional Image Reconstruction in Radiology andNuclear Medicine, Newport, RI,USA (2015). P. J. Schreier and L. L. Scharf,
Statistical Signal Processing of Complex-Valued Data: The Theory of Improper and Noncircular Signals (Cambridge University Press, 2010). 38 able
I: List of symbols and their dimensions. L: length; E:energy.Symbol Dimension Description ∗ Complex conjugate † Conjugate-transpose A L Basis coefficient line integral˜ A L Fourier transform of Ad D L Fourier transform of dd s D s d s d + L − Sampled signal as delta-pulse train d ′ DQE E − Detective quantum efficiency matrix in energy basis
DQE B x , ∆ y L Detector pixel width and height E E Incident energy of photon ε E Registered energy in event F ˜ A L Fisher information matrix for ˜ A F Continuous Fourier transformation F DS Discrete-space Fourier transformation f L − Material basis function G h L − Point-spread function H H B L − Rescaled transfer matrix in material basis K d dK d + L − Cross-covariance matrix of d + K s d s L TF pre n = ( n x , n y ) 1 Discrete pixel coordinate NEQ L − E − Noise-equivalent quanta matrix in energy basis
NEQ B L − Noise-equivalent quanta matrix in material basisNyq Nyquist region of the (u,v) planeΦ L − E − Prepatient photons per area and energy q L − E − Incident photons per area and energy q tot L − Total incident photons per area Q E − Fourier transform of q r = ( x, y ) L Position on detector r n L Pixel center position ∆S L Task function in energy basis ∆S B L Task function in material basis T X to ( X , X ∗ ) T u = ( u, v ) L − Spatial frequency W s d s W d + L − Cross-spectral density of d + X Expected value of X
20 40 60 80 100 12000.20.40.60.81 E (keV) R e l a t i v e s pe c t r a l den s i t y a IncidentRegistered E (keV) E ’ ( k e V ) b D Q E ⋅ ∆ E D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task −1−0.500.51 R e l a t i v e t a sk f un c t i on DQE ∆ S/ ∆ S R e l a t i v e s pe c t r a l den s i t y d E (keV) E ’ ( k e V ) e D Q E ⋅ ∆ E D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task R e l a t i v e t a sk f un c t i on D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task R e l a t i v e t a sk f un c t i on D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task R e l a t i v e t a sk f un c t i on D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task R e l a t i v e t a sk f un c t i on Fig. E = 1 keV. (b) DQE matrix in the nonoverlapping spectrumcase. (c) diagonal of the
DQE matrix together with the task function and DQE task for a densityimaging task in the nonoverlapping spectrum case. (d) Incident and deposited spectra in theoverlapping spectrum case. (e)
DQE matrix in the overlapping spectrum case. (f-i) Diagonal ofthe
DQE matrix in the overlapping spectrum case, together with task function and DQE task forfour different tasks: (f) nonoverlapping density imaging task; (g) overlapping density imaging task;(h) nonoverlapping spectral imaging task; and (i) overlapping spectral imaging task. Since
DQE issingular, the plotted DQE curves in (c,f-i) show the coefficient in front of δ ( E − E ′ ).
20 40 60 80 100 1200123 E (keV) R e l a t i v e s pe c t r a l den s i t y a IncidentRegistered0 10 20050 E (keV) E ’ ( k e V ) b D Q E ⋅ ∆ E D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task −1−0.500.51 R e l a t i v e t a sk f un c t i on DQE ∆ S/ ∆ S D Q E ( de l t a f un c t i on c oe ff. ) DQE
Task R e l a t i v e t a sk f un c t i on Fig.
DQE matrix, discretized with ∆ E = 1 keV. Note that the colorscale is logarithmic, to visualize the low-intensity background of 6 . · − keV − . (c) diagonal of the DQE matrix together with the task function and DQE task for a density imaging task. (d) diagonalof the
DQE matrix together with the task function and DQE task for a spectral imaging task. Since
DQE is singular, the plotted DQE curves in (c-d) show the coefficient in front of δ ( E − E ′ ). PS F ( c oun t s / p i x e l )
30 keV
Bin 1 Bin 2 Bin 3 Bin 4 Bin 50.2 0.25 0.3 x (mm) PS F ( c oun t s / p i x e l ) -0.5 0 0.5x (mm)00.20.40.60.81
70 keV x (mm)
100 keV x (mm)
Fig. x axis of the simulated energy-bin point-spread functions in the CdTemodel measured as counts per pixel, ∆ x ∆ y h k ( ∆r , E ), for three incident monoenergies: 30, 70 and100 keV. The lower row contains close-ups of the region closest to the pixel border. (In the lowerrow, the data points are not connected by lines since the psf is not expected to be smooth.) T r an s f e r f un c t i on
30 keV
Bin 1 Bin 2 Bin 3 Bin 4 Bin 50 1 2 3 u (mm -1 ) -0.0500.050.10.150.2 T r an s f e r f un c t i on
70 keV u (mm -1 ) -0.0500.050.10.150.2 0 1 2 3-0.200.20.40.60.81
100 keV u (mm -1 ) -0.0500.050.10.150.2 Fig. H k ( u , E ) in the CdTe model for three incident monoenergies: 30, 70and 100 keV. These are plotted along the positive u axis up to three times the Nyquist frequency1 mm − . The lower row contains close-up views of the upper plots. C r o ss - s pe c t r a l den s i t y ( mm ) Bin 1
Bin 1 Bin 2 Bin 3 Bin 4 Bin 50 0.5 1 u (mm -1 ) C r o ss - s pe c t r a l den s i t y ( mm ) Bin 2 u (mm -1 ) Bin 3 u (mm -1 ) Bin 4 u (mm -1 ) Bin 5 u (mm -1 ) Fig. x ∆ y ( W d + ) kk ′ ( u ) in the CdTe model for broad-spectrum illumi-nation, along the positive u axis up to the Nyquist frequency (1 mm − ). Each column of plotsrepresents one k and each curve one k ′ . The lower row contains close-up views of the upper plots.Note that each cross-term is plotted twice since ( W d + ) kk ′ ( u ) is symmetric in k and k ′ . (keV) E ’ ( k e V ) a −2 keV −2 )0200400600 0 50 10005001000150020002500 E (keV) N E Q ( mm − k e V − ) bE (keV) E ’ ( k e V ) c −1 )00.020.040.060.08 0 50 10000.050.10.15 E (keV) D Q E ( k e V − ) d Fig.
NEQ matrix for the CdTe model. (b) Diagonal elements of the
NEQ matrix in (a). (c) Zero-frequency
NEQ matrix for the CdTe model. (d) Diagonal elements of the
DQE matrix in (c). u (lp/mm) N E Q B ( mm − ) a WaterBoneCross−term 0 1 2 300.20.40.60.81 u (lp/mm) D Q E B b WaterBoneCross−term0 1 2 300.20.40.60.81 u (lp/mm) D Q E t a sk f o r de t e c t i on c Max DQE task
Min DQE task D Q E T a sk f o r quan t i f i c a t i on d WaterBone
Fig.
NEQ B matrix with water and bone as basis functions. The plotswith markers apply to the CdTe simulation model whereas the plots without markers are thecorresponding values for an ideal detector. (b) Components of the DQE B matrix for the CdTesimulation model with water and bone as basis functions. (c) Maximum and minimum DQE task for any detection task, for the CdTe model. (d) DQE task for quantifying the amount of water andbone, respectively, in a two-basis decomposition with the CdTe model.for quantifying the amount of water andbone, respectively, in a two-basis decomposition with the CdTe model.