Photon-counting detector step-wedge calibration enabling water and iodine material decomposition

Photon-counting detector computed tomography (PCD-CT) has demonstrated improvements in conventional image quality compared to energy integrating detector (EID) CT. PCD-CT has the advantage of being able to operate in conventional and spectral mode simultaneously by sorting photons according to selected energy thresholds. However, to reconstruct spectral images PCD-CT requires extensive calibration and specifically fine-tuning a spectral response. This response is then used to perform material decomposition (MD). We propose a step-wedge phantom made of water and iodine to calibrate a prototype PCD-CT system. Four methods were tested and compared based on calibration accuracy and CT image quality. The exhaustive PCD response was not well calibrated, but a reduced model was defined that was able to perform accurate water/iodine MD and to reduce the ring artifact intensity. The impact of the number of energy bins (from two to seven) was also studied. The number of bins did not affect the spectral accuracy. However, compared to the two energy bin configuration, the seven bin configuration decreased the noise by 10% and 15% in the water and iodine maps, respectively. The model was tested on ex-vivo tissue samples injected with iodine to demonstrate the results of the water/iodine MD on biological materials.


Introduction
X-ray computed tomography (CT) is a powerful diagnostic imaging tool widely used in clinical practice [1].Conventional CT provides images of gray-scale contrast which is able to reveal human anatomy but does not enable accurate quantitative analysis.Spectral CT is an evolution of conventional CT that uses spectral information to provide accurate quantitative measurements [2].The spectral information comes from the energetic dependence of the linear attenuation coefficient (LAC) of the object of interest.Spectral CT is based on the measurement of at least two energies in order to study this dependence.This can be performed with either a dual-energy CT (DE-CT) system equipped with energy integrating detector (EID) or a photon-counting detector CT (PCD-CT) system.DE-CTs have now been commercialized over the last twenty years [2] and the first PCD-CT system has been recently commercialized [3] based on a CdTe semi-conductor material.Several other PCD-CT prototype systems are under development using different semi-conductor materials, such as CdZnTe (CZT) [4][5][6], or Si [7].Scintillation detection technology using silicon photomultipliers (SiPM) -1 -which provide similar energy resolution and count rate capability have also been proven effective for photon-counting CT in medical imaging applications [8].
PCD-CT has been shown to provide better image quality in conventional imaging compared to EID-CT [1,9,10].PCD-CT has the ability to reduce noise compared to EID-CT at the same spatial resolution level.Moreover, as an EID can only measure one energy, DE-CT requires a modification of the conventional CT system to measure two energies.This can be done by either generating incident spectra at two different energies or through the use of a dual-layer EID [2].However, PCDs can directly sort the incoming photons according to their energy to provide several measurements in one acquisition.The number of measurements (or energy bins) available is defined by the number of thresholds provided by the PCD.The energy bin number is commonly between two and eight for PCDs designed for medical applications [11].These thresholds can be exploited for several spectral applications such as metal artifact reduction [12], K-edge imaging [13] or material quantification [2].
One of the main clinical applications of quantitative CT is CT angiography which involves the injection of iodine-based contrast agent.In this case, spectral CT enables the separation of the materials in the image into three principal categories [2]: soft-tissue, bone and iodine signal.After the materials have been classified into these three categories, additional spectral images can be generated such as iodine concentration maps, virtual non-contrast (VNCo) images, or virtual non-calcium (VNCa) images [14].There are several ways to perform this material decomposition (MD).The first is the reconstruction of a low-and high-energy CT image.Due to the attenuation properties of iodine and bone, these materials can be separated from soft-tissue and then identified.The second method is the computation of a MD relying on the Alvarez-Macowski technique [15].The result of the MD is two basis material density maps (in absence of a K-edge material).Commonly a low-and high-attenuating basis materials are selected to ensure the stability of the MD.From the MD, each material can either be directly quantified if it is one of the basis materials or represented as a linear mixture of the basis materials.
Different strategies have been developed to perform MD and are categorized based on whether they are computed before, after, or during CT image reconstruction.Pre-reconstruction (or projectionbased) methods [16] decompose the detected sinograms into basis material sinograms.In the second step, these basis material sinograms are reconstructed.This group of methods requires perfect co-registration between energy bins.Co-registration means that all energy bins match and are acquired with the same geometry.Thus pre-reconstruction methods are not adapted to all DE-CT configurations (e.g rapid switching kilovoltage peak (kVp) DE-CT systems do not have a perfect co-registration between energy bins).On the contrary, post-reconstruction (or image-based) methods [17] first reconstruct each of the energy bins.A MD is then computed from the reconstructed images.The issue with this approach is that the artifacts present in the reconstructed images, such as ring or beam hardening artifacts, will also be present and enhanced in final basis material density maps.The last approach performs the MD and CT reconstruction in the same step, so-called one-step [18] (or joint-based) algorithms that have a long computation time.
Pre-reconstruction methods are a natural choice for PCD-CT, as the PCD provides perfect co-registration of the energy bins.Decomposing the energy bin sinograms into basis materials requires knowledge of the incident spectrum and detector response.The detector response model is the probability of an incoming photon to be detected in an energy bin according to its energy.This model depends on the detector properties such as the semi-conductor material thickness, the pixel size, -2 -and the electric field intensity [11].Because of crystal impurities, this model varies from one pixel to another.Sophisticated simulation tools can provide a first order estimation of the detector model, but MD performed with such a model will not generate accurate quantitative results.The model needs to be refined from experimental results, such as based on synchroton source acquisitions [16] or the use of a step-wedge phantom [19].
A step-wedge phantom is a calibration tool composed of known lengths and materials [20].Such a phantom provides the ability to acquire a series of transmission measurements that can be exploited to adjust the detector model.Step-wedge phantoms have been studied to determine x-ray source spectra [20] or detector response models [19,21].In both cases it has been demonstrated that the calibration process requires a good initial estimate in order to converge to the true result.
Step-wedge phantoms are commonly made of two materials: a highly attenuating material such as aluminum or polyvinyl chloride (PVC) and a low attenuating material such as polymethyl methacrylate (PMMA) or polyethylene (PE).
In this work we propose an innovative step wedge phantom made of different water lengths and iodine concentrations.All of the measurements were performed with a PCD with up to seven energy bins.Four different projection-based MD methods were investigated, all calibrated with the step wedge phantom acquisitions.The first method is a polynomial based method that directly estimates the material sinograms [28].The other three methods calibrate a forward spectral model that is then inverted to perform MD.The first forward model is an exhaustive physical model which estimates the system's spectral response.The second model is a model based on a reduced number of parameters [27].The last method is a specific case of the reduced model where the calibrated model is linear [22].After the step-wedge calibration, the methods were applied to CT acquisitions and resulted in a water/iodine MD.The different models were compared by measuring various iodine concentrations as well as the noise in CT images in the seven energy bin configuration.The best performing model was then studied further to investigate the impact of a decreased number of energy bins on CT image quality.

Experimental system
The system used in this work [12,13,29], is prototype system composed of an x-ray tube (MXR 160/22, Comet Technologies, San Jose, CA, U.S.A.), a rotation stage, and a PCD panel (figure 1).The source to iso-center distance (SID) and source to detector distance (SDD) are 322 mm and 578 mm, respectively.The PCD is manufactured by Redlen Technologies (Saanichton, BC, Canada) and the semi-conductor crystal is made of CdZnTe (CZT).The semi-conductor material thickness is 2 mm and the pixel pitch is 330 μm.The detector is composed of 16 sub-modules, each containing 24 × 36 pixels, arranged to create an array of 24 × 576 pixels.The lowest threshold of the detector is set to 24 keV and six additional thresholds can be set to any value > 24 keV, leading to seven energy bins.The x-ray source was operated at 120 kV and filtered with 3 mm aluminum and 2 mm titanium.These values were selected to reproduce medical imaging conditions [30,31].Indeed, the beam spectrum was simulated with the Spekpy module [32] to compute its mean incident energy of 68.2 keV and its equivalent aluminum first half-value layer (HVL) of 9.5 mm.The source includes a tungsten based anode and was operated at 1 mA.In addition to the source filtration, a narrow 0.2 mm -collimation was -3 - set on the source to irradiate an active area corresponding to 4 × 576 pixels, leading to a field-of-view of 100 × 100 × 0.6 mm.Given the low current, the filtration, and the narrow beam collimation, non-linearities such as the pile-up effect or scattering could be neglected in this setup.The PCD threshold calibration [33] was performed shortly before all the acquisitions measured in this work.

Pixel heterogeneity response
CZT-based detectors suffer from spectral distortion which is caused by physical effects such as charge-sharing or x-ray fluorescence losses.The effect of this spectral distortion is displayed in figure 2. It illustrates the spectral distortion between the incident source spectrum and the measured spectrum, as well as the heterogeneities between pixels.This non-uniformity between pixels is even more preponderant at lower energies (< 35 keV).The detector response was simulated with a Monte Carlo code [34] based on the detector characteristics.The simulated model has already been demonstrated to enable accurate energy threshold calibration [33].It provides a first order approximation of the detector response.However, it suffers from two flaws.First, it needs to be refined to better match experimental measurements (as depicted by the difference between the red and black continuous curves), and second, it needs to be adjusted pixel per pixel, each one having a slightly different response.An accurate detector response enables the spectral distortion to be inverted to provide accurate MD.

Thresholds calculation
In a first section of the results, four methods are compared using data exploiting all the thresholds available, i.e. seven energy bins.The method presenting the best results will then be tested with a reduced number of bins.Aiming to perform a fair comparison, the optimal threshold values are computed for each energy bins number.To reduce the number of possibilities, only configuration with a regular spacing between thresholds are investigated.For each threshold number and sampling is computed the Cramer-Rao lower bound (CRLB) [22,35,36].The CRLB depends on the selected statistical estimator and physical model.Here the CRLB was computed for a Poisson likelihood approach and the exhaustive model (eq.(2.3)), with the simulated detector response displayed in figure 2. The CRLB was evaluated for the MD basis material lengths of 10 cm of water and 2.1 cm of iodine at 5 mg/ml.The configuration providing the lowest CRLB value is then selected.

Material decomposition
The Alvarez-Macowski technique [15] states that a material LAC, referred as , can be expressed as the combination of two selected materials.In this work, the two materials are water and iodine: where () is the LAC to decompose,   and   are the water and iodine densities, respectively; () and   () are the water and iodine attenuation coefficients, respectively.In x-ray transmitted imaging, the attenuation is the integration of the LAC along a x-ray path : where  is the volume location parameter,   and   are the water and iodine lengths, respectively.Given the tube voltage and the source filtration (9.5 mm Al HVL), the iodine K-edge at 33.2 keV can not be exploited and the iodine is not considered here as a K-edge material.In this case it can be selected as one of the two basis materials.

Step-wedge calibration matrix
In this work, the step-wedge phantom is composed of combination of water and iodine lengths.Different effective iodine lengths were obtained by placing different iodine concentrations in vials with the same thickness: 2.1 cm (figure 4)).Thereafter the iodine concentration refers to the iodine length divided by the vial thickness.The step-wedge phantom is characterized by the number of -5 -combinations   .In a step-wedge calibration, different combinations of known material lengths are measured.The resulting number of counts (noted ) for bin  and material combination  is then: where  is the discretized energy index,    is the spectral response per bin combining the source spectrum and the PCD response,    and    the water and iodine lengths of the   ℎ combination, respectively.
The step-wedge calibration can then be expressed as a linear system: where c  and S  are the concatenated counts and spectral response vectors, respectively, and M is the step-wedge calibration matrix.This matrix is represented in figure 3 (left) for a 1 keV energy step and 50 calibration steps: 5 water lengths (2, 4, 6, 8 and 10 cm) and 10 iodine concentrations (0, 1, 2, 3, 4, 5, 7, 8, 9 mg/ml).Figure 3 (right) was constructed to help define the optimal number of material combinations for the step wedge.The x-axis shows the number of combination indices chosen at random from those in figure 3 (left) to form a sample matrix.The y-axis tracks the rank of the matrix, with the ideal rank being equal to the number of rows and the effective rank being the rank of the sample matrix generated.The rank increases linearly with the number of combinations until 10 steps.At a higher number of combinations, the effective rank is no longer linearly increasing and above 30 steps the rank is constant.This means increasing further the number of combination will provide redundant information.In this work, a number of steps equal to 20 was selected because it represents a compromise between the spectral information contained and the acquisition time required.

Step-wedge calibration measurement
The step wedge calibration were made of 20 combinations with 4 water lengths (2.1, 4.2, 7.2, 10.2 cm) and 5 iodine concentrations (0, 1.1, 2.2, 4.4, 8.7 mg/ml), placed in a 2.1 cm thick vial (figure 4).The resulting grid is represented in figure 4.These combinations were used as calibration points for the -6 -models described in the next section.Each combination was acquired as one long acquisition with a time of 200 s to ensure low noise.In addition to the calibration points, five test points were measured to evaluate and compare the various methods (figure 4).These points were also measured during 200 s but with a frame rate of 1 Hz, leading to 200 acquisitions of the same test points (1 s acquisition time).These 200 measurements were used to measure the variance of each method.

Material decomposition methods
In this work, a projection-based approach was investigated [16].In this process the first step is the MD computed pixel per pixel.It aims at obtaining two basis material length sinograms (water and iodine) from the counts measured in the energy bins.This step needs to take into account the spectral distortion described in figure 2. In a second step these sinograms will be reconstructed.
The first method is a direct decomposition method by the evaluation of a polynomial on the measured data [28].
A. Direct polynomial.A second order polynomial directly decomposes the attenuation  into a basis material  length (  ): where   and  0  are the measured counts with and without the object, respectively,  1 and  2 the polynomial coefficients for each bin   is the number of energy bins.These coefficients need to be calibrated.
The other following methods compute the MD with a different scheme.A maximum likelihood algorithm is implemented to find the water and iodine lengths [16].The algorithm requires the definition of a forward model.
-7 -B.Exhaustive model.For a PCD measurement, the exhaustive forward model is: where   (  ;   ) is the estimated count number in bin ,    is the spectral response per bin combining the source spectrum and the PCD response;  is the discretized energy index, and    and    are the water and iodine LAC at energy .The spectral response was initialized with the simulation model displayed in figure 2. It is then adjusted pixel per pixel solving the linear system (2.4) with the maximum likelihood -expectation maximisation (ML-EM) as implemented in Sidky et al., 2005 [19, 20].The algorithm convergence was studied and stopped after 550 iterations to minimize the MD error.For this method a 1 keV sampling between 20 and 130 keV was implemented.With this energy sampling,   has 110 values per bin to be fine-tuned.Models with a reduced number of parameters have been proposed to stabilize the calibration process.
C. Reduced model.The number of incident energies in the method B (eq. (2.6)) can be reduced to a number   [27]: where   ,    and    are the coefficients to be calibrated;   is the number of components to be selected.In this model, the coefficients   ,    and    do not have a physical significance.The advantage of this model is that the equation still reflects the Beer-Lambert law and can be inverted with the same algorithm implementation as method B, utilizing distinct input data for the algorithm.The impact of the number   was studied for   between to two and five.For each case the number of parameters per bin is equal to: 3 ×   .For each component number the coefficients were optimized by maximizing the Poisson likelihood with the Powell algorithm [37] implemented in the Scipy module (version 1.11.3) in Python.This algorithm's stopping criterion is the tolerance level.The convergence of this model and the optimal number of components   were studied.For this purpose, three criteria were computed for different tolerance and component number values.The first one is the normalized mean counts error percentage (CEP) calculated for the calibration points: where   is the number of combinations.
The second and third criteria are the MD error and variance evaluated for the test point measurements (figure 4).D. Linear model.The last investigated method is a specific case of eq.(2.7) where   = 1: In this case the model is linear [22,25,38] and can be written with the attenuation  (eq.(2.5)): where    and    are coefficients calibrated with the same algorithm as method C.This model has two parameters per bin to be calibrated.The convergence of the algorithm in this case will also be investigated with the same methodology as method C.
For the CT acquisitions, the forward models defined in B, C and D are inverted by the Poisson likelihood (PL) maximization [16]: The Gauss-Newton algorithm is implemented to maximize the defined Poisson likelihood (eq.(2.11)) and to compute the basis material lengths.
As the detector panel module is a prototype, defective pixels need to be identified and the data locally interpolated.The step-wedge calibration was also used to identify these defective pixels in the detector panel.For all four methods, the MD was computed on the test points and 10% of pixels with the highest MD error were assessed as defective.

Bias correction
Each method was also investigated for the impact of a bias correction [22,25,39].After a model was calibrated on the step-wedge phantom measurement, it was tested on 20 calibration points.The differences between the resulting value and the theoretical value provided a regular 2D grid of errors.A 2D piecewise linear error map was computed from these points through interpolation.The interpolation performs a triangulation of the data and then a linear barycentric interpolation on each triangle.This is implemented in the Python function LinearNDInterpolator -Scipy version 1.11.3.The bias correction was implemented in each method and results with and without this correction were compared.The defective pixel masks were adjusted for each model with the bias correction.

Phantom CT acquisitions
Two phantoms were scanned to compare the four methods.The CT acquisitions were performed with the same source parameters and detector thresholds as the step wedge calibration measurements.Each phantom acquisition contained 720 projections over 360°that were decomposed into water and iodine sinograms.This projection number ensures sufficient geometric conditioning to not interfere with quantitative results.These sinograms were then reconstructed with a filtered backprojection (FBP) algorithm implemented in the Reconstruction Toolkit [40].The voxel size was set to 0.15 × 0.15 × 0.2 mm and three slices were reconstructed.
The first phantom (figure 5) is a 100 mm diameter cylinder made of high density PE (HDPE).The phantom held four tubes filled with water, seven tubes filled with iodine (Iomeprol, Iomeron 400, Bracco Imaging, Milan, Italy) at different concentrations, and seven tubes filled a calcium perchlorate tetra-hydrate (CaCO 4 -Sigma-Aldrich, Saint Louis, U.S.A.) compound diluted at different concentrations.For the iodine, three tubes were filled at 1.1 mg/ml, two tubes at 2.2 mg/ml, one tube at 4.4 mg/ml, and one at 8.7 mg/ml.For the calcium compound, three tubes were diluted at 52 mg/ml, two tubes at 104 mg/ml, one tube at 208 mg/ml, and one at 416 mg/ml.The iodine error was assessed as the absolute difference between the tube value in the reconstructed iodine map and the true value.
The second phantom presented a cylindrical symmetry (figure 5) and its center was placed perfectly at the rotational axis position.It was made of only two materials: a 100 mm diameter HDPE -9 - cylinder with a 50 mm diameter hole in the center.The hole was filled with iodine at 2 mg/ml.This phantom was scanned two consecutive times with identical acquisition parameters.The water/iodine MD of this phantom was computed with all four methods and both acquisitions.Then the consecutive basis material maps were subtracted one from another to obtain white noise images.The noise power spectrum (NPS) [41] was evaluated on these images for each material map and the noise level was defined as the NPS area.This technique enables to separate the noise measurement from the ring artifact intensity.The first acquisition was studied further to measure the ring deviation.For this purpose, profiles from the center to the insert edge were radially averaged.Then a second order polynomial was subtracted from the profile to obtain a white noise profile.The ring artifact deviation was then defined as the standard deviation along this profile [33].As this deviation was impacted by the image noise, the ring artifact deviation was normalized by the noise value.

Impact of the number of energy bins
The method which performed the best among the four methods presented in section 2.7 was selected to study the impact of the number of energy bins on CT image quality.The thresholds for each number of bins were selected according to the methodology described in section 2.3.The results were compared on the test calibration points for each number of bins and only the phantom with tubes was scanned.The iodine concentration error was evaluated and the noise was defined as the standard deviation in the background of both basis material maps.

Lamb rib CT acquisitions
In order to illustrate the water/iodine MD for clinical applications, an ex-vivo ovine rib was scanned with the same acquisition parameters as the phantoms.The water/iodine MD was computed with the best performing method.Various spectral images used in clinical practice were then computed.The first one was the iodine concentration that was directly obtained from the iodine density maps.Virtual mono-energetic images (VMIs) were computed at 40 and 150 keV with the linear combination of the two material maps (eq.(2.1)).Virtual non-contrast (VNCo) image was computed as the normalization of the water density map in Hounsfield units (HU).Additionally an overlay of the VNCo with the iodine map is presented.Finally, a virtual non-calcium (VNCa) image was computed.For this image, -10 -the bone was segmented by thresholding the iodine above 2.5 mg/ml and the water below 1000 mg/ml.Then a VMI was computed (eq.(2.1)), but instead of adding the iodine signal, it was subtracted in the bone area.All the generated images correspond to different clinical applications.

Thresholds calculation
The results of the CRLB computation for each number of bins is presented in figure 6 left).These results are displayed as a function of the highest threshold to enable comparison across different numbers of energy bins.As the threshold sampling remains constant within each configuration, the highest threshold uniquely characterizes each setup.The sampling rate can be computed by subtracting 24 from the highest threshold and dividing by the number of bins.In figure 6 left), the CRLB has been normalized by dividing it by the minimum value observed in all energy bins.For each number of energy bins, an optimal value exists and the bin configuration corresponding to this value was selected.
The optimal threshold configuration corresponding to each number of bins is displayed in table 1.Interestingly, new threshold configurations with a smaller number of energy bins can be measured as a re-binning of another configurations.For example, the four energy bins configuration can be measured by binning the seven energy bin configuration.The second result of the threshold optimization was that the optimal CRLB value decreased with the number of bins.This results has already been demonstrated by Alvarez, 2011 [22], but with an idealized detector response.In order to quantify the difference between the number of energy bins, the optimal CRLB values are displayed in figure 6 right).The largest CRLB increase occurs between the number of bins   = 3 and   = 4.The ratio between the two energy bin and the seven energy bin configurations is around 1.13 (depending on the material).This ratio can be superior to the ratio between two threshold configurations for the same bin number (figure 6 left) and emphasizes the importance of the threshold choice on the MD noise for a given number of energy bins.

Model parameters
For the models with a reduced number of parameters (methods C and D), two parameters could impact the material decomposition: the number of components and the stopping criterion for the algorithm optimizing the coefficients.The results of this parameter study are displayed in figure 7. The number of components   = 1 represents the linear model (method D).This model converged faster than the models with a number of components   > 1.For these models, the MD error difference between -11 - the number of components was negligible when the convergence was reached.This demonstrates a difference compared with the results observed by Mechlem et al., 2017 [27], who observed that a number of components   = 4 led to a significantly more accurate model.In this work, the number of components   = 2 was selected.For this case, lowering the tolerance level below 10 −7 did not improve the MD error, but did increase the variance for the iodine material.In order to avoid over-fitting of the model, a tolerance level equal to 10 −7 was considered optimal.For the linear model (  = 1), a tolerance level equal to 10 −6 was selected as the model converged for this tolerance level.

Model comparison
The evaluation of the criteria for each of the four methods without and with bias correction is presented in tables 2 and 3, respectively.The quantitative results demonstrate that the polynomial method A provided the best spectral accuracy on the calibration test points and in the CT images.However, it suffers from a high noise level.Method B did not reach satisfactory spectral accuracy which led to incorrect material quantification.Method C and D had similar results with an iodine quantification error below 0.4 mg/ml in CT images.Method D presents a slightly lower noise level than method C but method C reduced ring deviation compared to method D.
-12 - The results with the bias correction implemented are presented in table 3 with a color code to illustrate if the correction improved the result.The correction had almost no impact on the polynomial model.It significantly improved the exhaustive model spectral accuracy but this led to increased variance and noise properties.Even with the bias correction, the exhaustive model presents severe artifacts and a high iodine concentration error of 0.85 mg/ml.Finally, the bias correction led to calibration improvement for method C and D, but minimal difference in the CT images.
The water/iodine CT images for each method computed with the bias correction are presented in figure 8.It is visible that method A provided images with a higher noise level.The exhaustive method B presented extensive artifacts and the heterogeneity between the detector sub modules is evident.Method C and D generated images with similar properties, but the linear Method D caused more ring artifacts.The water/iodine images corresponding to the cylindrically symmetric phantom case are present in the supplemental material.
Based on these results, method C with 2 components had the best trade-off between satisfying iodine quantification accuracy and noise/ring artifact level in CT images.That is why this method was considered the best one overall and further evaluated.The bias correction for this method had a minor but positive effect and was therefore implemented in the next results.Moreover, once it is calibrated, this bias correction has a very low computational cost.

Impact of the number of bins
The impact of the number of bins study for method C are displayed in table 4. The results for the calibration test point accuracy are the same for all studied scenarios.The iodine concentration error increased for the number of bins equal to two and three.In order to illustrate the noise increase, the noise values are displayed in figure 9.The curves are not perfectly monotonic but the noise is decreasing overall with the number of bins and the smallest value is measured for the seven energy bins configuration.The noise value in the CT images are in accordance with the measured calibration -14 -variance values.A noise increase between 10% (water) and 15% (iodine) is observed between the cases with two and seven energy bins.Therefore the ex-vivo samples were reconstructed from a seven energy bins measurement, which presented the best performance.

Lamb rib
The spectral images of the lamb rib injected with iodine are displayed in figure 10.The iodine signal was successfully subtracted from the conventional image (map A) to produce the VNCo (map E).The iodine was injected at 30 mg/ml and the mean measured concentration was 8.0 mg/ml in the iodine map (map C).The iodine attenuation properties are well represented with the contrast agent appearing very bright in the VMI at 40 keV (map B) but diminished at 150 keV (map F).As in current clinical CT images, the bone signal is divided between the iodine and water material (maps C and G).This enables the identification and removal of the calcium portion of the object to generate the VNCa (map H).
The lamb rib images prior to iodine injection are provided as supplemental material.

Discussion
The exhaustive model calibration failed to provide satisfactory results with large spectral errors, even with bias correction implemented.This method strongly relied on the input of the simulated detector model.The detector model has been refined by means of multiple polyenergetic spectrum -15 - acquisitions such as those presented in figure 2, but this has not led to model improvements.The model could be improved through mono-energetic measurements such as those from a synchrotron source.Unfortunately accessing such an x-ray source and refining the model is very complex.Instead, a model with a reduced number of parameters has been successfully calibrated.The step wedge phantom design made from water and iodine enabled an iodine concentration accuracy below 0.4 mg/ml in CT images, similar to the error values reported in photon-counting preclinical or clinical systems with equivalent source filtration [3,5].The imaging system was operated at low flux with narrow collimation in order to reduce the pile-up effect and the x-ray scatter.In this situation, the< exhaustive model can be modeled by equation (2.6).The method C and D with a reduced number of parameters are also subject to these conditions and could be negatively impacted by a higher incoming flux or scattered x-rays.The detector behavior depends on the incoming flux and the model should be adjusted to this flux level [42].The scatter could be reduced by the utilization of an anti-scatter grid (ASG) that is not available on our system.The dependence of the model to these conditions was not evaluated.This step-wedge calibration can be performed directly in a CT clinical examination room.However, a different calibration needs to be acquired and processed when the source parameters are modified (voltage or current).
The presented method C shares the same equations as the exhaustive method B. The difference lies in the inputs with a set of learned parameters that replace the physical significance of the basis material attenuation coefficients and spectral response.The number of components is a parameter of this model, and it has been demonstrated that the ideal number is   = 2, leading to the calibration of six values per bin.The limitation of this model is the restriction of the basis material choices.As the calibrated parameters do not have a physical significance, they are only valid for a water/iodine MD.However, the implementation of bias correction with any model would restrict the basis material choice.Moreover, this material pair is relevant for clinical applications as the iodine maps enable direct contrast agent concentration quantification and the computation of the VNCo map [2].Our innovative step wedge phantom facilitates direct iodine concentration measurements in CT images by computing the iodine basis material density map.The lamb rib images (figure 10) thoroughly demonstrated the -16 -relevance of this MD choice, enabling the segmentation between bone and iodine and the generation of VNCo and VNCa maps.The water/iodine MD does not enable the segmentation of soft tissue (separating muscle from adipose tissue).However, this can be performed in the conventional image where the adipose tissue is characterized by negative CT numbers as presented in figure 11.The water/iodine maps can then provide material quantification of these tissues as presented in table 5.It is noted that the adipose tissue has a negative spectral footprint in the iodine map.Methods A and D presented interesting results.The polynomial method A is a relevant gold standard comparison for material quantification.For this method, the bias correction had almost no impact on the quantitative results.Moreover, this model has the lowest computation time as the MD is performed only by evaluating the calibrating polynomial on the measured data.Its drawback is a high noise level.This method could be used as an initial state of a more sophisticated reconstruction technique such as one-step MD algorithms or machine learning algorithms that perform efficient noise reduction.The linear method D (or one component model) presented results close to method C, only suffering from more pronounced ring artifacts, visible in figure 8.A model reducing ring artifacts was prioritized here.However, the ring artifacts could be dealt with by implementing a ring artifact correction method.The detector panel is also a prototype version and future detector with a better pixel response homogeneity could result in a reduction of ring artifact.Compared to method C, the linear model has two minor advantages: a slightly lower noise level and the use of a linear system optimization.
-17 -The bias correction significantly improved method B, which suffered form high quantitative deviations.However, it was not sufficient to reach a sufficient MD error.Moreover, a high level of bias correction leads to an increase of image noise.The bias correction enabled small improvements of method C and D, but did not result in improvement of iodine concentration measured or ring artifact reduction in CT images.This correction is performed with an interpolation of the step wedge calibration points.In this case, increasing the number of calibration combinations could refine the grid sampling and further improve the bias correction.This requires a significantly longer calibration time and was not studied in this work.The total acquisition time for the studied number of combinations was around 80 minutes.The optimization of this time was not studied.This would depend on the selected incoming flux as a sufficient number of counts is required to avoid a noise bias in the calibration point.The long term stability of PCD under clinical constraints is still under investigation.It remains uncertain at what frequency the calibration would be required to be updated, even for the same source parameters.This issue will need to be further monitored.
The study investigating the number of bins has not demonstrated a reduced spectral accuracy at lower bin numbers for calibration or CT images, but it did demonstrate a higher noise level.With optimized threshold configurations, the noise with two energy bins is 10% (in water map) and 15% (in iodine map) higher than with seven energy bins.These values are in accordance with the simulation (figure 6 right)) and past results [5,22].The limitation of this study is the non fully realistic detector response model used for the CRLB computation and it could explain the non monotonic aspect of the curve presented in figure 9.As displayed in figure 2 and demonstrated by the exhaustive method results, this method is imperfect and can not fully match the detector physics on a pixel-by-pixel basis.However, it provides a good first order estimation of the model and it was shown to result in an accurate threshold calibration [33].As presented in figure 6 right), the minimum CRLB for the cases   ≥ 3 is very close and the cases   = 4 and   = 6 presented in figure 6 left) were not perfectly quadratic.This demonstrates that the optimal threshold configuration computation for each number of bins is sensitive.Moreover, only configurations with a regular spacing were investigated.Future work will investigate methods to refine the optimal threshold configuration.

Conclusions
We have performed accurate iodine quantification in CT images with a bench top PCD-CT system.For this purpose, a reduced model that does not require a detector response was calibrated with step wedge phantom measurements.The calibration include a bias correction that slightly improved the results.The water/iodine material decomposition was tested on a lamb rib injected with iodine to illustrate the clinical relevance of this basis material choice.It was demonstrated that compared to seven energy bins, two energy bins provide equivalent material decomposition accuracy in CT images but a 10% noise increase.

Figure 1 .
Figure 1.Top and side view of the prototype system.

Figure 2 .
Figure 2. Measurements of a 120 kVp polyenergetic spectrum for 136 adjacent pixels of the detector, each represented by a different color (excluding pixels at the edge of a sub-module).Std: standard deviation.

Figure 3 .
Figure 3. Left: the step-wedge calibration matrix M represented for 50 combinations: 5 water lengths and 10 iodine concentrations.Right: rank of the matrix  as a number of combinations.The combinations were randomly selected in the matrix displayed on the left.

Figure 4 .
Figure 4. Left: calibration points scheme.Right: calibration grid presenting the 20 calibration combinations and the 5 testing points.The x values are the iodine concentrations in the 2.1 cm thick vial.The red arrows indicate the variance that will be assessed at the test points.

Figure 5 .
Figure 5. Left: conventional image (Hounsfield Unit -HU) of the HDPE phantom with an iodine insert.Right: conventional image (HU) of the HDPE phantom with tubes of water, iodine and calcium perchlorate tetra-hydrate.

Figure 6 .
Figure 6.Left: CRLB for the water (dashed lines) and iodine (continuous lines) length estimations as a function of the highest threshold.Right: optimal CRLB values per bin normalized by the minimal value.

Figure 7 .
Figure 7. Parameter studies for the method C and D, using reduced models.

Figure 8 .
Figure 8. Water/Iodine density maps for the four tested methods and the bias correction activated.The green and yellow circles point to the tube with the highest concentration of  4 and iodine, respectively.

Figure 9 .
Figure 9. Noise values measured in the phantom #1 background normalized by the value for the case   = 7.

Figure 10 .
Figure 10.Spectral images for a lamb rib injected with iodine.

Figure 11 .
Figure 11.Lamb rib conventional image (left) and the four materials segmentation masks.

Table 1 .
Optimal threshold configuration per number of energy bins.

Table 2 .
Comparison of the four methods without bias correction.The bold figure represents the best value per criterion.

Table 3 .
Quantitative results with bias correction.The bold figure represents the best value.The cell is colored in red, yellow, or green if the correction worsens, did not change, or improved compared to table 2, respectively.

Table 4 .
Quantitative results for the different number of bins.The bold figure represents the best value.

Table 5 .
MD quantification for each material measured in the areas presented in the left image.