Model-based pulse pileup and charge sharing compensation for photon counting detectors: A simulation study.

PURPOSE: We aim at developing a model-based algorithm that compensates for the effect of both pulse pileup (PP) and charge sharing (CS) and evaluates the performance using computer simulations.
METHODS: The proposed PCP algorithm for PP and CS compensation uses cascaded models for CS and PP we previously developed, maximizes Poisson log-likelihood, and uses an efficient three-step exhaustive search. For comparison, we also developed an LCP algorithm that combines models for a loss of counts (LCs) and CS. Two types of computer simulations, slab- and computed tomography (CT)-based, were performed to assess the performance of both PCP and LCP with 200 and 800 mA, (300 µm)2  × 1.6-mm cadmium telluride detector, and a dead-time of 23 ns. A slab-based assessment used a pair of adipose and iodine with different thicknesses, attenuated X-rays, and assessed the bias and noise of the outputs from one detector pixel; a CT-based assessment simulated a chest/cardiac scan and a head-and-neck scan using 3D phantom and noisy cone-beam projections.
RESULTS: With the slab simulation, the PCP had little or no biases when the expected counts were sufficiently large, even though a probability of count loss (PCL) due to dead-time loss or PP was as high as 0.8. In contrast, the LCP had significant biases (>±2 cm of adipose) when the PCL was higher than 0.15. Biases were present with both PCP and LCP when the expected counts were less than 10-120 per datum, which was attributed to the fact that the maximum likelihood did not approach the asymptote. The noise of PCP was within 8% from the Cramér-Rao lower bounds for most cases when no significant bias was present. The two CT studies essentially agreed with the slab simulation study. PCP had little or no biases in the estimated basis line integrals, reconstructed basis density maps, and synthesized monoenergetic CT images. But the LCP had significant biases in basis line integrals when X-ray beams passed through lungs and near the body and neck contours, where the PCLs were above 0.15. As a consequence, basis density maps and monoenergetic CT images obtained by LCP had biases throughout the imaged space.
CONCLUSION: We have developed the PCP algorithm that uses the PP-CS model. When the expected counts are more than 10-120 per datum, the PCP algorithm is statistically efficient and successfully compensates for the effect of the spectral distortion due to both PP and CS providing little or no biases in basis line integrals, basis density maps, and monoenergetic CT images regardless of count-rates. In contrast, the LCP algorithm, which models an LC due to pileup, produces severe biases when incident count-rates are high and the PCL is 0.15 or higher.

INTRODUCTION

Photon counting detector (PCD)‐based X‐ray computed tomography (CT) has great potential in many clinical applications, , , and prototype systems have shown excellent performances in phantom and clinical studies. , , , One of the challenges with PCDs, however, is spectral distortion due to pulse pileup (PP) and charge sharing (CS). , It is critical to address the effect of the spectral distortion, because many clinical applications rely on accurate spectral information PCDs are expected to output. It is impossible for the current PCDs using pulse height analysis to address both PP and CS simultaneously by tweaking their design parameters and specifications, because two major parameters, the pixel size and the pulse shaping time, have an opposite effect on PP and CS. , , For example, a smaller pixel size mitigates PP but worsens CS. A desirable strategy is to (i) design a PCD that balances the effects of PP and CS and (ii) employ an algorithm to compensate for the effect of the (remaining) spectral distortion.

Several algorithms, both model‐ , , , and data‐based, , were developed to explicitly address the CS problem with low count‐rates, and a few data‐based methods , , were developed to address both PP and CS implicitly. To our knowledge, however, there is no model‐based algorithm that can address both PP and CS altogether. We suspect that the main reason for the absence is the complexity of models. CS is better characterized than PP. A few CS models were developed and showed good agreement with physical PCDs or Monte Carlo simulators. , , , One of them, available to academic researchers, allows us to compute the expected spectrum using a simple matrix for the spectral distortion and a vector for the incident spectrum. Several PP models showed good agreement with physical PCDs or Monte Carlo simulators , , , , ; however, they are nonlinear and shift‐variant and a lot more complex than the CS models. We have developed both PP and CS models and are ready for the challenge.

Both model‐ and data‐based methods have different strengths and weaknesses. For example, model‐based methods can compute PCD data at any desirable conditions; however, they may not match the measured data completely if the model is inaccurate. Data‐based methods, on the other hand, better represent a PCD with a specific configuration and predict measured data at a specific condition; however, the accuracy of computed data at unmeasured conditions is unknown. It makes sense to combine both of the approaches eventually, and therefore, it is desirable to have a model‐based algorithm that can address both PP and CS to prepare for the integration.

The purpose of this study was to develop a model‐based algorithm that compensates for the effect of both PP and CS and evaluate the performance using computer simulations. The paper is structured as follows. In Section 2, we outline the proposed algorithm and simulation methods. We present the results in Section 3, discuss relevant issues in Section 4, and conclude the paper in Section 5. Acronyms are listed in Table 1.

TABLE 1

Acronyms

Acronyms, variables	Meaning
CS	Charge sharing
CT	Computed tomography
LC	A loss of counts
LCP	LC and CS compensation algorithm
nSD	Normalized standard deviation
PCD	Photon counting detector
PCL	Probability of count loss
PCP	PP and CS compensation algorithm
PP	Pulse pileup
XCAT	extended NURBS‐based cardiac‐torso phantom
x(r,E), x(E)	A linear attenuation coefficient of a voxel at a position r and energy E. Bold letters indicate vectors.
Φ(E)=[Φ_1(E),Φ_2(E)]	A set of basis functions used in material decomposition
w_1(r), w_2(r), w 1, w 2, w=[w_1,w_2]^T	A set of characteristic coefficients of basis functions Φ_1(E) and Φ_2(E)
v 1, v 2, v=[v_1,v_2]^T	Line integrals of characteristic coefficients (i.e., dimension‐less relative densities) of basis functions Φ_1(E) and Φ_2(E)
y_n	Noisy PCD data (counts) with multiple energy windows

, , w 1, w 2,

METHODS

We outline the proposed algorithm and an algorithm to compare within Section 2.1 and the assessment schemes in Section 2.2.

PP and CS compensation (PCP) algorithms

The proposed algorithm integrated the PP and the CS models and used an optimization framework to estimate line integrals of basis functions from measured PCD data while compensating for the effects of PP and CS. For the comparison purpose, we also combined models for a loss of counts (LC) and the CS and developed an algorithm that can address not the spectral distortion but the dead‐time loss (i.e., scaling) due to pileup (the spectral distortion due to the CS was modeled correctly). Figure 1 presents the expected spectra of the PP–CS, the LC–CS, and the CS models only, at two different count‐rates (hence, at two different probabilities of count loss [PCLs] due to dead‐time loss or PP; see Appendix for the precise definition of PCL). The algorithm that used the PP–CS model was called PCP for PP and CS compensation, and the one that used the LC–CS model was called LCP for LC and CS compensation. In the following, we outline the three key elements for both the PCP and the LCP: the system model, the cost function, and the optimization algorithm.

FIGURE 1

The 140‐kVp spectra at the probability of count loss (PCL, due to dead‐time loss, defined by Equation A11) of 0.15 (a) and 0.30 (b) after attenuation by 10 cm of adipose and charge sharing (CS). The spectra were computed by (i) the pulse pileup (PP)–CS model that takes into account the effects of PP, CS, and attenuation on both the X‐ray intensities and spectral distortion; (ii) the loss of count (LC)–CS model, which takes into account an LC due to PP and spectral changes due to CS and attenuation; and (iii) the CS model, which takes into account the attenuation and CS. The difference between (ii) and (iii) denotes the dead‐time loss (or an LCs) due to PP, whereas the difference between (i) and (ii) shows the spectral distortion due to PP. More details of the simulation settings are provided in Sections 2.2.1 and 2.2.2; both PP–CS and LC–CS models are outlined in Section 2.1.1. Results presented in Section 3 show that biases were significant with LC–CS model when PCL was higher than 0.15.

Object and system modeling

We start with modeling an object. Let x( , E) denote a linear attenuation coefficient (1/cm) of a voxel at a position and an energy E. Bold letters denote vectors in this paper. Using the material decomposition, x( , E) can be expressed as a linear combination of basis functions: where and are characteristic coefficients (i.e., dimensionless relative densities) of basis functions and (1/cm), respectively, at the position . A set of bases can be chosen from either physics phenomena, such as photoelectric effect and Compton scattering, or representative materials such as water and bone. Two basis functions suffice if the object does not contain materials with K‐edge within an X‐ray energy range and three with a K‐edge material. The line integrals of the object,  can then be computed by where are called basis line integrals (cm), and .

Now, we model PCDs. In previous studies, we developed the CS model and the PP model separately and showed that cascading multiple models for different processes could approximate the expected PCDs spectra accurately. For this study, we used the PP , and the CS models under an assumption that flat‐field X‐rays were incident onto PCDs. Let h be a function that outputs a vector of expected PCD counts for multiple energy windows, , with the effects of both PP and CS given a set of line integrals of basis functions, :  = h( ). A function h that outputs a vector of expected PCD counts with both LC and CS was also constructed by combining the LC and CS models. Both h and h will reduce to CS only at an extremely low incident count‐rate. We assumed that there was no correlation between neighboring PCD pixels nor between multiple energy windows within the same pixel and that noisy PCD data were Poisson‐distributed: where is the true basis line integrals. More details are provided in Appendix.

Cost function

The cost function for PCP and LCP, respectively, is the Poisson log‐likelihood (LL) of the corresponding PCD data with the expectation computed by the PP–CS and the LC–CS models, respectively:

Optimization algorithm

With PCP and LCP, we wish to find the maximizer that maximizes the corresponding function:

The function h (or its partial derivative) for the cascaded systems model PP–CS is computationally too expensive to evaluate on the fly and be used as a part of an iterative optimization algorithm. Thus, we precomputed the expected counts at grid points with an increment of Δv 1 and Δv 2 along v 1‐ and v 2‐axes, respectively, and employed a computationally efficient exhaustive search that takes the following three steps (see Figure 2 for pictorial description):

FIGURE 2

The proposed optimization algorithm using a multistep exhaustive search. Step 1 is to find the maximizer among s along the diagonal axis and Step 2 is to find the maximizer among s in a banana‐shaped region Ω( ) (a); and Step 3 is to find the maximizer among s in a rectangular region Σ( ) (b). Note that Δv 1 and Δv 2 are the original sampling pitch used in Steps 1 and 2, whereas Δv 1/10 and Δv 2/10 are used in Step 3.

Step 1. Candidate points along the diagonal line in the v 1–v 2 plane (Figure 2a) were used to compute LL values. The maximizer, , was then selected.

Step 2. Candidate points were those in Ω( ), a “banana‐shaped” region, which was a group of points in the v 1–v 2 plane that results in as the maximizer in Step 1. Ω( ) was determined using noisy PCD data. The maximizer of the LL, , was selected.

Step 3. Candidate points were those from Σ( ), a small rectangular region centering at with 10 times denser samples than those used in Step 2 (i.e., with an increment of Δv 1/10 and Δv 2/10, respectively, over a range of (−1.2Δv 1, 1.2Δv 1) and (−1.2Δv 2, 1.2Δv 2), respectively; see Figure 2b). The point that maximizes the LL was selected and called the global maximizer.

The banana‐shaped region Ω( (1)) used in Step 2 was preconstructed as follows. First, we generated 100 noisy data for each grid point in the v 1–v 2 plane. Second, we performed an exhaustive search along the diagonal line in the v 1–v 2 plane, which is the same as Step 1 outlined earlier. Let us call the maximizer (diag). Third, we performed an exhaustive search of the global maximum using all of the grid points in the entire v 1–v 2 plane. Let us call the maximizer (global). Fourth, we added the point (global) to the set Ω′( (diag)). By repeating the previous process for all of the noisy data, we get a set of Ω′( (diag)) for each (diag). Finally, we treated Ω′( (diag)) as a binary region in the v 1–v 2 plane and performed a morphological dilation operation to enlarge Ω′( (diag)) and obtained Ω( (diag)).

Computer simulation study

We outline the common settings and two assessment schemes in the following: slab‐ and CT‐based schemes.

Common settings

We used a cadmium telluride PCD with a pixel size of (300 μm)2, a thickness of 1.6 mm, four energy thresholds at (20, 45, 70, and 95 keV), charge cloud size of 48 μm in full‐width‐at‐half‐maximum, electronic noise of 2.0 keV, and non‐paralyzable detection with a dead‐time of 23 ns. The X‐ray spectrum was 140 kVp. Using the CS model and the other parameters presented in Ref. [19] and the PP model in Ref. [22], the expectation of the PCD outputs for both PP–CS and LC–CS models was computed for a set of basis materials attenuating the X‐rays. See Figure 1 for two spectra at two different count‐rates. A set of basis materials was adipose for representing soft tissue materials and iodine for high‐Z materials in this study; the densities were 0.92 g/cm3 for adipose and 4.94 g/cm3 for iodine; for Steps 1 and 2, the range of thicknesses was (−18 cm, 126 cm) with an increment of 0.45 cm for adipose and (−0.180 cm, 1.800 cm) with an increment of 6.0 × 10−3 cm for iodine. Negative thicknesses were included to minimize the boundary effect of the search range. When the true is on or near a boundary, estimation results will be biased because an exhaustive search, such as PCP or LCP, functions as constrained optimization, truncates the probability distribution and produces biases (with reduced noise). For Step 3, the expected counts were linearly interpolated to create 10 times denser samples (i.e., with an increment of 0.045 cm for adipose and 6.0 × 10−4 cm for iodine). It took 2.5 days in total per one tube current setting to compute h( ) for the PP–CS model for all of the sample points for Steps 1 and 2 using a 3‐GHz 6‐core Intel Core i5 2018 CPU chip with 64‐GB memory. The computation of h( ) for the LC–CS model was significantly faster. The expectation for Step 3 was obtained by linearly interpolating h( ) or h( ) generated for Steps 1 and 2.

As will be discussed later, we simulated a scenario with 4 × 4‐superpixel processing in this study, which added the outputs of 16 pixels to create one large PCD pixel. Some prototype research PCD‐CT systems (SOMATOM CounT system; Siemens Healthineers; Forchheim, Germany) output superpixel data by default, and it reduces a computational burden significantly.

Slab‐based assessment

We generated noisy data at off‐grid points not sampled in Section 2.2.1, used tube current values of 200 and 800 mA, a time duration of 400 μs per reading, and repeated the measurement 160 000 times for each . Sixteen noise realizations were then added, resulting in 10 000 noisy data for each . We performed both PCP and LCP to estimate and assessed the bias and standard deviation over multiple noise realizations. Biases larger than 2.0 cm of adipose and 2.0 × 10−2 cm of iodine were considered significant, because it was found later in the CT‐based assessment that they produced noticeable biases and artifacts. For a reference estimation noise level, the Cramér–Rao lower bound was computed for each condition using a formula for a multivariate normal distribution with off‐diagonal covariance elements being zeros. Standard deviations >10% larger than the square root of the Cramér–Rao lower bound is considered significantly large.

CT‐based assessment

We used the four‐dimensional extended NURBS‐based cardiac‐torso (XCAT) phantom version 2.0 , and generated CT images of the chest and head‐and‐neck areas, (E), at energy E = 40, 50, …, 140 keV. Material decomposition was then applied to each voxel to compute a density of adipose and iodine, , for each voxel. Cone‐beam projections (i.e., line integrals of ) were then computed by the ASTRA Toolbox , using a PCD with 4000 channels, 640 rows, and 1250 projections over one gantry rotation, to which basis line integrals of a bowtie filter, , was added to produce ( + ). The bowtie filter consists of Teflon and two bowtie shapes were used: The one used for the chest/cardiac scan was thinner and the thicknesses increased slowly toward the peripheral rays, whereas the one for the head‐and‐neck scan was thicker and the thicknesses increased quickly. The bowtie filter designs were similar to those used in clinical CT systems for body and head scans, respectively, except for the use of Teflon. Teflon was chosen for this study to make the beam hardening effect similar to soft tissue materials in contrast to metals, for example, aluminum.

The true PCD data were computed by h(  +  ) with the tube current of 800 mA and the acquisition time of 200 μs per reading (simulating 0.25 s/rot scan) for the chest/cardiac scan and 200 mA and 400 μs (simulating 0.50 s/rot) for the head‐and‐neck scan. Noisy PCD data were then calculated using the Poisson probability function with h(  +  ) as the expectation. For the reasons previously outlined, a 4 × 4‐pixel binning was employed on to create PCD data with 1000 channels, 160 rows, and 1250 projections.

We estimated (  +  ) by performing PCP with and then computed by subtracting the (known) line integrals of the bowtie filter . We reconstructed the basis function density maps by performing filtered backprojection on using the ASTRA Toolbox and synthesized CT images (E) at E = 40, 70, and 130 keV.

For comparison, we performed LCP and obtained , , and (E). In addition, we binned the true cone‐beam projection , reconstructed the true basis function density maps, , and synthesized the true CT images (E). All of the datasets included the effect of the system's sampling and resolution, which allowed for a side‐by‐side comparison.

RESULTS

Slab‐based assessment

Figure 3a,b presents biases of adipose thicknesses (v 1) estimated by PCP and LCP, respectively, with 800 mA. Very little bias was present with PCP even when the attenuation was smaller (e.g., v 1 < 20 cm) and the incident count‐rates were higher (curved arrow, Figure 3a), and PCL was as high as 0.8 (arrow, Figure 3g). In contrast, negative biases were present with LCP under the conditions (arrow, Figure 3b). This demonstrated the difference between the PP model that fully modeled the spectral distortion due to PP and the LC model that modeled an LCs only. It appeared that biases in v 1 were significant (>±2 cm) with LCP when PCL was higher than 0.13–0.15 (see Figure 3b,g), depending on the amount of iodine. We originally anticipated that the “threshold” would be ∼0.30, because the spectral distortion due to PP appeared to be minor at PCL of 0.15 (see Figure 1). It demonstrated that the impact of PP spectral distortion on spectral tasks was stronger than it appeared in spectra. Both PCP and LCP had significant biases (>±2 cm) when the expected counts per 16‐pixel binned datum were 10–120 or fewer (Figure 3h), depending on the amount of iodine (i.e., signal strengths). We do not think that these biases were related to PP and believe that they were attributed to the fact that the maximum likelihood did not approach the asymptote due to fewer counts, a finite number of energy windows, and degraded signal‐to‐noise ratios (thus, needing more photons and energy windows).

FIGURE 3

Results with 800 mA. Very little bias was present with pulse pileup (PP) and charge sharing (CS) compensation (PCP) (arrow, a) even when the attenuation was smaller (v 1 < 20 cm), the incident count‐rates were higher, and probability of count loss (PCL) was higher than 0.8 (arrow, g). In contrast, negative biases were present with loss of count (LC) and CS compensation (LCP) under the conditions (arrow, b). Biases in v 2 (c,d) had similar results, and noise (e,f) were comparable. In (h), the corresponding count‐rates were 27.8 × 106 counts/s/mm2 for 103 counts per reading and 277.8 × 106 counts/s/mm2 for 104 counts. nSD, standard deviation normalized by square root of Cramér–Rao lower bound

Figure 3c,d presents the biases in the iodine thickness estimation (v 2s). We made observations similar to the adipose estimation results previously described, except that biases at high count rates were positive for v 2, not negative. The kink observed in Figure 3d was due to clipping at the search boundary (i.e., v 1 = −18 cm).

Figure 3e,f shows the normalized standard deviation (nSD) values, that is, the standard deviation of adipose thickness estimation normalized by the Cramér–Rao lower bound at the corresponding condition. The results for the iodine estimate were very similar to the adipose estimation, thus, not presented. The nSD values of PCP were between 0.94 and 1.08 (except for two exceptions discussed later) when the expected counts were larger than 400 events, which indicated that the PCP algorithm was statistically efficient. The nSD value was erratic with (v 1 ≤ 10 cm and v 2 = 0 cm) (Figure 3e). This may be attributed to a severe spectral distortion due to PP, but other conditions with higher PCL had nSDs close to 1 (e.g., the nSD was 0.96 for a PCL of 0.56 with v 1 = 0 cm and v 2 = 0.05 cm, whereas the nSD was 1.42 for a PCL of 0.43 with v 1 = 10 cm and v 2 = 0.00 cm). We are investigating the reason for this observation. The nSD values were larger than 1 when counts were fewer; biases were present under these conditions, and PCP was not an unbiased estimator. The computed Cramér–Rao lower bounds, which denote the minimum variance of unbiased estimator, were not meaningful under these conditions when the presence of unbiased estimator was in doubt.

Figure 4 shows the results with 200 mA. The PCP had no or little biases at higher count‐rates, whereas the LCP had significant biases, albeit smaller than with 800 mA due to lower count‐rates (compare Figure 4a for PCP with Figure 4b for LCP; Figure 4b for 200 mA with Figure 3b for 800 mA). With the PCL < 0.02 at the same adipose thicknesses, the bias with both PCP and LCP was larger with 200 mA than with 800 mA due to four times fewer counts per reading.

FIGURE 4

Results with 200 mA. Very little bias was present with pulse pileup (PP) and charge sharing (CS) compensation (PCP) (arrow, a) even when the attenuation was smaller (v 1 < 20 cm), the incident count‐rates were high, and probability of count loss (PCL) was high (g). In contrast, negative biases were present with loss of count (LC) and CS compensation (LCP) under the conditions (arrow, b). Biases in v2 (c,d) had similar results, and noise (e,f) were comparable. In (h), the corresponding count‐rates were 27.8 × 106 counts/s/mm2 for 103 counts per reading and 277.8 × 106 counts/s/mm2 for 104 counts. nSD, standard deviation normalized by square root of Cramér–Rao lower bound

CT‐based assessment

Chest/cardiac scan

Figure 5 shows the true and estimated basis line integrals of a projection from the chest/cardiac scan. Both for adipose and for iodine estimated by PCP had little or no biases, whereas both and estimated by LCP had biases when rays had high count‐rates (arrows, Figure 5c,g). The corresponding PCLs were ≥0.15 (Figure 5i) for the X‐ray beams with large biases. We did not observe biases due to fewer counts in the CT scans because most data had sufficient counts with basis line integrals being v 1 ≤ 25 cm and v 2 ≤ 0.02 cm even with the bowtie filter. The noise levels of PCP and LCP were comparable to each other.

FIGURE 5

Results of the chest/cardiac scan. The true and estimated basis line integrals, v 1 (a–c) and v 2 (d–f). The profiles of the center row of v 1 (g), v 2 (h), and the probability of count loss, PCL (i). The pulse pileup (PP) and charge sharing (CS) compensation (PCP) algorithm had no visible biases (b, e, g, and h), whereas the loss of count (LC) and CS compensation (LCP) produced negative biases in v 1 in lung regions and near the body contour (c and g, arrows) and positive biases throughout the v 2 image (f and h). The PCL values for the X‐ray beams with biased LCP estimates were higher than 0.15 (i). The counts were in the range of 2.7 × 103–6.5 × 104 per datum and the maximum PCL was 0.46. The window width/center was 26/10 (cm) for v 1 and 0.04/0.00 (cm) for v 2, respectively.

Figure 6 presents basis density maps with a slice thickness of 2.8 mm reconstructed from the corresponding basis line integrals. Both the adipose ( ) and the iodine maps ( ) of PCP were very accurate for the entire imaged area, whereas those of LCP had biases almost throughout the images. The adipose map of LCP (Figure 6c,g) shows that even though large biases in projections were present for the limited areas only (arrows, Figure 5c,g), they were propagated to the entire image during the image reconstruction process.

FIGURE 6

Results of the chest/cardiac scan. The true and estimated relative density maps of basis functions, w 1 (a–c) and w 2 (d–f). Horizontal profiles of 35 mm above the center for w 1 (g) and w 2 (h). The pulse pileup (PP) and charge sharing (CS) compensation (PCP) algorithm had no visible biases, whereas the loss of count (LC) and CS compensation (LCP) images had positive biases in the w 1 map (c and g) and negative biases in the w 2 map (f and h). The window width/center was 1.0/1.0 (d.l.) for w 1 maps and 0.50 × 10−3/0.15 × 10−3 (d.l.) for w 2 maps. Biases (d.l.) with PCP were (b) 8.3 × 10−4 and (e) −1.0 × 10−5 for ROI 1 [indicated by a circle in (a)] and (b) 5.1 × 10−3 and (e) −4.0 × 10−5 for ROI 2. Biases (d.l.) with LCP were significantly larger and were (c) 1.5 × 10−1 and (f) −1.2 × 10−3 for ROI 1, and (c) 3.3 × 10−1 and (f) −3.0 × 10−3 for ROI 2. Standard deviations (d.l.) with PCP were (b) 6.7 × 10−2 and (e) 5.4 × 10−4 for ROI 1 and (b) 6.8 × 10−2 and (e) 5.2 × 10−4 for ROI 2. Those with LCP were comparable to PCP and were (c) 6.7 × 10−2 and (f) 5.5 × 10−4 for ROI 1, and (c) 7.3 × 10−2 and (f) 6.0 × 10−4 for ROI 2, respectively.

Figure 7 shows monoenergetic CT images synthesized from the basis density maps. The PCP images appeared very similar to the true images with no visible bias nor artifacts except for noise and streaks in the posterior wall at 40 keV, which appeared to be caused by fewer photons detected in the lateral views. Subtraction images did not present any unnatural patterns (not presented). In contrast, the LCP images had severe biases throughout the 40‐ and 130‐keV images (Figure 7c,i), and shading artifacts were observed near the body contours in the 70‐keV image (Figure 7f).

FIGURE 7

Results of the chest/cardiac scan. Monoenergetic chest/cardiac computed tomography (CT) images synthesized from the true density maps (a, d, g), the pulse pileup (PP) and charge sharing (CS) compensation (PCP)‐estimated maps (b, e, h), and the loss of count (LC) and CS compensation (LCP)‐estimated maps (c, f, i). The synthesized energies were 40 keV for (a–c), 70 keV for (d–f), and 130 keV for (g–i). The window width/center was 600 HU/0 HU. Biases (HU) with PCP were (b) −2.2, (e) −0.2, and (h) 0.5 for ROI 1 (indicated by a circle in (a)) and (b) −12.3, (e) −0.7, and (h) 3.4 for ROI 2. Biases (h) with LCP were significantly larger and were (c) −381.0, (f) −26.6, and (i) 97.8 for ROI 1 and (c) −929.8, (f) −82.8, and (i) 214.5 for ROI 2. Standard deviations (HU) with PCP were (b) 167.3, (e) 20.0, and (h) 45.1 for ROI 1, and (b) 158.8, (e) 21.5, and (h) 47.1 for ROI 2. Those with LCP were comparable to PCP and were (c) 170.5, (f) 20.6, and (i) 45.0 for ROI 1, and (c) 189.3, (f) 25.9, and (i) 49.1 for ROI 2, respectively.

It took 68.7 ± 4.0 min per scan for PCP to compute using a 3‐GHz 6‐core Intel Core i5 2018 CPU chip with 64‐GB memory. Biases and noise measured over multiple noise realizations and subtraction images essentially yielded no new findings, thus, not presented.

Head‐and‐neck scan

Figure 8 presents the true and estimated basis line integrals of a projection from the head‐and‐neck scan. Both and estimated by PCP had little or no biases. In contrast, those estimated by LCP had biases just outside the neck when the count‐rates were high and PCL was higher than 0.15. Note that the tube current was modest at 200 mA, and the thick bowtie filter suitable for the head scans was used for this scan. We did not observe biases due to fewer counts in the CT scans because most data had sufficient counts with basis line integrals being v 1 ≤ 20 cm and v 2 ≤ 0.01 cm even with the bowtie filter. The noise levels of PCP and LCP were comparable to each other.

FIGURE 8

Results of the head/neck scan. The true and estimated basis line integrals, v 1 (a–c) and v 2 (d–f). The profiles of the 50th row of v 1 (g), v 2 (h), and the probability of count loss, PCL (i). The pulse pileup (PP) and charge sharing (CS) compensation (PCP) algorithm had no visible biases (b, e, g, and h), whereas the loss of count (LC) and CS compensation (LCP) produced negative biases in v 1 in lung regions and near the body contour (c and g, arrows) and positive biases throughout the v 2 image (f and h). The PCL values for the X‐ray beams with biased LCP estimates were higher than 0.15 (i). The counts were in the range of 6.8 × 103–7.4 × 104 per datum and the maximum PCL was 0.27. The window width/center was 26/6 (cm) for v 1 and 0.01/0.00 for v 2, respectively.

Figure 9 shows basis density maps with a thickness of 3.4 mm reconstructed from the corresponding basis line integrals. Both the adipose and the iodine maps of PCP were very accurate throughout the imaged area, whereas those of LCP had biases almost throughout the images. The adipose map of LCP (Figure 9c,g) shows that large biases outside the neck in projections (Figure 8c,g) were propagated to the entire image during the image reconstruction process, which was consistent with the chest/cardiac scan (Figures 5 and 6).

FIGURE 9

Results of the head‐and‐neck scan. The true and estimated relative density maps of basis functions, w 1 (a–c) and w 2 (d–f). Horizontal profiles of 35 mm above the center for w 1 (g) and w 2 (h). The pulse pileup (PP) and charge sharing (CS) compensation (PCP) algorithm had no visible biases, whereas the loss of count (LC) and CS compensation (LCP) images had positive biases in the w 1 map (c and g) and negative biases in the w 2 map (f and h). The window width/center was 1.0/1.0 (d.l.) for w 1 maps and 0.50 × 10−3/0.15 × 10−3 (d.l.) for w 2 maps. Biases (d.l.) with PCP were (b) 8.9 × 10−4 and (e) −1.0 × 10−5 for ROI 1 (indicated by a circle in (a)) and (b) 8.5 × 10−4 and (e) 0.0 for ROI 2. Biases (d.l.) with LCP were significantly larger and were (c) 2.5 × 10−2 and (f) −1.9 × 10−4 for ROI 1 and (c) 1.3 × 10−1 and (f) −9.2 × 10−4 for ROI 2. Standard deviations (d.l.) with PCP were (b) 2.3 × 10−2 and (e) 1.4 × 10−4 for ROI 1, and (b) 2.2 × 10−2 and (e) 1.3 × 10−4 for ROI 2. Standard deviations (d.l.) with LCP were somewhat larger than those with PCP and were (c) 3.0 × 10−2 and (f) 1.9 × 10−4 for ROI 1, and (c) 2.8 × 10−2 and (f) 1.9 × 10−4 for ROI 2, respectively.

Figure 10 shows monoenergetic CT images obtained by PCP and LCP. Similar to the chest/cardiac scan, the PCP images had neither biases nor artifacts, whereas the LCP images at 40 and 130 keV (Figure 10g,i) had biases throughout the head and inconsistent shading/whitening artifacts near carotid arteries on the posterior side. To our surprise, the 70‐keV LCP image displayed very little biases (Figure 10f). This was a coincidence as positive biases in adipose density images and negative biases in iodine density images canceled out each other via a weighted summation when the 70‐keV CT image was synthesized.

FIGURE 10

Results of the head/neck scan. Monoenergetic chest/cardiac computed tomography (CT) images synthesized from the true density maps (a, d, g), the pulse pileup (PP) and charge sharing (CS) compensation (PCP)‐estimated maps (b, e, h), and the loss of count (LC) and CS compensation (LCP)‐estimated maps (c, f, i). The synthesized energies were 40 keV for (a–c), 70 keV for (d–f), and 130 keV for (g–i). The window width/center was 600 HU/−100 HU. Biases (HU) with PCP were (b) −0.7, (e) 0.3, and (h) 0.7 for ROI 1 (indicated by a circle in (a)), and (b) −0.5, (e) −0.4, and (h) −0.7 for ROI 2. Biases (h) with LCP were significantly larger and were (c) −541.7, (f) −1.8, and (i) 18.6 for ROI 1, and (c) −267.4, (f) −1.7, and (i) 92.5 for ROI 2. Standard deviations (HU) with PCP were (b) 38.3, (e) 6.7, and (h) 17.1 for ROI 1, and (b) 37.7, (e) 5.8, and (h) 15.9 for ROI 2. Those with LCP were somewhat larger than with PCP, and they were (c) 54.3, (f) 6.8, and (i) 24.0 for ROI 1, and (c) 54.2, (f) 6.0, and (i) 20.8 for ROI 2, respectively.

DISCUSSION

Using the PP–CS model, there was no model–data mismatch in the PCP algorithm. The PCP compensated the effect of both PP and CS successfully as long as the number of detected events was larger than 10–120 per datum (pixel), producing no or very little biases in basis line integrals even though count‐rates were high and PCL was high (e.g., >0.5). As most of the rays in the CT‐based assessment satisfied the condition, no measurable biases were present in both basis density maps and monoenergetic CT images. In contrast, the use of an LC–CS model with the LCP algorithm had model–data mismatch and, therefore, resulted in biases when incident count‐rates were higher, the PCL was ≥0.15, and the mismatch was more significant. The “threshold” for 2.0 cm of adipose bias came at the PCL of 0.15, which was lower than we originally anticipated (which was 0.30). Biases may be present with a small fraction of projection data only; however, the biases were spread over the entire image via the image reconstruction process. One could use a nonlinear image reconstruction method to suppress the spread; however, a better and more robust solution is to eliminate a model–data mismatch by using the PP–CS model (hence, the PCP algorithm).

Modeling the spectral distortion due to PP is more challenging in general than modeling for CS. Consequently, accurate PP models are more computationally expensive to evaluate than accurate CS models. It makes it extremely challenging to develop a model‐based iterative PCP method by integrating two models and using it as a part of the forward imaging process during iterations. The proposed PCP algorithm allows us to use the complex PP–CS model, efficiently performs an exhaustive search with three steps, and maximizes the Poisson likelihood of PCD data while compensating for the effect of both PP and CS. The PCP algorithm is a statistically sound method when pre‐sampling intervals are sufficiently small. Because the maximizer is one of the sampled points in Ω and Σ, the PCP is essentially the nearest neighbor operation and, therefore, adds the effect of discretization to the estimation results. It is essential to use sufficiently small sampling intervals.

One can improve the computational efficiency of the PCP. A one‐step exhaustive search would have required Poisson LL evaluation at 10 890 000 data points [h( )] for every noisy dataset, which would have required 160 days per CT scan. The three‐step PCP computes the Poisson LL at ∼3200 data points only, resulting in 68 min per CT scan. Nonetheless, 3200 data points may still be large. One could decrease the number of data points by limiting the search range and eliminating unnecessary points in both Ω and Σ. Alternatively, one could construct an iterative maximum likelihood method that computes Poisson LL at any point by interpolating precomputed data at grid points. Even though evaluating the cost function may take long time for each iteration, the number of iterations may be small and the overall computational cost for the iterative method may be less than that for the current three‐step exhaustive search.

We included the effect of a bowtie filter in the simulation using the scheme outlined in Section 2.2.3: estimating basis line integrals for both the object and the bowtie filter (  +  ), subtracting the (known) line integrals of the bowtie filter , and obtaining those for the object only, . We think it is a clever way to use single PP–CS model and include PCD pixel‐dependent incident X‐ray intensities and spectra due to the bowtie filter, instead of creating a pixel‐specific PP–CS model for thousands of pixels. When the proposed PCP is applied to a physical PCD‐CT system, pixel‐to‐pixel variations and condition‐specific deviation from the expected counts (due to, e.g., sensitivity variations) may become an issue. One may need to use pixel‐specific model parameters and develop a wrapper that absorbs such variations and converts pixel‐ and condition‐specific outputs to a standard pixel's outputs, similar to sensitivity normalization methods used for many sensors ; then single PP–CS model for the standard pixel will be applied to all of the (converted) pixels, including the effect of a bowtie filter.

The study has a few limitations. First, we did not use physical PCD data nor PCD‐CT system data. Developing a method with controlled data and applying it to an actual system poses two different challenges and one needs to accomplish each work carefully. The PP–CS model used in this study showed excellent agreement with a few physical PCDs and Monte Carlo simulation programs in the previous studies. As discussed in the previous paragraph, we will need to use pixel‐specific model parameters and develop an effective wrapper to absorb pixel‐to‐pixel variations and condition‐specific deviations. It requires a substantial amount of effort based on our previous experiences; we shall leave it for the future work. Second, tube current modulation was not employed during CT scans for simplicity. To use PCP with the tube current modulation, one will need to generate multiple PP–CS models at different tube current values, perform inter‐model interpolation to compute a PP–CS model for the tube current for each projection, and run PCP using the projection‐specific PP–CS models. Third, Poisson data with no correlation were used and the correlation between neighboring pixels and multiple energy windows of the same pixel was not simulated. As a consequence, this study could not assess the potential noise penalty of the PCP algorithm that used Poisson LL and ignored the correlation. We believe that the assessment of biases due to PP was valid, because the expectations of PCD data were accurate and recorded counts were very high because X‐rays were intense (e.g., >40 000 counts per superpixel per reading). To generate PCD data with such complex correlations, we would need to use a Monte Carlo simulator. A limited speed of Monte Carlo simulators would, in turn, limit the number of PCD pixels and projections that could be used in a study and that would have made it impossible to simulate a CT scan with a large number of PCD pixels and projections. We decided to use Poisson data with no correlation because we were interested in studying the biases and artifacts in basis density maps and monoenergetic CT images. We plan to perform the Monte Carlo simulation study as the next step.

CONCLUSIONS

We have developed the PCP algorithm that uses the PP–CS model. The PCP algorithm successfully compensates for the effect of the spectral distortion due to both PP and CS and provides little or no biases in basis line integrals, basis density maps, and monoenergetic CT images even though the PCL is higher than 0.8 with very intense X‐rays in some cases. In contrast, the LCP algorithm, which models an LC due to pileup, produces severe biases when incident count‐rates are high and the PCL is 0.15 or higher.

CONFLICT OF INTEREST

The authors CP and KS are with Siemens Healthineers. The authors have no additional relevant conflict of interest to disclose.