Framework for Photon Counting Quantitative Material Decomposition

In this paper, the accuracy of material decomposition (MD) using an energy discriminating photon counting detector was studied. An MD framework was established and validated using calcium hydroxyapatite (CaHA) inserts of known densities (50 mg/cm3, 100 mg/cm3, 250 mg/cm3, 400 mg/cm3), and diameters (1.2, 3.0, and 5.0 mm). These inserts were placed in a cardiac rod phantom that mimics a tissue equivalent heart and measured using an experimental photon counting detector cone beam computed tomography (PCD-CBCT) setup. The quantitative coronary calcium scores (density, mass, and volume) obtained from the MD framework were compared with the nominal values. In addition, three different calibration techniques, signal-to-equivalent thickness calibration (STC), polynomial correction (PC), and projected equivalent thickness calibration (PETC) were compared to investigate the effect of the calibration method on the quantitative values. The obtained MD estimates agreed well with the nominal values for density (mass) with mean absolute percent errors (MAPEs) 8 ± 11% (9 ± 15%) and 4 ± 6% (9 ± 14%) for STC and PETC calibration methods, respectively. PC displayed large MAPEs for density (27 ± 9%), and mass (25 ± 12%). Volume estimation resulted in large deviations between true and measured values with notable MAPEs for STC (40 ± 90%), PC (40 ± 80%), and PETC (40 ± 90%). The framework demonstrated the feasibility of quantitative CaHA mass and density scoring using PCD-CBCT.


I. INTRODUCTION
S PECTRAL computed tomography (CT) enables discrimination of different materials and improved tissue contrast by utilizing polychromatic X-ray attenuation information [1]- [3]. Consequently, it has been applied for several diagnostic tasks such as bone and calcium removal from CT angiography [2], characterization of gout [4], and assessment of myocardial blood flow [5], [6]. In health care, spectral CT is currently performed using dual-energy CT (DECT) scanners that employ spectrally indiscriminate energy-integrating detectors (EIDs). The maximum number of separable tissues in DECT material decomposition (MD) with EIDs is three [7]. Furthermore, DECT suffers from spectral overlap which reduces the accuracy of MD. To increase the MD accuracy in DECT, either the peak kilovoltage (kVp) difference between low and high energy spectra could be increased, filtration could be changed between acquisitions to improve energy separation or the patient dose could be increased. As photon counting detectors (PCDs) can discriminate photons of different energies through pulse-height analysis, they enable multi-energy spectral CT [8], quantitative MD [9], and K-edge imaging [10], [11]. Furthermore, PCDs with more than three bin counters allow decomposition of more than three materials [12].
In MD, multi-energy data are decomposed into basis materials with distinct attenuation properties. MD in image space, i.e. reconstruction space, utilizes the effective energies of a multi-energy measurement to separate the CT reconstruction into basis materials [9], [13]. In contrast, projection based MD models the X-ray attenuation process, allowing the multi-energy projection data to be decomposed into material-specific projections [11], [14]. Performing MD in projection space is preferable as the modeling of X-ray attenuation theoretically eliminates beam hardening [15]. MD based on multi-energy data has been successfully used for quantification of iodine [8], [9], gadolinium [9], and lipid [16] content, as well as the amount of Calcium Hydroxyapatite (CaHA) [7], [17], which is required for the assessment of bone mineral density. These studies have shown strong correlations with measured and true concentrations, confirming that spectral CT can be used for quantitative MD and tissue discrimination.
With PCDs, several phenomena, such as radiation scattering, beam hardening effects, limited detective quantum efficiency, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ detector dead time, inaccuracies with spectral models, and pulse pile-up, deteriorate the energy-resolving capabilities and the accuracy of tissue quantification [18]. It is therefore essential to develop both accurate models for detector characteristics [19] (e.g. charge transport properties [20]), and suitable correction techniques, such as calibration material measurements dedicated for PCDs [21], [22], to counter the factors degrading spectral discrimination. On the other hand, calibration techniques may influence MD accuracy [23], and care should be taken when choosing the appropriate methods.
For energy-independent EIDs, conventional flat-field correction using an air scan and dark field image is often sufficient, since it corrects fixed-pattern noise, such as gain variations and scintillator screen dust, commonly encountered with EIDs. With clinical CT systems, however, additional complicated correction methods beyond simple offset and air correction are needed to achieve high-fidelity imaging performance over the range of tube voltage settings necessary for clinical imaging. Due to the nonlinearity and variation of inter-pixel detection efficiencies with respect to energy, flat-field correction is not adequate for PCDs either [24]. Commonly PCD calibration techniques rely on X-ray transmission measurements of calibration materials that are used to correct the measurement data. For example, the Signal-to-equivalent Thickness Calibration (STC) [24], [25] method has been widely used for PCDs, while the Polynomial Correction (PC) has been used for decomposing calibration materials from dual-energy measurements [26].
On the image reconstruction side, the commonly used analytic CT reconstruction technique, filtered backprojection (FBP), may not provide sufficient image quality in projection based spectral CT due to increased noise in the material-decomposed projection data [10]. Therefore, the choice of reconstruction technique has a significant impact on the image quality of tissue-specific reconstructions. Instead of the FBP method, penalized weighted least squares (PWLS) approaches are often used for CT reconstruction [27] and tissue-specific reconstruction [10], [28] as they accommodate noise statistics and regularization into the reconstruction scheme.
Currently, the severity of coronary artery calcification is evaluated using CT calcium scoring method known as Agatston score [29]. The score is given by a slice-specific product of calcification area and a weighting factor that depends on the maximum Hounsfield Unit (HU) value of the calcification determined from a non-invasive CT scan [29]. The total Agatston score of a coronary artery is the sum of the slice-specific scores. Agatston score is associated with the risk of coronary events and its negative predictive value for coronary artery disease is high [30], [31]. Additionally, calcium hydroxyapatite (CaHA) mass [32] and calcification volume [33] are commonly evaluated metrics in cardiac CT studies, with the CaHA mass being more reproducible than Agatston score [34]. Agatston score does not quantify the concentration of coronary calcification exactly but rather estimates the concentration through HUs. However, beam hardening artifact and the selected peak kilovoltage affect the HU estimate and hence the Agatston score [35], generating a need for a more accurate quantitative methodology to estimate calcium concentration and mass. MD in projection space models the beam hardening process and theoretically, it could provide a robust alternative for calcium scoring.
This study aims to establish a framework for MD using dual-energy PCD measurements and apply it for quantifying coronary artery calcium (CaHA) in an experimental setting. Furthermore, the effect of the calibration method on the framework is investigated by comparing three different calibration techniques; the STC, PC, and a new Projected Equivalent Thickness Calibration (PETC). Finally, the quantitative accuracy of the developed framework is validated by characterizing CaHA insert densities, masses, and volumes in a cardiac rod phantom.

A. Experimental Setup and Imaging Phantom
Cylindrical CaHA inserts (n = 12) with varying densities (50 mg/cm 3 , 100 mg/cm 3 , 250 mg/cm 3 , and 400 mg/cm 3 ), diameters (1.2 mm, 3.0 mm, and 5.0 mm) and a fixed height of 7 mm were placed inside a cylindrical (9 cm diameter) cardiac rod phantom (008C, CIRS, Inc., Norfolk, VA) (Fig. 1). The 5.0 mm insert was positioned in the right coronary artery, the 3 mm insert in the apex, and the 1.2 mm insert in the left anterior descending coronary artery of the phantom (Fig. 1).
The tungsten X-ray source of a C-arm (Philips BV29, Philips Healthcare, Netherlands), with inherent filtration equivalent to 3.0 mm of Aluminum at 75 kVp, was operated at 100 kVp and 3.0 mA. The focal spot size of the source was 0.6 mm. A PCD (Flite FX15, XCounter AB, Danderyd, Sweden) was used as the radiation detector, and low tube current was used to mitigate the effects of dead time and pulse pile-up in the PCD. The PCD and the X-ray source were mounted on an optical board and combined with a motorized rotation stage (NR360S/M, Thorlabs, Inc., Newton, New Jersey) on which the cardiac rod was placed, this setup enabled PCD cone beam CT (PCD-CBCT) (Fig. 1). Sourceto-object and source-to-detector distances were set to 54.0 cm and 85.48 cm, respectively. The angular velocity of the rotation stage was set to 3 • /s, resulting in 1 mAs exposure per 1 • of rotation with a total scan time of 120 s. The frame rate of the detector was set to 24 frames/s to avoid count saturation. Eight consecutive frames were averaged, resulting in 360 projections over a 360 • angle of rotation and a total tube exposure of 360 mAs. No anti-scatter grid was used. For each CaHA density, measurements were repeated three times, and inserts were repositioned into the rod between measurements.
The used PCD was a flat-panel direct conversion detector utilizing Cadmium Telluride (CdTe) of 750-μm thickness as a conversion medium. The detector consisted of 24 separate tiles, each containing a 256 × 128-pixel matrix with 100-μm pixel pitch resulting in a tilewise active area of 2.56 cm × 1.28 cm. These tiles formed a matrix with a 100-μm gap between tiles. The total active area of the PCD was 5.13 cm × 15.47 cm. Prior to MD, the detector matrix was binned by a factor of two, resulting in an isotropic 200-μm pixel size. Taking geometric magnification into account, the isotropic voxel size of the reconstruction was 126 μm × 126 μm × 126 μm.
The PCD had two 12-bit counters per pixel with adjustable energy thresholds that differentiated the incident photons into two energy bins as follows; total energy (TE) bin containing photons with higher energy than the lower threshold, and high energy (HE) bin including photons whose energy exceeded the higher threshold. Low energy (LE) image was obtained through subtraction LE = TE -HE. The energy thresholds of the PCD were set to 10 keV and 50 keV. The purpose of the low energy threshold was mainly to remove electronic noise, and the high threshold was chosen so that the counts between the energy bins were approximately equal. Furthermore, a charge sharing correction implemented in the PCD software was also used to improve spectral separation [36].

B. Beam Hardening and Scatter Correction
Three different calibration methods (STC, PC, and PETC) were implemented for beam hardening and scatter correction in this study. Each of these methods transforms the measured counts of the raw data into calibration material thicknesses (Fig. 2 a) -b)). After this transformation, the thicknesses are projected into corrected counts (Fig. 2 c)).

1) Signal-to-Equivalent Thickness Calibration (STC):
A common beam hardening and scatter correction approach is the STC method [24], [25], where a series of calibration measurements, usually performed with polymethyl methacrylate (PMMA) blocks or Aluminum (Al) sheets, are used to find the calibration parameters (o, b, and a) for each detector pixel by fitting an exponential function to the calibration data. We used PMMA-based STC for calibrating the measurements because PMMA has more similar scattering properties to soft tissues compared to Al. This method translates the known calibration material thicknesses (t i c ) into measured counts where c denotes calibration, i denotes the energy bin, and the coefficient a is considered as an effective attenuation coefficient incorporating an experimental correction for scatter and beam hardening. The attenuating polychromatic photons were approximated to follow the monochromatic exponential law be at [24]. The calibration parameter o relates the counts measured due to scattering within a specific calibration thickness interval. With this approximation, we obtain the attenuation function (1). Once the calibration parameters have been solved for every detector pixel, the phantom measurement (y i m ) can be corrected into PMMA equivalent thickness (t i m ) by computing This approach is repeated for each energy bin i separately.

2) Polynomial Correction (PC):
Another calibration method is the polynomial correction (PC) technique [26]. With this method, combinations of two calibration materials (Al, PMMA) are measured in the calibration process. The second order polynomial approximation for two energy bins i is , and y i c,0 is the flat-field image. The polynomial constants (c and d) are determined from a least-squares fit of the calibration measurements. Phantom measurement is calibrated into thicknesses of PMMA and Al by directly applying (3) and (4).

3) Projected Equivalent Thickness Calibration (PETC):
In this study, we suggest a modification of the STC method, by accounting for two calibration materials (PMMA and Al) in an attempt to incorporate the differences in beam-hardening and scattering properties of calcifications and soft tissues more accurately. To do so, we solved the corrected count estimate y i, in energy bin i through a process where the STC estimate at first calibration basis (Al) Al t Al (5) was projected to PMMA equivalent thickness t i The superscript i in PMMA equivalent thickness highlights the fact that although the same basis material thickness (t) was used for each energy bin, the projected equivalent thickness (t i Al→P M M A ) deviated between energy bins due to the different attenuation properties of the calibration materials. The corrected count-estimate (y i, ) was obtained by combining the current thickness-estimate of PMMA and the projected equivalent thickness (projection from Al to PMMA) In order to decompose the measurement data into a combination of calibration basis materials (t = (t Al , t P M M A )), the obtained calibration maps for PMMA and Al (obtained through the STC technique) were used in the minimization problemt where y is the measured counts and y denotes the current calibration count estimate, yielding the decomposed combination of basis materials for a specific pixel. y and y are vectors containing the counts of different bins. The Gauss-Newton (GN) method [37, Ch. 9] with 1000 iterations was used to minimize (8). y was updated during each iteration through the subsequent utilization of (5) - (7). This method is visually illustrated in Fig. 3 and referred in the rest of the study as Projected Equivalent Thickness Calibration (PETC) method. For each calibration method, the calibration material thicknesses were selected to cover the thickness variations of the rod. Therefore, PMMA slabs with thicknesses of 0 cm, 5.3 cm, and 10.7 cm were measured and used in the STC and PETC methods, and Al sheets with thicknesses of 0 cm, 2.1 cm, and 3.1 cm were measured for the PETC method. For PC, combinations of PMMA thicknesses (0 cm, 5.2 cm, and 10.4 cm) and Al thicknesses (0 cm, 0.05 cm, 0.1 cm, . . ., 0.6 cm) were used. As PC utilizes combinations of Al and PMMA thicknesses in calibration, thinner Al plates were used compared to PETC. PC also fits more calibration parameters than PETC, and thus Al thicknesses were measured more showing the STC map of PMMA. Using the calibration material estimate t = (t Al , t PMMA ), a) translates the Al thickness (t Al ) to estimated counts using (5), b) illustrates translation to PMMA basis, in c) the count estimate from Al basis is projected to PMMA equivalent thickness (6), d) adds the projected equivalent thickness and thickness estimate of PMMA (t PMMA ) and e)-f) gives the final count estimate y i, (t ) using (7). densely in PC compared to PETC. Calibration measurements were performed using the same imaging parameters as in the phantom scans.
To assess the scatter-correction mechanisms of the calibration techniques used in this study (STC, PC, and PETC), the post-correction scatter-to-primary ratio (SPR) with the cardiac rod phantom was evaluated for each method. The term "postcorrection SPR" is used to distinguish between the actual physical SPR of the measurement and the post-calibration SPR. For reference, we also evaluated the physical SPR from a flat-field corrected measurement. The measurements for evaluating the SPR are described in the Appendix A. 1 In summary, the thickest part of the cardiac rod was imaged with and without a lead (Pb) collimator with a rectangular 1.4 mm × 2.8 mm × 115 mm opening. The uncollimated rod measurement was corrected with the calibration techniques (STC, PC, and PETC). The collimated measurement (m 1 = P) contained only primary photons (P), whereas the corrected uncollimated counts (m 2 = S + P) were a composition of primary and scattered (S) photons. The SPR was calculated as

C. X-Ray Attenuation and Detector Model
The X-ray attenuation through M basis materials and the resulting detector counts in detector pixel (k) and energy bin i (λ i k ) was modeled using the attenuation law where R i (E) is the detector energy response of energy bin i at energy E. (E) and D(E) are the incident source spectrum and detection efficiency for energy E. μ m is the mass attenuation coefficient of material m. A k,m is the projected mass density (PMD) along path l k where ρ m is the density of m th basis material in basis reconstruction voxels x m . The energy-dependent mass attenuation coefficients (μ m (E)) were obtained from the NIST database [38]. Source spectrum ((E)) was simulated using the spektr toolbox [39] for MATLAB (v. 9.4, The Math-Works Inc., Natick, MA, 2018) and the detection efficiency of CdTe was computed from the attenuation law. The energy response of the PCD was modeled using a particle transport simulation package (vGATE 7.2) [40] and a charge transport model [41] that propagates the induced charge in the detector crystal (CdTe). A particle transport simulation package utilizing the Geant4 toolkit [42], [43] with PENetration and Energy LOss of Positrons and Electrons (PENELOPE) code system [44] was used to model the particle interactions within the PCD. A pixelated (100 × 100 pixels) CdTe detector with 100 μm × 100 μm pixel and 750-μm CdTe thickness was created. Compton and Rayleigh scattering, photoelectric effect, bremsstrahlung, ionization, and atomic de-excitation (i.e. fluorescence and Auger effect) were taken into account in the model. The monochromatic energy responses of the PCD were modeled using γ -source simulations within the energy range of 10 keV -100 keV with 1 keV increments. The types of the interactions, interaction times, energies, and locations were tracked for each individual primary photon and daughter particle generated during the interactions of the primary photons within the CdTe.
The propagation of charge cloud generated by photon hits was modeled with MATLAB using a charge transport model similar to [41]. The bias voltage was set to −400 V [41] and the charge transport properties of CdTe, e.g. electron and hole mobility, resistivity and pair creation energy, were obtained from [45]. An eight-neighbor charge sharing correction with 500 ns signal integration time was incorporated in the model. The signal integration time was chosen based on prior research for a CdTe detector [45, p. 57].
The accuracy of the resulting energy responses obtained from the particle and charge transport models were experimentally validated by sweeping the detector high energy threshold from 11 keV to 80 keV with 1 keV increments and measuring the response of the PCD for an Americium-241 ( 241 Am) (E γ = 59.5 keV, activity 3.7 GBq) radioisotope [36], [46]. Full width at half-maximum (FWHM) and main peak location were evaluated from the 241 Am response to study the energy resolution of the PCD. Peak-to-valley ratio and the ratio between peak area and valley area were also computed to assess the robustness of charge sharing correction and the effect of the K-escape photons of Cd and Te on the measured and simulated energy responses. These parameters were separately calculated for each detector tile and the simulated response model. The simulated and measured 241 Am responses were filtered using a moving average filter with a denominator coefficient of 1 and window length of 5. For the evaluation of FWHM and the ratio between peak area and valley area, a Gaussian function was fitted to the filtered main response peak. The peak area was defined as the area of the Gaussian fit. Valley area was determined as the area outside the Gaussian fit, i.e. the area of the low-energy tail of the response. Similar to [47], the X-ray source spectrum (100 kVp, 3.0 mA) was swept with 2 keV increments using the C-arm to validate the overall accuracy of the energy response model over the entire source spectrum. The lower energy threshold was fixed at 10 keV in both 241 Am and C-arm sweep measurements.

D. Material Decomposition
After correcting the phantom measurements to combinations of calibration material thicknesses (for STC: PMMA, for PETC and PC: PMMA and Al) for each pixel k, they were propagated to corrected counts λ corr (10), where λ corr,1 and λ corr,2 are the corrected counts of LE and HE bins, respectively (Fig. 2 c)). Basis materials (PMMA and CaHA) were then decomposed (Fig 2 d)) by obtaining the PMD ( A k = (A k,P M M A , A k,Ca H A )) that solves the problem A GN method was used to find the MD as reported previously [48]. The GN approach relies on a two-step minimization process with the update rule where the current decomposition estimate at pixel k is updated with a GN step ( p where the Jacobian (J k ) is updated during each GN iteration ( j ) with T denoting the transpose. Gradient descent algorithm [49] was used to solve the least-squares problem (12). 100 GN iterations were required for convergence in decomposition. The Jacobian in (14) was block-diagonalized as in [48] to decompose each pixel simultaneously.
As the dual-material calibration methods (PC and PETC) decompose the measurement data into two basis materials (PMMA and Al), we could alternatively map the (PMMA, Al) content directly to (PMMA, CaHA) content with a linear transformation [50]. With this method, the mass attenuation coefficient of CaHA is estimated as a linear combination of the mass attenuation coefficients of PMMA and Al over the used X-ray energy range. This approach would remove the need for detector model, count propagation, and material decomposition steps (Fig. 2. c) -d)). With STC, however, this approach would not be possible as it decomposes the counts into thicknesses of one basis material (PMMA) in different energy bins. The count propagation and material decomposition steps are, therefore, required with STC to obtain the CaHA content. As we wanted to have a consistent material decomposition workflow between the calibration methods for a reasonable comparison between their calcium quantification accuracy, we did not use the linear mapping from (PMMA, Al) basis to (PMMA, CaHA) content in this study. However, we included a comparison between the density quantification accuracies of our framework and the linear translation method for PC and PETC in Appendix B. 2

E. Image Reconstruction and Quantitative Calcium Scoring
A PWLS algorithm with the Huber penalty [51], [52] and variance based weighting was applied to reconstruct CT data from the decomposed CaHA basis projection images (Fig. 2 e)). As the material decomposition algorithm is computationally time-consuming, the projection data was binned by a factor of two and eight subsequent frames were averaged to speed up the material decomposition process to a reasonable level. The variance was evaluated from these averaged pixels, and the weighting matrix ( −1 b m ) was constructed by diagonalizing the reciprocals of these variance-estimates for each binned pixel.
To estimate the variance of the decomposed projections, CaHA basis projections at full 513 × 1547 resolution and for each measured frame were required. As iteratively decomposing the projection data at full resolution was considered too time-consuming, the non-iterative areal bone mineral density (aBMD) [53] was used to approximate the material decomposition. aBMD can be calculated from the dual-energy measurement using , (15) where A Ca H A denotes the aBMD and μ is the mass attenuation coefficient. The energies for LE and HE bins (E 1 and E 2 ) were calculated as weighted average energies of the LE and HE spectra. Variance for each iteratively material decomposed pixel was estimated from the aBMD decomposed frames within a 2 × 2 binning area and from eight subsequent frames. These pixels were included in the variance-estimation as they are averaged prior to MD.
The PWLS minimization problem for the m th basis material in matrix form was the same as in [54]: where the reconstructionx m was obtained through our own implementation of the gradient descent method with Barzilai-Borwein step size update with 100 iterations. The projection matrix W was constructed using the ASTRA (v. of reconstructed voxels, and R is the Huber penalty where penalty parameter β and the threshold δ were set to 400 and 50 mg/cm 3 , respectively. As the Huber penalty adapts the regularization based on pixel differences, it can simultaneously provide robust edge preservation and efficient noise-removal [52]. Reconstructions were computed using MATLAB. The calcifications were segmented using the seeded region growing algorithm [57]. A detailed description of the method can be found in the Appendix C. 3 Briefly, the algorithm contains three major steps: 1) Seed point selection using a fixed threshold. The threshold was determined by multiplying the mean background density with an experimentally selected fixed multiplication factor (1.25).
2) The choice of tolerance for the region growing algorithm using a histogram-analysis of the local neighborhood of the seed pixel. 3) Region growing algorithm initiating from the seed pixels and using the tolerance obtained from 2).
Calcification volumes were evaluated from the segmentations and the mass of the calcification was determined as a product of mean density and volume. In the evaluation of mean calcification density, a morphological erosion was performed on the three-dimensional segmentation mask. A sphere of five pixels radius was used as the structuring element. Morphological erosion removed partial volume effect and blurring near calcification boundaries arising from the limited spatial resolution of the system. Finally, the mean density within the eroded mask was determined. Morphological erosion was not used in the evaluation of mass or volume.
The quantitative accuracy for each calibration method was evaluated. Mean value, coefficient of variation (CV), and standard deviation (SD) were determined from three repeated measurements for mass, density, and volume. Regression slope, Pearson correlation (r 2 and p-values), boxplots with percent error, Bland-Altman (BA) plots, and mean absolute percent error (MAPE) were used to illustrate correlation and deviation between nominal and measured values. MAPE was calculated as where X t is the measured value and Y t is the nominal value.
The analyses were carried out using MATLAB.

A. Validation of the Detector Model
When the sweep measurement using the PCD was compared with the simulated gamma photon response of 241 Am (E γ = 59.5 keV), a clear drift to higher energy was observed in the location of the main response peak (Fig. 4 and Table I).  The drifts were between 7 keV -14 keV and depended on the detector tile. No large deviation was observed in the pixelwise responses within specific tiles. Similar observations were also made from the simulated and measured tilewise 100 kVp C-arm spectra (Fig. 4). The simulated and measured low energy spectra (10 keV -60 keV) agreed well, whereas substantial differences were evident at high photon energies (60 keV -100 keV). In addition to the observed high energy differences between simulated and measured energy responses, large tile-to-tile variations were observed in the energies of the 241 Am γ -peaks (Table I and Fig. 5). Apart from the varying locations of the main response peaks, the measured low-energy tails (10 keV -43 keV) of different tiles agreed well (Fig. 5) and the low SD (0.07) for the ratio between peak areas and valley areas indicated similar characteristics between responses of each tile. These variations, FWHMs, peak locations, peakto-valley ratios and peak-area-to-valley-areas for measurement and simulation are summarized in Table I.

B. Evaluation of Scatter-to-Primary Ratio
The calculated post-correction SPRs for different calibration methods are tabulated in Table II. The physical SPRs for the flat-field corrected measurements were higher when compared to the post-correction SPRs of the calibration methods used in our study. Notably, the STC and PC techniques resulted in negative post-correction SPR values for each energy bin. Both negative and positive values were observed with our proposed PETC method.
To illustrate the insufficiency of the flat-field correction with PCDs in these scattering conditions, we visually compared it with STC, PC, and PETC, and observed that the detector tiles were clearly visible in the flat-field corrected images (Fig. 6).

C. Quantitative Material Decomposition
Overall, the decomposition framework provided robust separation between the basis materials (PMMA and CaHA) for each calibration method used in this study (Fig. 7). Due to the limited spatial resolution of the imaging setup, the reconstructions of 1.2 mm inserts appeared blurred (Fig. 8), and the MAPEs of the density estimates for 1.2 mm inserts were 37 ± 3%, 52 ± 2%, and 38 ± 4% for STC, PC, and PETC, respectively. Thus, 1.2 mm inserts were excluded from density analysis. However, they were included for volume and mass analysis. Additionally, all 50 mg/cm 3 inserts and the 100 mg/cm 3 insert with the diameter of 1.2 mm were not separable from the cardiac rod (Fig. 8), and therefore could not be analyzed.
The measured densities for the calibration methods (STC, PC, and PETC) are illustrated with Table III and BA plots (Fig. 9 a) -c)). Strong linear correlations were observed between measured and true CaHA densities for both STC (regression slope = 0.98, r 2 = 0.97, p < 0.0001) and PETC (regression slope = 0.87, r 2 = 0.99, p < 0.0001). PC displayed moderate correlation with the density estimates (regression slope = 0.66, r 2 = 0.99, p < 0.0001). Although the STC method showed a strong correlation, the decreased coefficient of determination and higher CV (7.9%) indicated increased variability with the density estimates when compared to PETC (CV = 6.2%). However, STC outperformed the PC method, which portrayed an even higher CV (19.9%). These interpretations were confirmed by the BA plots (Fig. 9 a) -c)). Additionally, the obtained densities with the STC method delineated clearly between different insert diameters (3 mm and 5 mm) (Fig. 9 a)), suggesting a relationship with insert diameter and STC calibrated density. For PETC and PC,  no similar relationship was observed, and the measured densities were similar between different insert diameters ( Fig. 9 b) -c)). However, calcium density could not be accurately resolved with the PC method ( Fig. 9 b)). PETC resulted in the lowest MAPE (4 ± 6%) in the quantification of calcium density, when compared to 8 ± 11%, and 27 ± 9%, for STC and PC, respectively.
Distinct variations in quantitative accuracy were observed in the evaluation of CaHA mass between the methods.  Strong linear correlations were identified for both STC (regression slope = 1.10, r 2 = 0.99, p < 0.0001) and PETC (regression slope = 0.96, r 2 = 0.99, p < 0.0001), but a slight increase in CV was observed with STC (CV = 10.7%) when compared with PETC (CV = 7.1%). Similarly to density estimation, PC showed the weakest correlation with the nominal masses (regression slope = 0.68, r 2 = 0.99, p < 0.0001). Furthermore, for a specific insert diameter (5 mm and 3 mm), the boxplots for STC illustrated increasing percent mass error estimation with increasing density (Fig. 10). This trend was distinguishable for each density with diameters 5 mm and 3 mm, whereas the inverse effect was observed for PETC. No similar phenomenon was seen with the PC method. The overall mass quantification errors between PETC (MAPE = 9 ± 14%) and STC (MAPE = 9 ± 15%) were highly similar. Visual assessment of the BA plots ( Fig. 9 g) -i)) and the CVs (10.7%, 36.4%, and 7.1% for STC, PC, and PETC), however, confirmed PETC to be the most robust in mass evaluation. PC exhibited high quantification error in the evaluation of mass (MAPE = 25 ± 12%).

IV. DISCUSSION
In this study, we established and validated a framework for MD with our experimental PCD-CBCT system and STC, PC, and PETC calibration techniques using quantification of coronary calcium as a case study. Strong correlations were observed in the evaluation of CaHA density and mass for STC and PETC methods, with PETC showing the best accuracy in the quantification of CaHA inserts. PC did not provide accurate quantification of CaHA density, volume or mass.
The STC, PC, and PETC methods reduced scatter compared to flat-field correction as observed in the reduced post-correction SPR values. The SPRs with the flat-field correction were comparable with breast CBCT and a 10 cm 50% glandular -50% adipose composition phantom using an air gap of 41.5 cm at 80 kVp (SPR = 21%) [58]. In these scattering conditions, the STC, PC, and PETC methods will reduce scattering effectively. For CBCT of the heart, however, the scattering from the torso further increases the SPR, thus complicating the imaging process, but this assessment was left for subsequent studies.
In the presence of extensive scatter in CBCT, using only one calibration material in STC led to lower quantitative accuracy when compared to PETC. This was most likely due to the assumption in STC that all tissues have the same beam hardening and scattering properties as the calibration material. As CaHA has a higher effective atomic number compared to PMMA, it has an increased probability to attenuate through photoelectric effect compared to PMMA. Consequently, for the same amount of measured counts PMMA images contain more scatter than the CaHA image. As the STC method corrects scattering using an estimate obtained from PMMA measurements, the CaHA images of the cardiac rod containing less scattering were overcorrected. This is supported by the negative post-correction SPR values with STC indicating an overcorrection of scatter with the cardiac rod. The overcorrection manifested as thickness dependence in the estimation of calcification density, and thickness and density dependence in the quantification of mass as observed in Fig. 9 and Fig. 10.
Although PMMA is a widely used calibration material for soft tissues due to both having similar attenuation and scattering properties, calcified tissues behave differently from PMMA. Al differs significantly from PMMA as a calibration material and affects the measured spectrum more through photoelectric effect and beam-hardening. As hypothesized in the introduction of the PETC method, PMMA accurately characterized the scattering and beam hardening from soft tissue equivalent rod and a combination of PMMA and Al effectively modeled the increased fraction of photoelectric effect with CaHA. The simultaneous utilization of two calibration materials in the PC and PETC methods, therefore, reduced the diameter dependence of density estimates (Fig. 9 b) -c)).
Of the two dual-material calibration techniques, PC could not accurately quantify the calcifications but underestimated the CaHA densities (Table III). We hypothesize that this is related to the differences in calibration measurements of PETC and PC; in PETC we measured the calibration material thicknesses separately without combination measurements, whereas combinations of calibration materials were measured in PC. The combination measurements in PC warrants further investigation since the order of the calibration plates might also affect the results. For instance, Al sheet will act as a filter for the scattered photons originating from the PMMA plate if the Al sheet is positioned after the PMMA plate. However, future research focusing specifically on the calibration process is required to reliably explain these observations.
The importance of the choice of calibration technique can also be observed through a visual comparison between the calibration methods (Fig. 6). The detector tiles were clearly visible in the flat-field corrected images, whereas STC, PETC, and PC achieved comparably good tile uniformity. These results are consistent with a more comprehensive comparison study on tile-variation between flat-field correction and STC correction, where STC outperformed the flat-field correction in tile uniformity [59].
Overall, a more accurate thickness estimation might have been obtained by measuring more calibration thicknesses for the STC and PETC methods. The choice of three calibration measurements for STC was based on prior research with an XCounter detector [59], and the instructions given by the PCD manufacturer. Three calibration thicknesses were also used for PETC to make the results comparable with STC. Therefore, three calibration points per calibration material in PETC and STC were chosen. As the PC method contains more fit parameters than STC and PETC, we measured more calibration points (36 thickness combinations) for PC. In the energy response evaluation of the PCD, we observed excellent tilewise agreement between low energy tails (energies 10 keV -43 keV) of the 241 Am response. These tails originate from the interactions between detector crystal and the L α , L β and L γ -escape X-rays of 241 Am daughter particle 237 Np (13.9 keV, 17.8 keV, and 20.8 keV, respectively) and K-escape photons of Cd (26.7 keV) and Te (31.8 keV) (Fig. 4, Table I). Due to these photon escape effects, the low-energy tails did not solely arise from charge sharing. Therefore, the charge sharing correction of the PCD software could not remove the low-energy tails entirely.
Clear inter-tile variation and high-energy drifts were observed in the tilewise energy responses of the PCD (Fig. 5). High-energy drifts might have been caused by a pulse-height discrimination calibration that was not optimized for higher energies. Moreover, the inter-tile variation was likely caused by the lack of tile-specific transfer functions, since they have a critical role in the reliable conversion of pulse-height in the PCD. We did not account for these drifts as a comprehensive drift-model characterization would have required several monochromatic radiation sources (e.g. from synchrotron) within the used energy range, and this was considered to be out of the scope of this study. Disregarding these drifts, the response model and measurements had similar peak areato-valley areas (simulated model = 2.19, measured = 2.04 ± 0.07) ( Table I) and exhibited the same FWHM (14.2 keV). Nevertheless, the substantial deviation between the simulated and true spectra at high energies generated a major source of error for quantitative MD.
As a case study for our quantitative MD framework, we applied it for the quantification of coronary calcium using a simple cardiac rod phantom. Currently, calcium scoring is performed with single-energy CT. However, the utilization of dual-energy information in CT angiography has improved the calcium quantification accuracy compared to single-energy CT [60]. Although our results indicated that MD using PCDs enables accurate quantification of coronary calcifications, these results are not yet applicable to diagnostic cardiac CT. Current PCDs are severely limited in their ability to process high pho-ton fluxes encountered with clinical CT protocols complicating their use in clinical cardiac imaging [18]. Furthermore, cardiac motion needs to be addressed for long scan times with low flux. This aspect needs to be investigated with motion correction schemes e.g. using ECG gating [61] in future studies.
Several limitations should be addressed when interpreting the results. Since the cardiac rod fitted into the small field of view of our PCD, it was imaged instead of a torso phantom. This was done to avoid the interior tomography geometry [62], which can affect the quantitative values. Also, the lack of suitable anti-scatter grid or post-collimation affected the measured counts through increased scattering. Addition of beam-collimation would most likely have improved the accuracy of the framework for each calibration technique.
Furthermore, the non-optimal positioning and the large focal spot of the X-ray tube introduced geometric inaccuracy to the imaging system. This inaccuracy manifested as blurring of the calcification boundary in the reconstructions (Fig. 8). Although the tube could not be positioned exactly perpendicular to the detector surface and through the rotation axis, the resulting error was mitigated by mounting the C-arm to the optical table using a lift jack and mounting poles. The blurring of the calcifications due to the large focal spot could also have been mitigated by reducing the magnification by decreasing the object-to-detector distance. However, our main concern was to reduce the effects of X-ray scattering. Therefore, we decided to utilize a large air gap between the object and the detector [63] to reduce the effects of scatter at the detector.
Finally, the image reconstruction method may have introduced bias to the obtained quantitative results. Particular care was therefore taken during image reconstruction and analysis; the same regularization parameter, chosen after a visual comparison between reconstructions obtained with different regularization parameters, was used for each measurement. We also observed ringing artifacts in tissue-specific reconstructions (Fig. 7). These artifacts were caused by dead pixels located in the tile borders.
Regarding noise estimation in the PWLS reconstruction, more accurate noise-estimate could have been obtained by estimating the noise directly from the material decomposed frames. However, this procedure would have resulted in high computational demands due to the requirement of iterative decomposition of each measured frame at full resolution (513 × 1547). Therefore, variance-estimation was obtained from aBMD decomposed projections to achieve reasonable computation times. A noteworthy observation is that (10) reduces to the aBMD estimate (15) when assuming a monochromatic dual-energy measurement with the same energies as the weighted average energies of the polychromatic LE and HE measurements. Therefore, aBMD was a justifiable approximation to use in the estimation of variance in the material decomposed projections.
The choice of a particular segmentation technique can introduce inaccuracies to the quantitative results. Commonly in Agatston scoring, a fixed segmentation threshold is used. Due to the large variation between the used insert densities (100 mg/cm 3 -400 mg/cm 3 ), using one fixed threshold did not provide visually sufficient segmentation accuracy, and would have required manual intervention to correct the segmentations. Thus, we decided to use a locally adaptive region growing method to properly account for the varying insert densities. The segmentations were visually assessed to confirm the accuracy of the segmentations. We believe that the main reason for the observed high errors in volume quantification was likely related to the limited spatial resolution of the imaging system. Future work will address a validation of the framework using a torso phantom in interior tomography geometry. This validation is required for proving a clinically relevant demonstration of scatter-resilience with our framework. Future research should also investigate the accuracy of the framework in the quantification of contrast agents such as gadolinium and iodine which has been successful in prior research [9]. Finally, the minimum radiation dose required for accurate MD combined with the optimal low-dose reconstruction technique warrants further investigation.

V. CONCLUSION
We decomposed PCD measurements successfully to CaHA and PMMA bases and reconstructed them into material-specific images. Moreover, we were able to obtain accurate quantification of CaHA density and mass in a cardiac rod phantom with the introduced multi-energy framework. In light of our results, we propose that MD using PCDs can be used for accurate coronary artery calcium mass and density quantification. More studies are, however, needed to verify these results in a clinical setting for cardiac CT.