A framework for performance characterization of energy-resolving photon-counting detectors

A framework for performance characterization of energy-resolving photon-counting detectors

Mats Persson matspers@stanford.edu Departments of Bioengineering and Radiology, Stanford University, Stanford 94305, California, USA
(Corresponding author. E-mail: matspers@stanford.edu)
   Paurakh L. Rajbhandary Departments of Electrical Engineering and Radiology, Stanford University, Stanford 94305, California, USA    Norbert J. Pelc Departments of Electrical Engineering, Bioengineering and Radiology, Stanford University, Stanford 94305, California, USA
July 14, 2019
Abstract

Purpose: Photon-counting energy resolving detectors are subject to intense research interest, and there is a need for a general framework for performance assessment of these detectors. The commonly used linear-systems theory framework, which measures detector performance in terms of noise-equivalent quanta (NEQ) and detective quantum efficiency (DQE) is widely used for characterizing conventional X-ray detectors but does not take energy-resolving capabilities into account. The purpose of this work is to extend this framework to encompass energy-resolving photon-counting detectors and elucidate how the imperfect energy response of real-world detectors affects imaging performance.
Method: We generalize NEQ and DQE to matrix-valued quantities as functions of spatial frequency, and show how these matrices can be calculated from simple Monte Carlo simulations. To demonstrate how the new metrics can be interpreted, we compute them for simplified models of fluorescence and Compton scatter in a photon-counting detector and for a Monte Carlo model of a CdTe detector with pixels.
Results: Our results show that the ideal-linear-observer performance for any detection or material quantification task can be calculated from the proposed generalized NEQ and DQE metrics. Off-diagonal elements in these matrices are shown to be related to loss of energy information due to imperfect energy resolution. The Monte Carlo model of the CdTe detector predicts a zero-frequency dose efficiency relative to an ideal detector of 0.86 and 0.65 for detecting water and bone, respectively. When the task instead is to quantify these materials, the corresponding values are 0.34 for water and 0.26 for bone.
Conclusions: We have developed a framework for assessing the performance of photon-counting energy-resolving detectors and shown that the matrix-valued NEQ and DQE metrics contain sufficient information for calculating the dose efficiency for both detection or quantification tasks, the task having any spatial and energy dependence. This framework will be beneficial for the development of photon-counting X-ray detectors.
Key words: X-ray imaging, photon counting, spectral imaging, linear-systems theory, detective quantum efficiency, noise-equivalent quanta

pacs:

I Introduction

The development of energy-resolving photon-counting detectors for medical X-ray imaging is an active research field.Taguchi and Iwanczyk (2013); Shikhaliev, Xu, and Molloi (2005); Bornefalk and Danielsson (2010); Shikhaliev (2008); Giersch, Niederlöhner, and Anton (2004) Systems with such detectors are currently commercially available for mammographyBerglund et al. (2014), and prototypes have been demonstrated for computed tomography (CT)Schlomka et al. (2008); Persson et al. (2014); Yu et al. (2016); Ronaldson et al. (2012). An energy-resolving photon-counting detector uses several electronic threshold levels to separate the registered counts into energy bins depending on the deposited energy in each event. This enables spectroscopic imaging, i.e. using the information contained in the energy distribution of the incoming spectrum to obtain more information about the imaged object. The energy information can be used for improving the the signal-difference-to-noise ratio (SDNR) by optimal weightingCahn et al. (1999); Schmidt (2009), for removing beam-hardening artifacts and for generating material-selective images.Alvarez and Macovski (1976); Roessl and Proksa (2007); Schlomka et al. (2008)

The ongoing development of improved energy-resolving photon-counting X-ray detectors raises the question of how detector performance should be measured. Using a relevant metric is important in order to optimize and compare detector designs with different detector materials and pixel sizes or with different properties of the readout electronics such as the energy 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).Sharp (1996) This metric, which is usually expressed as a function of spatial frequency, measures how well the detector is able to use the incoming photon statistics for each spatial frequency compared to an ideal detector. A closely related metric is the noise-equivalent quanta (NEQ), which measures the number of photons that an ideal detector would require in order to obtain the same signal-difference-to-noise ratio as the studied detector. An overview of the theory of DQE and NEQ for ordinary (non-spectral) detectors can be found in Ref. 16.

Several studies have been published previously on how to model the frequency-dependent detective quantum efficiency of photon-counting detectors,Acciavatti and Maidment (2010); Tanguay et al. (2013); Xu et al. (2014); Stierstorfer () but these studies do not take energy information into account. The impact of detector nonidealities on spectral imaging tasks has been studied by several authorsChen et al. (2015); Cammin et al. (2016); Taguchi et al. (2011); Tanguay et al. (2015); Faby et al. (2016); Taguchi et al. (2016); Rajbhandary, Hsieh, and Pelc (2017). However, these publications only study the zero-spatial-frequency performance and do not investigate the frequency-dependent spectral performance. Shikhaliev et al. Shikhaliev, Fritz, and Chapman (2009) studied the effect of characteristic X-rays on the spectrum shape and on the spatial resolution. Richard and SiewerdsenRichard and Siewerdsen (2008) used linear-systems theory to study the performance of a dual-energy system. Also, the detectability for specific tasks with both spectral and spatial dependence has been studied by Fredenberg et al.Fredenberg et al. (2010), Yveborg et al. Yveborg, Danielsson, and Bornefalk (2015) and Chen et al.Chen, Danielsson, and Xu (2016)

A common way to combine information from several spectral channels is to make a weighted sum of the energy-selective images, with weights chosen to be optimal for the considered imaging task.Tapiovaara and Wagner (1985); Cahn et al. (1999); Schmidt (2009) A good measure of the performance of an energy-resolving detector should therefore take optimal weighting into account. Taking the spatial-frequency dependence of the signal and noise into account when calculating optimal weights can improve detectability in an energy-weighted image, either by calculating an optimized weight for each imageBornefalk (2011) or by letting the weights themselves be functions of spatial frequencyYveborg, Persson, and Bornefalk (2015). This is equivalent to applying different spatial filters to the different energy bin images before forming the weighted sum, an approach studied in the dual-energy case by Richard and Siewerdsen.Richard and Siewerdsen (2008)

Another way to utilize spectral X-ray measurement is to perform basis material decomposition, which builds on the observation that the linear attenuation coefficient of any substance in the human body can be approximated well as a linear combination of a small number of basis functions, typically two, or three in the presence of a high-atomic-number contrast agent with a K-edge in the diagnostic energy range. Basis material decomposition then amounts to using the energy-resolved X-ray data to estimate the amount of each of these basis materials along each projection line, thereby completely characterizing the energy-dependent attenuation of the object.Alvarez and Macovski (1976); Roessl and Proksa (2007) This is useful in particular for CT, to remove beam-hardening artifacts and characterize object composition.

The purpose of this work is to present a general framework for comparing the performance of energy-resolving imaging detectors, encompassing both energy-weighting and material-decomposition performance. To this end, we derive an expression for the detectability achieved by an ideal linear observer, taking into account the full spatial-frequency dependence of the available energy information. We also propose generalized NEQ and DQE metrics for measuring energy-resolving detector performance. To show how the proposed metrics can be used in practice, we apply the developed framework to simple models of detectors with nonideal energy response and to a more realistic simulation model of a photon-counting detector.

Ii Theory

ii.1 Matrix-valued NEQ and DQE

Symbol Dimension Description
Complex conjugate
Conjugate-transpose
Basis coefficient line integral
Fourier transform of
Presampling detected signal
Fourier transform of
Sampled signal as discrete sequence
Discrete-space Fourier transform of
Sampled signal as delta-pulse train
Optimal-linear-observer detectability
Detective quantum efficiency matrix in energy basis
Detective quantum efficiency matrix in material basis
Detector pixel width and height
Incident energy of photon
Registered energy in event
Fisher information matrix for
Continuous Fourier transformation
Discrete-space Fourier transformation
Material basis function
Large-area gain
Point-spread function
Transfer function
Rescaled transfer matrix in material basis
Cross-covariance matrix of
Cross-covariance matrix of
Cross-covariance matrix of
Basis transformation matrix element
Presampling modulation transfer function
1 Discrete pixel coordinate
Noise-equivalent quanta matrix in energy basis
Noise-equivalent quanta matrix in material basis
Nyquist region of the (u,v) plane
Prepatient photons per area and energy
Incident photons per area and energy
Total incident photons per area
Fourier transform of
Position on detector
Pixel center position
Task function in energy basis
Task function in material basis
Matrix mapping to
Spatial frequency
Cross-spectral density of
Cross-spectral density of
Expected value of
Table 1: List of symbols and their dimensions. : length; : energy.

In the following, we will use parentheses to denote that a quantity is a function of a continuous variable and use subscripts when the variable is discrete. Most of our notation is similar to Ref. 16 and 18. Consider a situation where a distribution of X-ray photons is incident on a detector. Let be a random process describing the number of incident photons per area and energy as a function of position and energy . This photon distribution is measured by an energy-resolving detector and registered as counts in energy bins where is the discrete coordinate of a detector pixel in the plane and is the index of the energy bin. To avoid boundary effects, we assume that the detector and the photon distribution are of infinite extent. Both and are random processes, and their expected values will be denoted by and , respectively. It is also useful to define a presampling detected signal , whose value at every point is the number of counts that would be registered in each energy bin by a fictitious pixel centered at . Therefore, where is the center of pixel and and are the pixel center-to-center distances in the and directions, respectively.

The continuous Fourier transform of with respect to is a random process given by , and is defined analogously as the transform of . The 2D discrete-space Fourier transform of with respect to is given by . We now assume that the system is linear and shift-invariant, which allows us to introduce the system point-spread function for the energy bin that relates the expected values of the the incoming photon distribution and the registered count distribution:

(1)

where is the maximum incident energy. Here, the pixel cell area has been broken out of the definition of , meaning that can be interpreted as a presampling distribution of registered quanta per unit pixel area in the energy bin when the input signal is a beam incident on a single point of the detector. thus contains information about both the quantum detection efficiency and the spatial resolution of the detector. Denoting the Fourier transform of with respect to and by ,

(2)

In analogy with conventional X-ray detectors, it is possible to define a presampling modulation transfer function for each incident energy and for each energy bin.

To calculate the performance of the detector for a given imaging task, we also need to know the covariance between different measurements . Assuming that the noise is wide-sense stationary, this is given by the cross-covariance matrix with elements

(3)

The cross-spectral density of can then be calculated as

(4)

where . For , is the noise power spectrum (NPS) in energy bin , whereas for describes the frequency dependence of correlations between different energy bins.

Now assume that the studied detector is used to discriminate between two cases: a object-absent (background) case with incident photon distribution and an object-present case where the incident photon distribution is , i.e. a signal-known-exactly/background-known-exactly (SKE/BKE) task. The difference in the expected presampled signal between the two cases is with the Fourier transform

(5)

where is the expected signal difference in the Fourier domain. The expected difference in the sampled signal is with the discrete-space Fourier transform .

Our measure of the detector performance for this task will be squared signal-difference-to-noise ratio (“detectability”) of the optimal linear observer, which tells us how well this model observer can discriminate between the two cases (object present or absent) given the available measurements in all pixels and all energy bins. Viewing and for all and as elements of vectors and (see Sec. 13.2.12 of Ref. 36), can be expressed as

(6)

By changing basis using the discrete-space Fourier transform , this expression is transformed into (see appendix A)

(7)

where denotes complex conjugate and denotes the Nyquist region .

The ideal-linear-observer performance is thus obtained as an integral over spatial frequencies of a quadratic form in , the difference vector of registered counts in all energy bins for all spatial frequencies . To attain this performance limit, the model observer needs to take data from all the energy bins into account, 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 which attains the same performance for that task.Yveborg, Persson, and Bornefalk (2015)

We will now express (7) in a way that separates parameters specific to the detector from parameters specific to the discrimination task. To do this, we need to assume that signal aliasing is negligible, so that . This is true if the pixel size is sufficiently small compared to other resolution-limiting factors, such as the detector point-spread function or the focal spot size, 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 will express the noise correlations in terms of , the cross-spectral density of the sampled pulse train signal . By observing that , where is the autocovariance of the presampling signal , and using (Ref. 16, eq. 2.108), we obtain . This allows (7) to be expressed as

(8)

To rewrite (8) in a form that is more similar to the corresponding expression for a conventional detector, we introduce the relative signal difference and define a frequency-dependent matrix with elements given by

(9)

This gives

(10)

which is analogous to the formula for conventional X-ray detectors (Ref. 16, Eq. 2.150, for non-pixelated detectors):

(11)

where (Ref. 16, Eqs. 2.193 and 2.209)

(12)

In this equation, is the large-area gain (average registered counts per pixel divided by ), is the presampling modulation transfer function and is a transfer function defined in analogy with .

is thus a natural generalization of the noise-equivalent quanta used to describe conventional detectors. It is a frequency-dependent matrix with continuous indices and encoding all the information about the detector that is relevant for evaluating the detection performance. It depends on the background case noise characteristics, and thus on the incident spectrum in that setting, but not on the discrimination task. The discrimination task function determines how the different frequency components of are weighted together to give the detectability. The matrix therefore defines a quadratic form in giving the contribution to at each spatial frequency.

The above description applies to a continuous-to-discrete imaging system, which is the relevant type of system for multibin photon-counting detectors. Note, however, that the analysis, apart from the sampling step, can be adapted to a continuous-to-continuous imaging system by replacing , and by continuous variables. An ideal energy-resolving imaging system is a continuous-to-continuous imaging system that registers the correct position and energy for each incoming photon. Its transfer function is therefore , where denotes the registered energy, and its noise cross-spectrum is white, i.e. is constant with respect to . For an ideal photon-counting detector, the input signal is Poisson distributed and the registered counts at different energies are uncorrelated, so that where is the expected number of incident photons per energy and area (Ref. 16, Sec. 2.6.2.3). For a pixelated but otherwise ideal detector, tends to as and tend to 0, which gives for an ideal detector.

In analogy with the conventional definition of DQE, , we can now define a matrix by normalizing by the ideal-detector performance. Therefore, let

(13)

or equivalently,

(14)

This definition gives for the ideal detector. The fact that no detector performs better than an ideal one is now formulated as where should be interpreted in the matrix inequality sense, i.e. if is positive semidefinite. Proving this inequality mathematically would require analyzing how the noise cross-spectral density is related to the signal transfer function, e.g. using cascaded systems analysis,Cunningham (2000) which is beyond the scope of this work. We therefore view it as a physically motivated requirement that any realistic detector model must satisfy.

As a simple example of a nonideal detector, consider a detector that has imperfect quantum detection efficiency but infinitely high spatial resolution and perfect energy response. For this detector, and . Our definition of can thus be seen as a generalization of the energy-dependent detection efficiency .

The diagonal elements of the matrix can be interpreted as the dose efficiency per energy for detecting a spectrum change at one single energy, whereas the off-diagonal terms reflect the degree to which different detected energies interfere with each other. This interference may be constructive or destructive, depending on whether the off-diagonal terms are positive or negative. To understand why the matrix has dimensions of inverse energy unlike the conventional DQE which is dimensionless, note that a task function which equals in a small interval around a single energy gives , which can be compared to for the ideal detector. is thus proportional to if is nonsingular. This is the case for binning detectors since these have difficulty detecting a change in a narrow energy band measured in a background of noise contributions from a wide energy interval.

ii.2 NEQ and DQE in material basis

The matrix contains information about the performance of the detector for any spectral discrimination task. In practice, this matrix would be cumbersome to report in e.g. a publication, especially if the energy scale is discretized into a large number of steps. For example, if the energy variable is discretized in 1 keV-steps from 0 to 140 keV, one needs to supply a matrix as a function of spatial frequency in order to characterize the detector performance. However, the matrix can be expressed more compactly by exploiting the basis material concept i.e. the approximation that all substances likely to be present in the field of view have attenuation coefficients in a low-dimensional space, spanned by a small number of basis functions : .Alvarez and Macovski (1976) As long as one is only interested in the detector performance for differentiating materials whose attenuation coefficients lie within this space, and as long as the signal difference is small enough that a linearized description may be used, one only needs to know the restriction of the matrix to this subspace to be able to predict the performance for the relevant tasks.

For a homogeneously illuminated object and in the absence of detected scatter from the object, we can assume that the photon distribution incident on the detector is given by

(15)

where is the number of photons per area and energy incident on the object and are the basis coefficients integrated along the X-ray beam path. We will study the task of differentiating between two situations, and , both small perturbations of a homogeneous baseline . In the small-signal approximation, , and similarly for . Here, is the expected photon flux at the detector, given by (15) for , and the derivative is evaluated for . The difference in expected photon flux between the two cases is approximately . Letting denote the continuous Fourier transform of ,

(16)

and

(17)

where we have defined with . can be regarded as an element of a transformation matrix from the basis of individual energies to the basis of differential spectrum changes corresponding to differential path length changes of basis materials . We have also defined the signal difference vector in the basis , , by . Eq. 10 can now be expressed as

(18)

where the matrix elements of the matrix in basis are given by

(19)

or, in a form that is easier to compute:

(20)

Here, is the transfer function, relating the magnitude of a small relative modulation of the input photon fluence corresponding to a change in at spatial frequency to the resulting change in registered counts per area in energy bin at that frequency.

Like the for non-energy resolving detectors, has units of area. Each diagonal component of the latter specifies the number of quanta per area that an ideal detector needs to measure to achieve the same detectability as the studied detector, when the task is to detect the addition of a small amount of basis material . The matrix in basis is defined in analogy with (13):

(21)

where our notation reflects that is independent of spatial frequency. The defined in this way is a dose-normalized quantity that encodes the performance of the detector for all different tasks. For any specific detection task, the performance relative to an ideal detector is given by the task-specific DQE, which is a scalar-valued quantity:

(22)

where we have used the notation for the contribution from frequency . Eq. 22 shows that for detecting a small difference at frequency in one of the basis materials used in the decomposition, basis material , is simply given by the diagonal element of . for detecting differences due to a material with another spectral response can be obtained by expressing the DQE matrix in any basis which includes the linear attenuation of that feature as a basis function. When transforming to another basis of differential spectrum changes, transforms according to the normal rules for coordinate changes of quadratic forms, but there is no simple transformation rule for since its components are obtained as ratios of components of two matrices both of which are basis-dependent.

For a given basis , the diagonal elements of tell us how efficient the detector is for detecting a change in only one of the basis functions, while the off-diagonal terms specify the degree of constructive or destructive interference when the task involves detecting a difference in two or more basis functions. A positive off-diagonal term corresponds to a positive interference effect on detectability when both basis material path lengths are increased simultaneously.

Eq. 22 can be used to calculate the maximum and minimum for any detection task. The change of variables transforms (22) into

(23)

with maximum and minimum values given by the eigenvalues of

.

ii.3 Material decomposition

We now leave the study of detection-task performance for a moment and turn to another question: how well can the material composition of an object be quantified using measurements with the detector? One could try to estimate directly, but a more fruitful approach is to estimate its Fourier transform . A limit on estimation performance is given by the Cramér-Rao lower bound (CRLB), which gives a lower limit of the covariance matrix of the estimated parameter for any unbiased estimator and has been shown to be a useful tool for analyzing estimation accuracy for spectral CT tasksRoessl and Herrmann (2009). Although the CRLB only provides a lower bound, in practice the variance of the maximum likelihood estimator usually agrees well with the CRLB for material decomposition tasks, as has been shown in simulationsRoessl et al. (2011). The authors are not aware of a tractable form for the CRLB for correlated Poisson variables in general, but as long as there are at least a few tens of photons in every energy bin and measurement, one can approximate the joint probability distribution of the bin counts as multivariate Gaussian. This approximation is valid both with and without pileup taken into account, but in the pileup-free case both mean and variance of each individual measurement will be equal to . The components of the Fisher information matrix for a real-valued vector parameter from a measured multivariate Gaussian random variable with mean and covariance matrix is given by: Kay (1993)

(24)

The CRLB now states that for any unbiased estimator of .

In our case, and are the mean and covariance of the measured counts vector , while the parameter to be estimated contains the real and imaginary parts of the vector obtained by concatenating the vectors for all spatial frequencies in the intersection of the right half-plane and the Nyquist region. Since is the Fourier transform of a real-valued function and hence conjugate-symmetric in , we only estimate it on half the Nyquist region in order to avoid getting a singular covariance matrix. In the absence of pileup, we note that both and , and therefore also the first term in (24), are proportional to the pre-patient photon flux, whereas the second term is independent of this quantity. As long as the photon fluence is high enough, we can therefore approximate the second term with zero. We then obtain the Fisher information matrix for as:

(25)

To obtain the CRLB for the material decomposition we introduce a change of variables by defining . Here is the block matrixSchreier and Scharf (2010) , where is the identity matrix on the vector space on which are defined. As shown in appendix B, this gives the CRLB as

(26)

leading to , for any unbiased estimator of . The components of are given by

(27)

Since different spatial frequencies are independent, the Fisher information matrix is block diagonal and the block corresponding to spatial frequency is given by a rescaled version of , where is the set of basis functions in which is expressed. can thus be inverted separately for each . This shows that , which describes performance for detection tasks, and , which describes the detector performance for material quantification tasks, are very closely related. Indeed, the Cramér-Rao lower bound can be expressed as

(28)

In analogy with the definition (22) of for a detection task we can define the for quantification of basis material as

(29)

where is the CRLB of the variance of when measured with an ideal detector.

The estimated basis images are useful as they are, e.g. as a map of contrast agent distribution. However, it is also possible to make a frequency-dependent optimal weighted sum of the decomposed basis images, just as this can be done with the original energy bin images. The lower bound to the achievable variance in the basis images given by the CRLB translates to an upper bound of the achievable detectability in such a weighted basis image. This bound is given by the optimal-linear-observer detectability:Barrett and Myers (2004)

(30)

The block-diagonality of the Fisher information matrix allows us to treat spatial frequency components independently, the first term in this sum is an integral over the half-plane . Since and are both conjugate-symmetric, the allows us to extend the integral to the entire plane:

(31)

By comparing (30) with (18), we see that the same detectability is obtained by combining the different basis images as from combining the original energy bin images, as long enough statistics is available that the estimator variance will be close to the CRLB. This is a generalization of the analogous result for zero frequency derived by Alvarez.Alvarez (2010)

Another way to combine the basis images is to generate a synthetic monoenergetic imageAlvarez and Seppi (1979); Lehmann et al. (1981) as . This is useful because the image values are easy to interpret. Note, however, that the basis images in this case are combined using frequency-independent weight functions, meaning that there is no guarantee that there is an 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 two simple 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.

iii.1 Fluorescence model

As a simplified model of fluorescence K-escape, we assume that an incident photon of energy will deposit either its entire energy, with probability , or an energy with probability . is thus the energy of the fluorescence photon, and the probability of fluorescence escape is a function of with for . For simplicity, we neglect detection of secondary (fluorescence) photons in other detector elements, e.g. the detector is assumed to have coincidence logic that eliminates secondary events but is not able to reconstruct the original photon energies. We also assume that the detector has infinitely many bins, infinitely small pixels and perfect energy resolution, i.e. the fluorescence escape is the only degrading factor. The transfer function is then and an incident spectrum results in a measured spectrum . Since each photon is only registered once, the measurements at different energies are independent random variables, and the cross-spectral density is given by . Inserting these results into (9) gives, after some algebra,