Linear fitting of multi-threshold counting data with a pixel-array detector for spectral X-ray imaging

Experiments and modeling are described to perform spectral fitting of multi-threshold counting measurements on a pixel-array detector. An analytical model was developed for describing the probability density function of detected voltage in X-ray photon-counting arrays, utilizing fractional photon counting to account for edge/corner effects from voltage plumes that spread across multiple pixels. Each pixel was mathematically calibrated by fitting the detected voltage distributions to the model at both 13.5 keV and 15.0 keV X-ray energies. The model and established pixel responses were then exploited to statistically recover images of X-ray intensity as a function of X-ray energy in a simulated multi-wavelength and multi-counting threshold experiment.

1. Introduction

Energy-selective X-ray imaging holds great promise at addressing major challenges in X-ray imaging and diffraction. Laue diffraction requires broad-bandwidth X-ray sources and energy assignments for each measured diffraction spot (von Laue, 1915; Schlichting et al., 1990; Clifton et al., 1997; Srajer et al., 1996). While this assignment is currently performed by analysis of the total diffraction pattern, improvements in assignments could be enabled through independent measurements of X-ray energies. In materials analysis, the transmittance of X-rays through different materials can be highly sensitive to the X-ray wavelength, particularly at wavelengths close to spectral band-edges. Consequently, spectral X-ray imaging provides contrast intimately connected to composition for materials analysis (Frojdh et al., 2013; Nik et al., 2011; Jakubek, 2007) and biomedical applications (Jakubek, 2007; Roessl & Proksa, 2007; Fredenberg et al., 2010; Boone et al., 1990; Butler et al., 2008). Multiwavelength anomalous dispersion measurements often perform poorly compared with single-wavelength measurements, due to the X-ray damage and 1/f noise attributed to serially measuring diffraction from each wavelength (Hendrickson, 2014). Many practical challenges historically associated with spectral X-ray imaging have the potential to be addressed in whole or in part through the development of photon-counting array detectors, in which a programmable counting threshold provides a means of performing energy-specific imaging (Procz et al., 2009; Llopart et al., 2002; Ballabriga et al., 2007; Broennimann et al., 2006; Henrich et al., 2009; Llopart et al., 2007).

The circuit operation of photon-counting devices compared with charge integration devices differs only in the way the amplified voltage signal is handled. Integration-type detectors integrate and report the amplified voltage over time, while photon-counting-type detectors simply count and report the number of times the voltage signal crosses a voltage threshold with a positive slope. A simplified circuit block diagram is shown in Fig. 1 for one of the pixels in a PILATUS photon-counting array [a more complete description of the per-pixel circuit diagram is provided elsewhere (Brönnimann et al., 2001)]. Since the initial number of photoelectrons produced following X-ray absorption is a random process, this spread in current/voltage per photon adds to the noise in integrating detectors that is removed upon photon counting. Even in the absence of charge sharing among pixels (Fig. 2), photon-counting-type detectors have a signal-to-noise advantage over integrating-type detectors at low photon fluxes, provided that the range of random voltages that a single photon generates is above the counting threshold.

Figure 1 After absorbing an X-ray, the resulting small voltage signal is amplified and threshold counted over a period of time before being transferred to a computer. The amount of voltage measured by the comparator per photon, as well as the amount of incident photons, depends upon a variety of parameters. Modeling these parameters allows quantitative spectral information to be extracted from multi-threshold data.

Figure 2 Model illustrating the X-ray sensing mechanism of the PILATUS detector, adapted from Kraft et al. (2009). There is no dead space in the silicon detector, causing the charge plume deposited by an X-ray photon to sometimes be fractionally detected across several charge-collecting pixels.

Multi-threshold detectors have an additional and unique advantage over integrating-type detectors in their ability to retain information about the photon energy distribution. However, accurately relating the measured counts back to X-ray photon energy remains challenging, since the amplified voltage that is measured by the threshold counter in Fig. 1 is dependent upon a large variety of factors. Specifically, the number of measured counts depends upon the voltage peak height distribution from a single photon absorbing in the active semiconducting region, the wavelength-dependent detector backscatter, the quantum efficiency for absorption, photon-counting paralysis at high count rates, the spatial distribution in the photoelectron plume over multiple pixels (Fig. 2), and the pixel-to-pixel variance in these parameters (Frojdh et al., 2013; Jakubek, 2007; Fredenberg et al., 2010; Procz et al., 2009). Consequently, the simplest approach of setting a threshold to detect one and not the other X-ray photon energy can be subject to significant errors.

In this work, each pixel in a photon-counting array was mathematically calibrated both within and beyond the normal operating parameters and fit according to an analytical model for the detector response for the purpose of performing a quantitative spectroscopic measurement. The detector parameters for the photoelectron plume spread, voltage amplitude and voltage standard deviation per incident photon were modeled on an individual pixel basis. After mathematically calibrating each of the six million pixels in the array as the model parameters, this information was in turn used to separate multi-threshold spectral images into individual contributions from different X-ray energies by fast linear fitting. In addition to successfully modeling fractional photon counts described in Fig. 2, the mathematical calibration process grants robustness to parameters such as signal amplification and offset settings, active area thickness and pixel size. The separated intensities may then be linearly rescaled for other effects, such as quantum efficiency and detector backscatter.

2. Theoretical framework

An analytical model for the peak-height distribution as a function of threshold setting was developed based on the assumption of a two-dimensional Gaussian spatial distribution in charge following X-ray absorption, which is consistent with previous models and simulations (Julien & Kadda, 2012; Trueb et al., 2012). The distribution has two contributions: one from X-rays in which the large majority of the plume lies within the area of a single detector pixel and another from plumes spreading over two or more pixels. In previous treatments, efforts to address the issue of charge distribution over multiple pixels have included setting the counting threshold equal to roughly half the mean voltage obtained from the plume centered on a pixel, such that only the pixel with the majority of charge would register a count (Broennimann et al., 2006; Mathieson et al., 2002). While highly successful in applications focused on detection of intensity for monochromatic radiation, this strategy becomes untenable in applications targeting energy discrimination, in which the mean voltage of a centered plume is also a variable in the analysis. Even in monochromatic X-ray detection, the half-maximum threshold approach can result in bias from corner effects, in which the plume is distributed over more than two pixels. As the dimensions of the array elements continue to be reduced in size, such effects are likely to become increasingly important.

The fraction of charge expected to be observed by a single pixel from any photon absorbed by the sensor due to this plume effect can be described by the two-dimensional surface integral across the x and y spatial coordinates of the silicon area above the pixel surface. For a pixel of width w, height h and a boundary distance of consideration b for photons hitting near the pixel active area, this surface integral and its solution is

The solution of (1) gives the fraction of charge F expected to be observed from a photon landing at point (μ_x, μ_y) with a Gaussian standard deviation plume of uncertainty of σ_psf. The spatial extent of the charge plume is affected by the detector thickness and reverse bias voltage, both of which vary with manufacturing tolerances. In addition, the per-photon plume width itself can vary depending on the depth within the sensor in which the photon is absorbed, which is only accounted for phenomenologically in this model through the effective measured value of σ_psf. Assuming a uniform X-ray intensity field over the pixel surface from a monochromatic light source, a probability density function (PDF) of fractional contributions was obtained by Monte Carlo simulation.

The peak height of the voltage transient observed by the threshold-counting electronics from a 100% contribution of a single photon's charge is proportional to the energy of the incident X-ray photon, and has a normally distributed peak current/voltage. For fractional contributions of photon energy, the fractional contribution multiplicatively weights the photon's energy contribution. The resulting PDF of voltage peak heights is then described by a multiplication of random variables, where a random weight is applied to a voltage of normal uncertainty. The resulting overall PDF for the peak voltage observed by the counting electronics from any direct or proximal X-ray photon strike can be obtained by Monte Carlo simulation or derived numerically using the product distribution integral. Generally stated, for Z = XY where X and Y are two independent random variables with PDFs f_x and f_y, the PDF of the product f_z is given by equation (2) (Springer & Thompson, 1966),

A representative peak height distribution is shown in Fig. 3, along with the complimentary cumulative density function (cCDF) corresponding to the anticipated measured probability of observing a count for a photon absorption event.

Figure 3 (Left) The PDF and complimentary CDF (cCDF) of the amount of equivalent charge observed by a single pixel from a single 13.5 keV X-ray photon. The Gaussian hump on the right of the PDF represents collection of the full deposited charge following X-ray absorption, and all probabilities to the left of this hump represent fractional observations of the X-ray deposited charge plume. This distribution has a total of three parameters, with the photon energy standard deviation and photon energy mean describing the position and width of the Gaussian hump, and the size of the point spread affecting the depth of the trough to the left of the Gaussian hump. (Right) The cCDF was directly measured by the PILATUS detector from a 13.5 keV monochromatic source at many comparator threshold voltage levels (Vth in Fig. 1) converted to equivalent units of incident X-ray energy, ranging from 7.5 keV to 21 keV. After averaging all pixel responses at every threshold, the data were fit to the model cCDF with excellent agreement.

3. Methods

All data were acquired at beamline 17-ID, IMCA-CAT, at Argonne National Laboratory. Diffuse scattering of vitreous ice was measured with a Dectris Pilatus-6M single-threshold detector at several detector comparator threshold levels (V_th in Fig. 1). The detector had a pixel size of 172 µm × 172 µm, a silicon active area thickness of 320 µm, and a bias voltage of 150 V. Five-second exposure times were taken at each detector threshold for both 13.5 keV and 15 keV incident X-ray energies in a standard lattice (153 ns between X-ray pulses). Absolute detector voltage thresholds at each pixel were automatically calibrated through Pilatus' internal voltage trim system to maintain threshold accuracy. A low-gain input amplifier setting (τ = 125 ns) was used for all measurements. The resulting internal voltage threshold levels are denoted here as equivalent thresholds in units of keV, which describes the equivalent X-ray energy that would deposit this mean level of voltage. The 13.5 keV incident energy measurements were serially taken with equivalent threshold energies from 7.5 keV to 21.0 keV in steps of 0.5 keV with a detector distance of 0.700 m. The 15 keV incident energy measurements were serially taken on a later day with a new ice sample with equivalent threshold energies from 7.5 keV to 20.9 keV in steps of 0.2 keV with a detector distance of 1.000 m. In all cases, the incident photon flux was kept low enough to ensure a low probability of pulse pile-up affecting counting results.

All data analysis was performed in MATLAB with custom software. Data files were read using the MATLAB macros package for cSAXS (Paul Scherrer Institute). ImageJ was also used to view data files using a plugin (CBF reader plugin, written by JLM).

4. Results/discussion

The validity of the mathematical model was assessed by fitting it to the detector response averaged across all pixels for weak scattering of ice from 13.5 keV photons, shown in Fig. 3. The novel utility of this mathematical approach is that, rather than satisfying the default intuition that each photon should be counted once as a full contribution, the unbiased intensity can be recovered by allowing for fractional photon counting (Coldwell et al., 2001; Burrell & Rousseau, 1995; Leydesdorff et al., 2013), in which the fraction of the photon at that pixel is included in the measured PDF and cCDF. Further, signal contributions from photon strikes near pixel corners that were previously lost by setting the threshold to the equivalent of 50% of the incident photon energy can be now be properly accounted for by knowledge of the cCDF. The corresponding signal-to-noise ratio increase will be particularly advantageous for counting arrays with small pixel sizes, which exhibit corresponding increases in the probability of the electron plume spreading over multiple pixels.

The counts were detected as a function of threshold for every pixel in a 6 Mpixel array detector, and nonlinear fits were performed for each pixel to determine best-fit values of the one-photon mean, one-photon standard deviation and point spread size, representative results of which are provided in Fig. 4. Raw data consisted of a stack of images of the counts acquired under different threshold settings from a diffusely scattering sample of vitreous ice. After performing the fitting, the best-fit values from each pixel were representable by a set of three images of the fitting parameters, also shown in Fig. 4. The nonlinear fits recovered the experimental cCDFs remarkably well, although with a significant diversity in the fitting parameters across the array. This process was repeated for both 13.5 keV and 15.0 keV X-ray energies.

Figure 4 Visual flowchart of the mathematical calibration process that was performed at each X-ray energy. (Top) Images of diffuse X-ray scattering of monochromic 13.5 keV photons from ice at several comparator threshold voltage levels (Vth in Fig. 1), converted to units of electron volts, ranging from 7.5 keV to 21 keV. One representative pixel was chosen from the image stack for displaying in a plot along with their fit to the cCDF model. Noise in the data trace is dominated by 1/f noise from the ∼2 min required to change detector voltage thresholds at all six million pixels. (Bottom) Each of the six million pixels were independently fitted to the cCDF model, with the three parameters of the fit displayed as images. These three pixel-specific calibration parameters define the instrument response of each pixel to X-ray color.

The noise in the data trace in the single-pixel plots of Fig. 4 is dominated by 1/f noise from the ∼2 min required to change detector thresholds. Inspection of the cCDF measurements acquired at adjacent pixels produced covariant noise features, suggesting slow drift in the overall instrument, including the position of the X-ray source in the sample, the X-ray beam intensity or the X-ray optics rather than dominance from white Gaussian noise in the detector. This is actually a favorable condition for two reasons. First, the 1/f drift was sufficiently low to still allow complete analysis of the array. Second, it suggests that significant noise reduction would be expected by assigning four independent thresholds for a 2 × 2 section, with corresponding improvements in the confidence in energy assignment.

Inspection of the three parameter images in the bottom of Fig. 4 reveals several interesting trends. First, the variation in response characteristics for each of the 12 rows and 5 columns of Pilatus 100k modules that make up the Pilatus 6M are clearly apparent from the recovered parameters. The subtle differences in the construction of each sensor component results in considerable diversity in the value of the one-photon mean, one-photon standard deviation and point spread size as a function of location across the detector array. In general, there was a strong correlation between the one-photon mean and one-photon standard deviation, and a strong anti-correlation for the point spread size. For example, the Pilatus 100k modules shown as the first column of expanded pictures in Fig. 5 from rows 3 and 4 of column 1 are shown as the whitest among the modules for both the one-photon mean and one-standard deviation, and shown as the darkest among the modules for the point spread. Some small regions in the detector slightly deviated from these trends as is shown in the second column of expanded images in Fig. 5, where pixels with a relatively low voltage per photon have relatively large standard deviation and point spread. The high sensitivity of the mathematical calibration was able to detect this relatively subtle deviance, given that all pixels remain within 2% of the average detector parameters.

Figure 5 (Left column) Blow-up of the 100k modules in column 1, rows 2 through 5, featuring the typical correlation between one-photon mean and one-photon standard deviation, and the anti-correlation with point spread. (Right column) Blow-up of a small section of pixels which slightly deviates from the trends depicted in the first column where a smaller voltage per photon resulted in a larger variance and point spread, from the 100k module in row 4 and column 3. The high sensitivity of the mathematical calibration was able to detect this deviance, despite all pixels remaining within 2% of the average detector parameters.

Overall, the values obtained from the fitting at 15 keV agreed well with those produced independently from the 13.5 keV data, indicating good reproducibility in the fitting approach. This agreement is particularly noteworthy given the large differences in the X-ray intensity measured across the array, indicating that the detector analysis can be reliably performed independently of detected intensity, provided a sufficient number of counts are recorded at each location to perform the nonlinear fitting.

Following calibration, the measured counts as a function of threshold can serve as a `hyperspectral' signature to separate out the multi-threshold image stack into independent images at each X-ray energy through simple linear fitting. In brief, the measured cCDF `spectrum' can be written as a linear combination of the cCDFs for the 13.5 keV and the 15.0 keV sources,

In equation (2), the measured cCDF is expressed as a linear combination of two other cCDFs, and their amplitude weights a₁ and a₂ are the parameters of the linear fit. This expression can be rewritten in matrix form by defining C to be a matrix of cCDF `spectra', then inverted to solve for the best-fit values of the parameters a. Images of each coefficient correspond to the best-fit estimates of each separated X-ray energy image,

Because the matrix C is known in advance following the initial nonlinear fitting calibration step, the combined inverse and transpose matrix can also be determined a priori. In the present case of discrimination between two known energies, [(C^TC)^-1C^T] is a 2 × n matrix, where n is the number of thresholds used to generate c_meas.

Experimental validation of the algorithm was performed by separating a composite image generated by summation of the counts acquired at 13.5 keV and 15.0 keV, shown in Fig. 6. The monochromatic images were recovered using counting data from just eight of the acquired threshold data. The target of eight thresholds was selected based on the potential ease of acquisition. Assuming two independently adjustable thresholds per pixel as with the Medipix3 technology (Pennicard et al., 2012), eight-threshold detection could be accomplished by pooling 2 × 2 pixel spots. Using equation (4) the best-fit component images were recovered. In general, the number of unique parameters related to the photon energy must be equal or less than the number of thresholds used to extract them. Further, the confidence in the values obtained from the fits will increase as the number of information-carrying threshold measurements increases. Although equation (5) is cast with respect to detection of two monochromatic responses for simplicity, the approach would work comparably well for linear functions of those responses (e.g. derived from principal component analysis or linear discriminant analysis).

Figure 6 (Top image stack) Multi-threshold measurement of two-color diffraction was simulated by summing the 13.5 keV and 15 keV diffraction image stacks at eight different pairs of matching detector comparator voltage thresholds. (Middle plots) The three parameters previously recovered for each color (Fig. 3) were used to re-separate the number of counts contributed by each of the two colors, demonstrated here for a single pixel. Comparison against the original monochromatic dataset shows good agreement. (Bottom images) The fit was repeated for all pixels and shown in the top row of images. Comparison against the original monochromatic dataset qualitatively shows good agreement.

Overall, the recovered images were in good agreement with the initial monochromatic images, as indicated in the error analysis shown in Fig. 7. All but a small fraction of the pixels recovered the correct amplitudes of the original monochromatic images to within 20% accuracy for a photon energy separation of just 11%. Pixels out of 15% tolerance in absolute accuracy are indicated in images to the right in Fig. 7. Comparing the error images with the monochromatic images in Fig. 6 suggests that the multiplex disadvantage is the major source of error in the decomposed images. In brief, pixels with high brightness in one monochromatic image will introduce Poisson counting uncertainty proportional to the square root of the number of photons counted. This noise will be distributed over both extracted images, potentially increasing the noise in the separated channel to a value significantly greater than the square root of the number of counts. From inspection of Fig. 6, the bright diffraction spot from the 13.5 keV contribution in the bottom right generated substantial noise in the 15 keV image, and vice versa for the bright halo in the center of the 15 keV contribution.

Figure 7 Quantitative assessment of the error of the fitted values compared with the original monochromatic dataset. 1/f noise associated with ∼2 min between threshold measurements contributes the majority of the error in this fitting, which would be substantially reduced by simultaneous acquisition of all measurement thresholds. Even with this noise present, the vast majority of all pixels are within 20% of their original monochromatic dataset values. Pixels that contain more than 15% absolute error were those subject to multiplex disadvantage (the color that contributes the most signal also contributes the most noise).

Although the focus of this work is on discrimination between two known energies, the consistency in the extracted array parameters from the nonlinear fits suggests that reasonable estimates for the anticipated parameters at intermediate photon energies could be obtained by interpolation. In this manner, calibration of the array using a relatively small number of judiciously selected photon energies could serve to enable analysis at arbitrary energies within the calibration window. With the exception of sharp features associated with absorption edges (e.g. the K-edge in Si around 2 keV), the X-ray responsivity of silicon photodiodes varies smoothly with photon energy (Fraser et al., 1994). The mathematical calibrations for many colors could then be interpolated and simultaneously fit.

From the information recovered in Fig. 3 for the one-photon distribution, the results of the fitting analysis can be used to assess the strengths and limitations of energy-dispersive spectral X-ray imaging for a broader class of applications. From inspection of Fig. 3, the one-photon spread in the voltage distribution corresponds to a relative standard deviation of ∼0.5–0.6 keV. Using a standard definition of resolution, R = in which μ_i and σ_i are the mean and standard deviation, respectively, of the ith peak, the two peaks are considered resolvable when R ≥ 1. This definition indicates the potential ability to resolve photon energies separated by ∼0.8 keV at the energies considered. This spectral resolution is lower but still reasonably close to that achievable by commercial single-channel energy-dispersive spectrometers used in X-ray fluorescence capable of resolving differences of several hundred eV. For comparison, wavelength-dispersive benchtop X-ray spectrometers routinely exhibit considerably better energy resolution on the order of a few eV. While the imaging approach described herein may be directly applicable to Laue diffraction, diffraction using multiple simultaneous monochromatic sources, and some X-ray fluorescence spectral imaging applications, the spectral resolution achievable in energy-dispersive X-ray spectrometers will complicate some of the more challenging applications only possible with high spectral resolution. Consequently, the trade-off between high-spectral/low-spatial-temporal resolution of a wavelength-dispersive system versus high-spatial-temporal/low-spectral of an energy-dispersive imaging system should be assessed on an application-specific basis. The primary advantage of the present approach is the potential ability to perform dual energy-dispersive imaging at high frame rates (current Pilatus detectors operate at kHz frame rates) with a Poisson-limited signal-to-noise ratio (SNR) exhibiting no read noise, on six million parallel channels.

5. Conclusions

Analytical models for mathematically calibrating the peak height distribution combined with fractional photon counting allowed accurate recovery of the measured counts as a function of comparator threshold, from which monochromatic X-ray images were independently recovered from a single combined dual-wavelength image. The model successfully included fractional photon counts as a means to compensate for photoelectron point spread, as shown by the fitting agreement between the model and averaged detector response in Fig. 3. The advantage of mathematically calibrating each pixel individually was observed to account for the wide variety of individual pixel characteristics as shown in Fig. 4. In general, the chromatic intensities of the pixels recovered from the fitting were in excellent quantitative agreement with the intensities within the initial monochromatic images shown in Fig. 6. Residual error in the fitting shown in Fig. 7 comes from 1/f noise during the raw data acquisition, and multiplex measurement disadvantage. Noise from 1/f will be substantially reduced for parallelized measurements.

Charge sharing was treated stochastically in this present analysis, meaning that the measured result is not completely deterministic. The statistical analysis described allows one to relate the measured counts as a function of threshold value back to the most probable value for the mean intensity at that pixel based on the observables. The information lost by charge-sharing in conventional photon counting is recovered through the statistics of the multi-threshold response.

The presented method of multi-threshold chromatic imaging can be further improved by compensating for X-ray wavelength-dependent effects, such as quantum efficiency and detector backscatter. Fortuitously, these effects can be accurately compensated for after the application of the presented method by linearly scaling the captured chromatic intensities with the relative magnitude of these effects at each wavelength. Because the detector response at each of the two wavelengths was directly used as the mathematical calibration points, no assumptions or interpolations were necessary regarding the responsivity of the detector as a function of wavelength. Mathematically, the number of X-ray energies that can be independently discriminated cannot exceed the number of thresholds used in the analysis, although in practice reliable discrimination is preferably overdetermined (more thresholds than values to be fit). However, additional independent information regarding the X-ray energy (e.g. via Laue diffraction) may reduce the total number of unknowns to allow analysis where the X-ray energies may include several values or even be continuous functions.

The present analysis was performed on a pixel-by-pixel basis, with no incorporation of the overall pattern of intensities or image analysis. Significant further improvements in the energy-discrimination capabilities can emerge from additional consideration of spatial arrangements. In diffraction analysis, Laue diffraction can provide independent assessment of the X-ray photon energy through the pattern of observed diffraction peaks. In X-ray imaging, the object resulting in image contrast typically contains rich additional information accessible through a host of different image analysis algorithms. These options are all beyond the scope of the present study, focused exclusively on mapping the per-pixel detector response. However, they could be integrated into application-specific measurements incorporating multi-threshold imaging detection with pixel-array detectors.