An Image-Based Prior Knowledge-Free Approach for a Multi-Material Decomposition in Photon-Counting Computed Tomography

Photon-counting CT systems generally allow for acquiring multiple spectral datasets and thus for decomposing CT images into multiple materials. We introduce a prior knowledge-free deterministic material decomposition approach for quantifying three material concentrations on a commercial photon-counting CT system based on a single CT scan. We acquired two phantom measurement series: one to calibrate and one to test the algorithm. For evaluation, we used an anthropomorphic abdominal phantom with inserts of either aqueous iodine solution, aqueous tungsten solution, or water. Material CT numbers were predicted based on a polynomial in the following parameters: Water-equivalent object diameter, object center-to-isocenter distance, voxel-to-isocenter distance, voxel-to-object center distance, and X-ray tube current. The material decomposition was performed as a generalized least-squares estimation. The algorithm provided material maps of iodine, tungsten, and water with average estimation errors of 4% in the contrast agent maps and 1% in the water map with respect to the material concentrations in the inserts. The contrast-to-noise ratio in the iodine and tungsten map was 36% and 16% compared to the noise-minimal threshold image. We were able to decompose four spectral images into iodine, tungsten, and water.

PCD CT systems, for the first time, enable the acquisition of more than two spectral images in a single scan and thus the decomposition into more than two base materials, hereafter called multi-material decomposition (MMD).In combination with the development of novel contrast agents (CAs) based on materials with high atomic number and thus high K-absorption edges (hereafter called K-edge CAs), this opens up new scanning techniques [20][21][22].As an example, a patient with a known history of hemochromatosis could be injected with a K-edge CA, which would allow an accurate iron quantification via MMD into iron, liver tissue, and the K-edge CA.To examine the lung ventilation, a patient could breathe in xenon and be injected with a K-edge CA to allow MMD into xenon, the K-edge CA, and water, which indicates soft tissue.
CT-based MMD has already been investigated [21,[23][24][25][26][27][28][29][30][31][32], including in prototype CT systems and simulations which provided at least three spectral datasets.It was also assessed in clinically used CT systems, which provided at most two spectral datasets.Consequently, the clinical system required repeated measurements or boundary conditions resulting from prior knowledge of the scanned object or assumptions as volume conservation or mass conservation for instance [23,33].
The new commercial CT systems involving PCD overcome these limitations.These systems allow a one-shot acquisition of four spectral images and thus prior knowledge-free MMD by a single scan.However, PCD CT systems with challenging effects such as pile-up and charge sharing require calibration of material CT numbers for an accurate material quantification [16,34].
In this work, we present and test an image-based prior knowledge-free method to perform MMD on the new PCD CT systems involving models to calibrate material CT numbers.We first describe the general method, applicable to a set of arbitrary materials with distinct spectral behavior and an equal or larger set of spectral images.To demonstrate the method's feasibility, we applied the approach to scans of an abdominal phantom using different phantom diameters, tube currents, and positions.Our material system included iodine, tungsten, and water.We acquired four spectral images per scan, which are energy threshold images in the used CT system.
We investigated the deviations of the material concentration estimation (CE) from the true value and the statistical uncertainty of the MMD results, analyzed by comparing the contrast-to-noise ratio (CNR) in the algorithm input and output.
In contrast to previous works [23][24][25][26][27][28][29][30][31][32][33]35,36], we used commercial hardware allowing one-shot acquisition of more than two spectral datasets.This is the first time a prior knowledge-free calibration and decomposition approach involving more than two materials was tested on a clinical PCD CT system using a four-threshold image acquisition protocol.In the literature, we have not found similarly detailed models as those we used to calibrate material CT numbers.We moreover created a model for the image noise covariance matrix (INCM) to reduce the random noise in the material maps.

Materials and Methods
Our approach is introduced in Section 2.1.We first present the method structure, outlined in Figure 1, and then discuss the method details involving model functions for material CT numbers and the INCM.All models are created analytically and do not require manually fine-tuned initial conditions to converge to a good optimum.The proposed model functions are further discussed in Appendix A. Two measurement series, one to calibrate the models and one to test the calibrated Two measurement series, one to calibrate the models and one to test the calibrated models, will be described in Section 2.2.We finally discuss the analysis of the precision and accuracy of the material images resulting from applying our approach to the test data as described in Section 2.3.The customized code for model creation, material decomposition, and image analysis was developed in MATLAB (Version R2021b, The MathWorks Inc., Natick, MA, USA).

General Method
The general method has two parts: First, the creation of CT number models of all materials of interest in the N spectral images series provided by the CT system and of a model of the INCM.Second, the decomposition of N input images into M ≤ N materials, called base materials (BMs), resulting in a vector → c ∈ R M of estimated BM concentrations for each image voxel.
• → y ∈ R N is the vector of CT numbers of one image voxel at fixed row and column in the different spectral images; • → x = (x 1 , x 2 , . . . are the parameters on which the models depend; Given the calibrated models and N spectral images to decompose, the material decomposition for a voxel with parameters → x and CT numbers → y is prepared as follows: → α with f the model of the r-th BM CT number in the p-th spectral image (p = 1, . . ., N), scaled to give the expected CT numbers at a concentration of 1 mg/mL for easy interpretation of the decomposition results.If the r-th BM is water, S pr is set to 1 for all p to ensure a linear water CE in the true water concentration for voxels containing at most water, interpreted as explained at the end of this section.

4.
The user then estimates the material coefficients → ρ ∈ R M using the method of gener- alized LS.This means solving the linear system of equations 5.
The result → ρ is the best linear unbiased estimator for the overdetermined system c is as follows: if the r-th material is a CA, ρ r equals c r in units of mg/mL.Note that the user could fill the r-th column in S with the expected CT numbers at any non-zero concentration q mg/mL.The interpretation is then that ρ r •q equals c r in units of mg/mL.
In the case of water, c r = (1000 + ρ r ) mg/mL, which means that ρ r = 0 corresponds to a water-filled image voxel plus possible CA and ρ r = −1000 means water absence.
By repeating this decomposition for all image voxels of the N input images, M material maps displaying material concentrations are generated.
Note that the choice of distinguishable materials is limited by their spectral separation.Since in CT, the absorption behavior of human body materials is dominated by the photoelectric effect and Compton scattering [34], discrimination of more than two materials requires all further materials to exhibit a unique K-edge absorption in the diagnostic energy window.

Model Parameters
All our models used an individual subset of the following calibration parameters:

•
Object water-equivalent diameter (WED) [mm] (x 1 ); • X-ray tube current [mA] (x 2 ); • Voxel distance to the isocenter [mm] (x 3 ); • Voxel distance to the object center [mm] (x 4 ); The WED is slice-wise computed from the image CT numbers.It equals the diameter of a disk of water providing the same area-integrated attenuation as the measured object [37].Each image voxel has associated distance measures x 3 −x 5 visualized in Figure 2a.The object off-centering is the distance between the isocenter and the center of mass of the measured object.

Material Model Creation
The material models consider several parameters which affect the measured CT number in a PCD CT system exemplarily shown in Figure 3.With an increasing object diameter, the X-ray spectrum becomes harder.This leads to a change in CA attenuation and hence the measured CT number.The PCD's sensitivity to the photon flux due to pile-up and charge sharing explains the non-zero correlation between CT number and tube current and between CT number and distance to the isocenter due to form filters on top of the X-ray source.

Material Model Creation
The material models consider several parameters which affect the measured CT number in a PCD CT system exemplarily shown in Figure 3.With an increasing object diameter, the X-ray spectrum becomes harder.This leads to a change in CA attenuation and hence the measured CT number.The PCD's sensitivity to the photon flux due to pile-up and charge sharing explains the non-zero correlation between CT number and tube current and between CT number and distance to the isocenter due to form filters on top of the X-ray source.

Material Model Creation
The material models consider several parameters which affect the measured CT number in a PCD CT system exemplarily shown in Figure 3.With an increasing object diameter, the X-ray spectrum becomes harder.This leads to a change in CA attenuation and hence the measured CT number.The PCD's sensitivity to the photon flux due to pile-up and charge sharing explains the non-zero correlation between CT number and tube current and between CT number and distance to the isocenter due to form filters on top of the X-ray source.We recommend modeling CA CT numbers in a spectral image defined by the following function: with constant numbers {α i ∈ R : i = 1, . . ., 10} left to be determined using G f .If water- filled regions do not all receive a CT number of 0 HU from the CT system, an intermediate computation, called water correction, is required.This consists of replacing the ground truth CT number y k by y ′ k with for each k = 1, . . ., K where w is the CT number of water measured at the same configuration → x k [38].According to the Hounsfield scale which is used for image-based material decomposition, water-filled regions should receive a CT number of 0 HU from the CT system.Regions with aqueous CA solutions should receive a CT number that is proportional to their concentration, at least in the thin-absorber limit [39].This is not everywhere the case in the images of the system we used.Applying Equation ( 5), we come closer to the preconditions while leaving material-free regions at −1000 HU.
Optimal model coefficients are then determined through an ordinary LS fit using the chosen model function f and the ground truth data {( . ., K}.We moreover recommend modeling water CT numbers using the function (6) and G g to calibrate the set {α i ∈ R : i = 1, . . ., 21}.The procedure is analogous to the CA models except for the water correction of the data points.Even if water is not chosen as a BM, the water model is created to allow water correction of input images for the decomposition algorithm.This makes CA CE linear in the voxel CT numbers.In particular, if one uses aqueous solutions as BMs, water-filled regions should be evaluated as having zero concentration.Aqueous solutions are a convenient choice for material decompositions in human soft tissue.For a voxel-wise water CT number prediction at arbitrary object size, position, and tube current, we need a water model in each threshold.
If water-filled regions receive a CT number of 0 HU from the CT system, the water model and any water correction is obsolete.

Image Noise Covariance Model Creation
The INCM model H serves for optimally solving Equation (2).For each combination of spectral indices p and q (p, q = 1, . . ., N), the matrix element H pq is the modeled image noise covariance between the p-th and the q-th spectral image.We split the model H pq into with I pq (p, q = 1, . . ., N) a model of the noise's correlation between the p-th and the q-th spectral image and j p (p = 1, . . ., N) a model of the noise's standard deviation in the p-th spectral image.Likewise, we extract correlations and standard deviations from the measured INCMs in G H to acquire ground truth for the models I pq and j p .
We recommend assuming the correlation not to depend on any of the parameters x 1 -x 5 and thus set I pq to the arithmetic mean of the available ground truth data of I pq .We moreover recommend modeling j p using the following function: Because the noise variance σ p 2 is proportional to the reciprocal of the tube current x 2 and increases exponentially with the absorption distance, and thus with the WED x 1 [40].See Appendix A for details on how j p was calibrated.

Experimental Setup and Image Acquisition
To apply the method proposed in Section 2.1, we performed two CT measurement series: calibration and test data acquisition.The calibration series provides ground truth data for the model calibration.The test data acquisition serves for inspecting a threematerial decomposition.

Material Description
We measured cylindrical WEP phantoms with heights of 100 mm and diameters of 150, 200, 250, 300, and 350 mm (QRM, Möhrendorf, Germany) for the calibration, cf. Figure 2c.The diameters were chosen to represent most clinical applications with different patients and body parts under investigation.These five phantoms were placed on the scanner table with the rotational axis pointing in the axial direction.
To acquire test data, we measured an abdominal phantom (QRM, Möhrendorf, Germany) with a height of 165 mm and an elliptical cross-section with a width of 300 mm and depth of 200 mm once without and once with an obese extension ring (QRM, Möhrendorf, Germany), cf. Figure 4, of width 400 mm and depth 300 mm (WED = 273 mm and 363 mm).The two configurations represent typical average adult patient sizes for abdominal scans.These two phantoms were positioned on the scanner table with the rotational axis of the central water-equivalent cylinder, visible in Figure 4, pointing in the axial direction.
moreover recommend modeling  using the following function: Because the noise variance  is proportional to the reciprocal of th and increases exponentially with the absorption distance, and thus with See Appendix A for details on how  was calibrated.

Experimental Setup and Image Acquisition
To apply the method proposed in Section 2.1, we performed two series: calibration and test data acquisition.The calibration series provi data for the model calibration.The test data acquisition serves for inspe terial decomposition.

Material Description
We measured cylindrical WEP phantoms with heights of 100 mm 150, 200, 250, 300, and 350 mm (QRM, Möhrendorf, Germany) for the cali 2c.The diameters were chosen to represent most clinical applications w tients and body parts under investigation.These five phantoms were pl ner table with the rotational axis pointing in the axial direction.
To acquire test data, we measured an abdominal phantom (QRM, M many) with a height of 165 mm and an elliptical cross-section with a wid depth of 200 mm once without and once with an obese extension ring (Q Germany), cf. Figure 4, of width 400 mm and depth 300 mm (WED = mm).The two configurations represent typical average adult patient siz scans.These two phantoms were positioned on the scanner table with t of the central water-equivalent cylinder, visible in Figure 4, pointing in t   In the abdominal phantom, we fill three drill holes, each with one insert: One in the liver-equivalent region, one in the spleen-equivalent region, and one in the soft tissue-equivalent region.The inserts were cyclically swapped between the scans, mentioned in Table 1, to measure each one in all the different regions.The central hole of diameter 100 mm is filled by water-equivalent plastics.All seven phantoms were positioned as described to achieve axially homogeneous images.Note that all phantoms contained drill holes for placement of syringes.
Our set of M = 3 base materials consists of the following: Iodine dissolved in water (c = 10 mg/mL), hereafter called iodine; • Tungsten dissolved in water (c = 10 mg/mL), hereafter called tungsten; • Water (c = 1000 mg/mL).
By filling 20 mm diameter syringes each with one of these liquids, we prepared three inserts for the phantoms.We chose tungsten to satisfy the requirement of one material clearly interacting with incident X-rays in CT by more than the photoelectric and Compton effect.The concentration of CA in the solutions was chosen to be among the clinically common iodine concentrations below 25 mg/mL [41][42][43].

General Scan Setup
All data were acquired by a dual-source photon-counting scanner (NAEOTOM Alpha, Siemens Healthineers, Forchheim, Germany).This scanner generates up to four energy threshold images, denoted T1, T2, T3, and T4.A threshold image is composed of all photons which deposit energy above the threshold energy in the detector pixel.
We used several fixed scan and reconstruction parameters throughout the experiments as shown in Table 2.In a four-threshold acquisition, the user can perform single-source scans at a tube voltage of 140 kV, rotation time of 0.25 s, and axial collimation of 96 × 0.4 mm.We chose the threshold energy quadruple (20, 55, 72, 90) keV expected to provide the best spectral separation between tungsten, with K-edge energy of 69.5 keV, and the other BMs [44].The reconstruction algorithm was a weighted filtered backprojection [45].Neither iterative reconstruction nor post-reconstruction filtering were applied.The chosen slice thickness and increment of 2 mm is used for several examinations of the head and spine, for example [43].We chose a large field of view to allow retrospective WED calculation.

Calibration
The calibration-specific scan setups are shown in Table 3 and an exemplary phantom in Figure 2c.We performed scans with the listed setups once for each of the three inserts located in the drill hole 30 mm away from the phantom center.We matched the tube currents to the phantom sizes to achieve the same five dose levels for each phantom.After extracting the ground truth for all models as described in the following paragraphs, we applied the recommended procedure described in Section 2.1.2.To extract ground truth for the material models, we averaged the CT numbers from cylindrical volumes with a diameter of 12 mm over the seven central slices in the reconstructed spectral images.To each of the ROIs, we assigned the distances x 3 − x 5 in relation to the circle center.
For the CA models, we averaged the CT numbers in the ROI within the insert of the scans with CA insert as described in Figure 2 and applied the water correction using averaged CT numbers of the corresponding ROI of scans with the water insert.This yields the threshold-wise data {( → x k ; y k ) : k = 1, . . ., K} for the LS optimization.For the water model, we determined the ground truth by averaging the CT numbers from the ROIs in WEP as described in Figure 2.
For the INCM compounds, the ground truth was extracted from the image voxels in the reconstructed images with a distance between 25 and 35 mm to the phantom center and 10 mm to the insert in each threshold image.We did not use all image voxels because a voxel sample with an inhomogeneous CT number background adds low-frequency components to the noise measurement, possibly biasing the measurement [46].With an incorrect estimated ratio of the noise between the thresholds, the solution of Equation ( 2) would no longer be the best estimator for the overdetermined LSE in Equation (3) among the linear unbiased ones and would thus suffer from an increased expectation estimation error [47].

Test Data Acquisition
Another measurement series with the parameters described in Table 1 and phantom shown in Figure 4 was performed to test the approach.We varied the phantom off-centering by using different table heights to simulate different levels of non-optimal object positioning on the table.
To investigate systematic errors, we repeated each setup according to Table 1.We matched the repetition numbers to the phantom sizes and tube currents to receive equal noise levels in the evaluation of each setup.

Statistical Analysis
After decomposing the seven central slices of all reconstructed series of the test data into iodine, tungsten, and water using the calibrated models, we analyzed the precision and accuracy of the resulting material maps.

Contrast-to-Noise Ratio
We evaluated the phantom ROIs enclosing inserts for analyzing the CE precision measured by the CNR as follows.
Due to the decomposition process described in Section 2.1, the material map noise was caused by the quantum noise in the threshold images and its propagation into the decomposed images.To understand the influence of the MMD on the material image noise level, we investigated the change in CNR, which typically decreases for materialdecomposed images.We defined the input CNR based on the T1 image because it is created by all incident photons on the detector.Thus, the input CNR a in in one ROI with average CT number quadruple → y = (y 1 , y 2 , y 3 , y 4 ) in one image slice reads where σ in is the image noise standard deviation in the T1 image measured in the voxels of the ROI.The output CNR a r, out in the same ROI of the r-th material map with average CE ρ r reads a r, out = ρ r √ σ r, out (10) where σ r, out is the image noise standard deviation measured analogously to σ in .The relative change in CNR, denoted CNR out /CNR in , of the r-th material was then defined as the ratio of the output and r-independent input CNR: Due to vanishing input CNR in CA-free ROIs and output CNR close to zero in ROIs with one CA in the map of a different BM, CNR out /CNR in according to Equation ( 11) is a meaningful number and was computed uniquely in the following cases:

•
Iodine CE in the iodine insert; • Tungsten CE in the tungsten insert.
We finished by taking the arithmetic mean over the central slices and possible scan repetitions to mitigate the contribution of quantum noise to the CNR quotient.Note that due to varying photon statistics, the noise in a in and a r,out changes with the tube current.As the material decomposition is a linear operation on the input image, cf.Equation (2), we do not expect a dependence of CNR out /CNR in in the tube current.

Estimation Accuracy
To investigate the accuracy of the concentrations resulting from the MMD, we compared them with the true concentrations.We expected a water concentration of 1000 mg/mL throughout the phantom.We expected an iodine concentration of 10 mg/mL in the iodine insert and 0 mg/mL otherwise.We expected a tungsten concentration of 10 mg/mL in the tungsten insert and 0 mg/mL otherwise.
We averaged the CT numbers within each of the three inserts.The resulting mean values of the seven slices in all repetitions were then averaged to mitigate the contribution of quantum noise to the CE errors.We then took the difference between the averaged means and the expected value.

Results
Exemplary material maps from two test data setups are visualized in Figure 5.The CA maps mainly show the insert containing the respective material.For each phantom, the tungsten map is noisier than the iodine map.The maps show the WEP and the water insert with a concentration of approximately 1000 mg/mL, lower concentration in the obese extension ring, and higher concentration in the CA inserts and abdominal material.The material maps for the phantom with extension ring are noisier than the ones for the phantom without.

Contrast-to-Noise Ratio
Figure 6 shows CNR out /CNR in for each CA.We found CNR out /CNR in in the iodine map within the range of ±3% around the arithmetic mean at 36% and in the tungsten map in the range of ±1% around the arithmetic mean at 16%.For each CA, the values tend to be higher for the larger phantom diameter.Figure 6.CNR out /CNR in between each contrast agent map and the T1 image, of the scans listed inTable 1.Each box displays the median of all contributing data points and is bounded by the lower and upper quartiles.The whiskers extend to the extreme values except for outliers, which are defined by being more than 1.5 interquartile ranges away from the top or bottom of the box.

Estimation Accuracy
Figure 7 shows the deviations of the CEs from the true values for the nine combinations of material maps (iodine, tungsten, and water) and inserts (iodine, tungsten, and

Contrast-to-Noise Ratio
Figure 6 shows CNR out /CNR in for each CA.We found CNR out /CNR in in the iodine map within the range of ±3% around the arithmetic mean at 36% and in the tungsten map in the range of ±1% around the arithmetic mean at 16%.For each CA, the values tend to be higher for the larger phantom diameter.

Contrast-to-Noise Ratio
Figure 6 shows CNR out /CNR in for each CA.We found CNR out /CNR in in the iodine map within the range of ±3% around the arithmetic mean at 36% and in the tungsten map in the range of ±1% around the arithmetic mean at 16%.For each CA, the values tend to be higher for the larger phantom diameter.Figure 6.CNR out /CNR in between each contrast agent map and the T1 image, of the scans listed inTable 1.Each box displays the median of all contributing data points and is bounded by the lower and upper quartiles.The whiskers extend to the extreme values except for outliers, which are defined by being more than 1.5 interquartile ranges away from the top or bottom of the box.

Estimation Accuracy
Figure 7 shows the deviations of the CEs from the true values for the nine combinations of material maps (iodine, tungsten, and water) and inserts (iodine, tungsten, and water).We found average false estimates of 0.4 mg/mL in the iodine and in the tungsten Figure 6.CNR out /CNR in between each contrast agent map and the T1 image, of the scans listed in Table 1.Each box displays the median of all contributing data points and is bounded by the lower and upper quartiles.The whiskers extend to the extreme values except for outliers, which are defined by being more than 1.5 interquartile ranges away from the top or bottom of the box.

Estimation Accuracy
Figure 7 shows the deviations of the CEs from the true values for the nine combinations of material maps (iodine, tungsten, and water) and inserts (iodine, tungsten, and water).We found average false estimates of 0.4 mg/mL in the iodine and in the tungsten map and 11.9 mg/mL in the water map.This means an average error of 4% in the CA concentration with respect to the insert concentration of 10 mg/mL and 1% in the water map with respect to the total concentration of 1000 mg/mL.Most of the false estimates in the CA inserts except for the tungsten map-iodine insert case were due to systematic effects, viewing propagated quantum noise of 0.01 mg/mL, 0.03 mg/mL, and 0.9 mg/mL on each data point in the iodine map, the tungsten map, and the water map, respectively.The CE in the water insert, with maximum errors below 2% with respect to the total concentration, was more accurate than in the CA inserts.In the iodine insert, the iodine concentration was overestimated with medians of 5 and 8% for low and high phantom off-centering, respectively.Meanwhile, the water concentrations were underestimated with a median of at most 1% independent of the phantom position.In the tungsten insert, the tungsten concentration was overestimated with medians of 6% and 7% for low and high phantom off-centering, respectively.Meanwhile, the iodine and water concentrations were underestimated with a median of 2% or less independent of the phantom position.The entirety of false estimates in the CA map-CA insert pairs tended to be overestimated.
For each map-insert pair, the 15 mm off-centered phantoms showed more accurate results than the 45 mm off-centered phantoms.The central quartiles did not show systematically higher false estimations.
We would like a conservative upper boundary for systematic CA CE errors at arbitrary CA concentration .As a voxel's absorption equals the sum of the absorptions of all included materials in the thin-absorber limit and as the CE is linear in the voxel's CT number quadruple, the largest errors measured at  = 0 (pure water) and  = 10 mg/mL, denoted as  and  , implied the following upper boundary: The CE in the water insert, with maximum errors below 2% with respect to the total concentration, was more accurate than in the CA inserts.In the iodine insert, the iodine concentration was overestimated with medians of 5 and 8% for low and high phantom off-centering, respectively.Meanwhile, the water concentrations were underestimated with a median of at most 1% independent of the phantom position.In the tungsten insert, the tungsten concentration was overestimated with medians of 6% and 7% for low and high phantom off-centering, respectively.Meanwhile, the iodine and water concentrations were underestimated with a median of 2% or less independent of the phantom position.The entirety of false estimates in the CA map-CA insert pairs tended to be overestimated.
For each map-insert pair, the 15 mm off-centered phantoms showed more accurate results than the 45 mm off-centered phantoms.The central quartiles did not show systematically higher false estimations.
We would like a conservative upper boundary for systematic CA CE errors at arbitrary CA concentration c.As a voxel's absorption equals the sum of the absorptions of all included materials in the thin-absorber limit and as the CE is linear in the voxel's CT number quadruple, the largest errors measured at c = 0 (pure water) and c = 10 mg/mL, denoted as ε 0 and ε 1 , implied the following upper boundary: For the phantom setup used during the test data acquisition, we expected the iodine and tungsten CE in a voxel of CA concentration c for a well-positioned object ( x 5 = 15 mm) to suffer from a systematic inaccuracy of less than 0.2 mg/mL + 0.11•c.

Contrast-to-Noise Ratio
For each CA, CNR out /CNR in varied by only a few percent across the different configurations.On average, the CNR declined to 36% and 16% for iodine and tungsten maps, respectively.
The stronger CNR reduction for tungsten and thus higher noise in the tungsten map compared to iodine, visible in Figure 5, is not surprising: Iodine's absorption curve, on the one hand, strictly decreases in the diagnostic energy range and leads to monotonously decreasing CT numbers with increasing threshold energy, which separates well from water whose CT number is assumed to be identical in all thresholds.Tungsten's absorption curve, on the other hand, contains a K-edge, which leads to a relatively weak variation of the threshold CT numbers and consequently to a worse separation from water.
We attribute the trend to a weaker CNR reduction for larger objects to stronger beam hardening and thus better spectral resolution of the system.However, for larger objects irradiated with equal X-ray tube flux, the noise in spectral input images is higher due to a lower photon flux on the detector.So, despite a slightly weaker CNR reduction, the material map is noisier for larger objects as can be seen in Figure 5.

Estimation Accuracy
Figure 7 shows mainly systematic CE errors resulting from applying the created models for a MMD of the test data into iodine, tungsten, and water.Lower false estimates in the case of low phantom off-centering support the best practice guideline to adjust patients as close to the isocenter as possible.As most absolute errors are on the same level as in the low off-centering case for all map-insert pairs, our approach shows certain robustness against object positioning.
We explain the observed tendency to overestimate CA concentrations in the CA inserts in the corresponding inserts by the ignored parameter x 4 in the CA models.We used x 4 = 30 mm throughout the calibration scans but used larger values of x 4 for the test data acquisition.Beam hardening induces a higher CT number for CA closer to the object boundary resulting in higher CEs.
We attribute the water concentration underestimation in the CA inserts and the underestimation in the tungsten map-iodine insert case to the surrounding materials equivalent to soft tissue, liver, spleen, and spinal column.Their CT number decreases with increasing threshold energy.Given x 1 , based on the CT numbers in T1, we thus expect a higher photon flux through the CA insert in the upper thresholds T2, T3, and T4 in the abdominal phantom than in a phantom of water-equivalent material with equal x 1 .This surplus is larger with larger photon energy.
The higher flux passing through the CA inserts decreases the CT numbers in the upper thresholds.This is compensated by an increased CA CE and a decreased water CE both for the iodine and the tungsten insert.The tungsten insert requires considering the K-edge which leads to a higher CT number in T3 than in T2 and a lower CT number in T4 than in T2.The increasing photon surplus with increasing threshold energy enlarges the CT number difference between T3 and T4 more than between T2 and T3.A negative iodine concentration compensates for this relatively lower difference between T2 and T3.
Thus, for targets whose absorption is not dominated by water, the WED is no longer a sufficient parameter to describe the pre-attenuation for an image voxel of interest.
A further increase in accuracy could be achieved by adding more parameters to the models.However, this would increase the number of required calibration scans, thus increasing the time needed to calibrate each material and hence limiting the method's practicability.
We assume the achieved estimation accuracy and CNR to be sufficient for applications aiming to differentiate image regions with different content, such as coronary CT angiography (CCTA) or stent imaging which have already been tested with tungsten-based CA [20,22].According to the results in Section 3.2, iodine and tungsten with a concentration of 3 mg/mL or above are detectable without a problem, despite systematic errors resulting from MMD into iodine, tungsten, and water.According to the results in Section 3.1, we predict the exemplary contribution of quantum noise to the CA CE in a volume of 16 mL, like a medium plaque [48], scanned in a CCTA with exposure 255 mAs [49].In a patient with a medium WED of 300 mm, the standard deviation on the CE in the averaged volume equals 0.3 mg/mL and 0.7 mg/mL for iodine and tungsten.This makes the medium plaque easy to distinguish from surrounding iodine-or tungsten-based CA.
Finally, Equation ( 12) may not be valid for arbitrarily high insert concentrations in the examined target due to a broken thin-absorber limit and thus non-proportional CT number and insert concentration, particularly if photon starvation reduces the incident signal below the detector sensitivity.

Conclusions
We presented a method for performing a MMD based on multiple energy-threshold input images from an energy-discriminating PCD CT system.The calibration does not require manually fine-tuned initial conditions.No prior knowledge is required to apply the calibrated models.All model computations are analytical, resulting in a fast and robust algorithm execution.It was designed for anthropomorphic geometries at arbitrary table heights and exposures.
We applied our approach to a commercial PCD scanner using a set of three BMs with suitable properties: Iodine, tungsten, and water.We evaluated its accuracy for several phantom configurations and scan parameters representative of typical clinical protocols and examined body parts.The investigation of more complex setups as a study involving biological specimens is planned as follow-up of this work.
We used a non-iterative image reconstruction to be as close to the statistical reality as possible.In clinical practice, an iterative reconstruction with lower quantum noise per image voxel and thus lower contribution of quantum noise to CE on each material map voxel is usually used [50].
The measured objects were axially homogeneous, whereas some clinical applications with axial inhomogeneity require a smaller slice thickness than what was used in this work.We scanned 20 mm diameter syringes with CA and water and evaluated volumes of 1.6 mL in the syringes to assess the algorithm-inherent performance.We expect the results for the algorithm's precision and accuracy to be valid also for smaller objects.Investigating very small objects such as narrow vessels with sizes close to the CT system's spatial resolution will introduce additional errors due to the partial volume effect.
As the images acquired to test the method was used to optimize the model types, both datasets are not entirely independent.The method should thus be tested with additional measurements in different setups for a more reliable validation.
The chosen threshold energy quadruple provides the best spectral separation between tungsten and the other BMs.Using another predefined threshold energy quadruple, we expect a stronger CNR decline in the material images and larger CE errors.Further studies should also investigate the reliability of the presented approach when introducing a BM with less optimal K-edge regarding the spectral bins or thresholds.
The method is not limited to a three-material decomposition.As the used system provides four spectral datasets per scan, a four-material decomposition is possible.Upcoming PCD CT scanners will likely be able to use more than four spectral datasets per scan.Decomposition into as many materials as there are spectral datasets is then possible.The quality of such MMDs would be highly dependent on the specific choice of BMs and scan parameters, especially on the optimal choice of spectral bins or thresholds.

Figure 1 .
Figure1.Summary of the general approach of obtaining  material maps from  spectral input images.Model creation: the covariance of the spectral images is measured to model the image noise covariance matrix.Averaged CT numbers of contrast agent-filled image regions are corrected by those of water before modeling.A model of the water CT numbers is built separately.Decomposition: given  input images, each individual voxel is corrected using the water model and then fed into the decomposition algorithm, whose output translates into material concentrations.

Figure 1 .
Figure 1.Summary of the general approach of obtaining M material maps from N spectral input images.Model creation: the covariance of the spectral images is measured to model the image noise

'
Best' means that among all linear functions in of Equation (2) has the lowest variance.Note that in the case of M = N, Equation (2) is identical to Equation (3) and can be solved without LS methods and the INCM.Furthermore, in the case of M < N, the solution of Equation (2) has a lower variance than the ordinary LS estimator S T S −1 S T → y

Figure 2 .
Figure 2. Object geometry.(a) Definition of the distance measures.(b) Exemplary evaluated circular region defined by its center's distance  4 to the phantom center and the angle  to the horizontal line through the phantom center.(c) Calibration phantom for the exemplary case of  = 200 mm.We evaluated a circular region at  4 = 30 mm and β = 180°, which is in the drill hole filled by one of the inserts.Moreover, at several angles β (−90° −45° 0° +45° +90°) and distances  4 to phantom center, a circular region in WEP was evaluated.For each scanned phantom diameter , we used  4 = 30 mm and at  ≥ 200 mm and  ≥ 300 mm , we additionally included  4 = 70 mm and 120 mm, respectively.

Figure 3 .
Figure 3. CT number of iodine solution (c = 10 mg/mL) positioned in a water-equivalent plastic cylinder of diameter  in the image using all photons above 90 keV of an energy-discriminating photon-counting detector CT system at 140 kV tube voltage and after water correction, described in Section 2.1.2,with fixed parameters as shown. denotes the distance of the solution to the isocenter.The error bars indicate the standard error on the iodine solution region mean (averaged over 1.6 mL) of the individual points due to quantum noise.To calibrate the models to different parameters, the user performs calibration scans with different parameter configurations  = { 1 …   } with  equal to or larger than the number of coefficients  of the model function ( ;  ), for each spectral image and CA in the BMs.This leads to a set of ground truth values = {(  ;  ) ∶  = 1 … } of

Figure 2 .
Figure 2. Object geometry.(a) Definition of the distance measures.(b) Exemplary evaluated circular region defined by its center's distance x 4 to the phantom center and the angle β to the horizontal line through the phantom center.(c) Calibration phantom for the exemplary case of D = 200 mm.We evaluated a circular region at x 4 = 30 mm and β = 180 • , which is in the drill hole filled by one of the inserts.Moreover, at several angles β (−90 • , −45 • , 0 • , +45 • , +90 • ) and distances x 4 to phantom center, a circular region in WEP was evaluated.For each scanned phantom diameter D, we used x 4 = 30 mm and at D ≥ 200 mm and D ≥ 300 mm, we additionally included x 4 = 70 mm and 120 mm, respectively.

Figure 2 .
Figure 2. Object geometry.(a) Definition of the distance measures.(b) Exemplary evaluated circular region defined by its center's distance  to the phantom center and the angle  to the horizontal line through the phantom center.(c) Calibration phantom for the exemplary case of  = 200 mm.We evaluated a circular region at  = 30 mm and β = 180°, which is in the drill hole filled by one of the inserts.Moreover, at several angles β (−90°, −45°, 0°, +45°, +90°) and distances  to phantom center, a circular region in WEP was evaluated.For each scanned phantom diameter , we used  = 30 mm and at  200 mm and  300 mm , we additionally included  = 70 mm and 120 mm, respectively.

Figure 3 .
Figure 3. CT number of iodine solution (c = 10 mg/mL) positioned in a water-equivalent plastic cylinder of diameter  in the image using all photons above 90 keV of an energy-discriminating photon-counting detector CT system at 140 kV tube voltage and after water correction, described in Section 2.1.2,with fixed parameters as shown. denotes the distance of the solution to the isocenter.The error bars indicate the standard error on the iodine solution region mean (averaged over 1.6 mL) of the individual points due to quantum noise.

Figure 3 . 1 ,
Figure 3. CT number of iodine solution (c = 10 mg/mL) positioned in a water-equivalent plastic cylinder of diameter D in the image using all photons above 90 keV of an energy-discriminating photon-counting detector CT system at 140 kV tube voltage and after water correction, described in Section 2.1.2,with fixed parameters as shown.L denotes the distance of the solution to the isocenter.The error bars indicate the standard error on the iodine solution region mean (averaged over 1.6 mL) of the individual points due to quantum noise.To calibrate the models to different parameters, the user performs calibration scans with different parameter configurations C = { → x 1 , . . ., → x K } with K equal to or larger than

Figure 4 .
Figure 4. Test scan setup with obese extension ring.In the abdominal phantom holes, each with one insert: One in the liver-equivalent region, one in the spleen and one in the soft tissue-equivalent region.The inserts were cyclically swapped mentioned in Table1, to measure each one in all the different regions.The cent 100 mm is filled by water-equivalent plastics.

Figure 4 .
Figure 4. Test scan setup with obese extension ring.In the abdominal phantom, we fill three drill holes, each with one insert: One in the liver-equivalent region, one in the spleen-equivalent region, and one in the soft tissue-equivalent region.The inserts were cyclically swapped between the scans, mentioned in Table1, to measure each one in all the different regions.The central hole of diameter 100 mm is filled by water-equivalent plastics.

Figure 5 .
Figure 5. Material decomposition input (T1 image) and output (material maps with contrast agent inserts visible as bright patch in the respective map) of the test data configuration with tube current 750 mA, 15 mm off-centering, tungsten in the left and iodine in the right insert (a) without obese extension ring and (b) with the obese extension ring.Peripheral image regions were cropped.

Figure 5 .
Figure 5. Material decomposition input (T1 image) and output (material maps with contrast agent inserts visible as bright patch in the respective map) of the test data configuration with tube current 750 mA, 15 mm off-centering, tungsten in the left and iodine in the right insert (a) without obese extension ring and (b) with the obese extension ring.Peripheral image regions were cropped.

Figure 5 .
Figure 5. Material decomposition input (T1 image) and output (material maps with contrast agent inserts visible as bright patch in the respective map) of the test data configuration with tube current 750 mA, 15 mm off-centering, tungsten in the left and iodine in the right insert (a) without obese extension ring and (b) with the obese extension ring.Peripheral image regions were cropped.

Diagnostics 2024, 14 , 1262 12 of 18 Figure 7 .
Figure 7. Errors in the estimated material concentrations of the test data using the models described in Section 2.1.There is one column for each of the inserts.The box construction is analogous to Figure 6.Positive errors indicate an overestimated concentration, while negative errors indicate an underestimation.Δ is the root mean square of all false estimations in units of mg/mL.

Figure 7 .
Figure 7. Errors in the estimated material concentrations of the test data using the models described in Section 2.1.There is one column for each of the inserts.The box construction is analogous to Figure 6.Positive errors indicate an overestimated concentration, while negative errors indicate an underestimation.∆ is the root mean square of all false estimations in units of mg/mL.
For each spectral image, a model of the CT numbers is created for the BMs and water even if water is not among the BMs as discussed in Section 2.1.2.Calibration measurements are performed to find ground truth at different parameters Having acquired the ground truth, the user calibrates H based on G H , g for each spectral image based on G g , and f for each remaining BM and each spectral image based on G f and G g , as discussed in Sections 2.1.2and 2.1.3.
the model of a CA CT number in one spectral image with coefficients → α to calibrate and a set of ground truth data G f ; • g = g → x ; → α is the model of the water CT number in one spectral image with coefficients → α to calibrate and a set of ground truth data G g ; • H = H → x ; → α is the model of the INCM with coefficients → α to calibrate and a set of ground truth data G H . → x for each model.

Table 1 .
CT scan setup for the acquisition of the test data with the phantom vis one insert configuration.In total, 360 scans were performed.

Table 1 .
CT scan setup for the acquisition of the test data with the phantom visible in Figure4for one insert configuration.In total, 360 scans were performed.

Table 2 .
Fixed parameters for the scans.

Table 3 .
CT scan setup for the calibration measurements of each base material with phantom configuration shown in Figure2.For each base material, 75 scans were performed in total.