Material decomposition maps based calibration of dual energy CT scanners for proton therapy planning: a phantom study

We introduce a new calibration method for dual energy CT (DECT) based on material decomposition (MD) maps, specifically iodine and water MD maps. The aim of this method is to provide the first DECT calibration based on MD maps. The experiments were carried out using a general electric (GE) revolution CT scanner with ultra-fast kV switching and used a density phantom by GAMMEX for calibration and evaluation. The calibration process involves several steps. First, we tested the ability of MD values to reproduce Hounsfield unit (HU) values of single energy CT (SECT) acquisitions and it was found that the errors were below 1%, validating their use for HU reproduction. Next, the different definitions of computed Z values were compared and the robustness of the approach based on the materials’ composition was confirmed. Finally, the calibration method was compared with a previous method by Bourque et al, providing a similar level of accuracy and superior performance in terms of precision. Overall, this novel DECT calibration method offers improved accuracy and reliability in determining tissue-specific physical properties. The resulting maps can be valuable for proton therapy treatments, where precise dose calculations and accurate tissue differentiation are crucial for optimal treatment planning and delivery.


Introduction
Proton therapy treatments have widely grown in recent years, mainly due to the clinical benefits and availability of this therapy.The clinical benefits of proton therapy are related to the radiobiology of dose deposition using proton beams and the reduction of the overall dose and side effects throughout treatment (Paganetti et al 2019).Traditionally, proton therapy has been recommended for a wide variety of cancers, including the pediatric, ocular, base of the skull, and re-irradiated tumors, among others (García 2019).Also, new ongoing studies and clinical trials are investigating clinical benefits for new tumor locations such as head and neck, and breast (Tambas et al 2020, Mutter et al 2021), among others.Regarding availability, the number of dedicated facilities and approved projects keeps growing worldwide.Particularly in Europe, there are more than ten new centers projected for the next years (Dosanjh et al 2018).
Proton therapy treatments are characterized by a high level of technology and precision in all the steps of the radiotherapy procedure, i.e. from beam generation and delivery to treatment planning and dose calculations, including patient immobilization.Additionally, daily patient positioning should be properly monitored via image-guided radiotherapy (IGRT), and adaptative radiotherapy (ART) strategies should be accounted for to ensure the correct delivery of planned dose distributions.The reason for these precautions is the uncertainties that affect proton treatments.Therefore, due to the finite range of the proton beams in water and human tissues, the reproducibility of the patient positioning and other aspects such as the filling of cavities and the constancy of the tissue characteristics in the proton beam path should be verified (Jones et al 2007).
In addition to the aforementioned IGRT and ART strategies for treatment delivery control, other useful tools such as the robust planning technique together with Monte-Carlo dose calculations are commonly used in proton treatments to minimize the impact of uncertainties associated with the dose calculation process itself (Paganetti 2012, Moyers et al 2020).In this case, there are two main root causes that explain this uncertainty: the stopping power ratio (SPR) characterization of the biological tissues and the calibration of the CT units employed to derive SPR in every voxel of the virtual simulation image (Schneider et al 1996, Ainsley andYeager 2014).It should be noted the dependency of the biological tissues' SPR in both the relative electronic density to water ( e W r ) and the chemical composition given by the effective atomic number (Z ef ) that allows obtaining the ionization potential of the tissue (I) (Bär et al 2018).
Achieving the greatest possible differentiation between tissues is a key point to reduce uncertainty in the calculation of the absorbed dose during treatment planning, but there exist many limitations that challenge the differentiation of materials, leading to confusion between tissues.In particular, for a given acquisition energy, CT numbers depend on both magnitudes, Z ef and e W r , so tissue differentiation is limited by the lack of information provided by single-energy CT (SECT) acquisitions.Performing CT acquisitions at multiple energies could decouple the influence of density and chemical composition, providing additional information (Johnson et al 2007).Indeed, Hounsfield already mentioned the potential of using two images to enhance tissues that cannot be differentiated using just one image (Hounsfield 1973).However, its implementation was not possible at that time given the lack of hardware development and the already high dose exposition to the patient during a single acquisition.
Dual energy CT (DECT) technology solves the uncertainties that the first CT generations introduced.Although the underlying principles of DECT are the same regardless of scanner type, single-source dual-energy CT (SS-DECT) scanners with fast kilovoltage switching (KVS-DECT) differ from dual-source dual-energy CT (DS-DECT) and dual-layered CT (DL-DECT) scanners both in data acquisition and in data processing methods.While KVS-DECT outputs are material decomposition (MD) maps and virtual monochromatic images (VMI), the other scanners offer High and Low Energy acquisitions (HU High and HU Low, respectively) together with VMI images.
Despite the major improvements of recent DECT scanners, inter and intra-scanner differences in output, that are transferred to measured HUs, still occur due to factors such as non-deterministic electronics and variability in tube output spectra, among others.Thus, calibration methods are required to account for variabilities and determine parameters to assess tissue-specific physical properties.Over the years, different calibration methods have been introduced for SECT and DECT calibrations.For instance, in the case of SECT scanners for proton therapy centers, stoichiometric calibration (Schneider et al 1996) is the current state of the art.
DECT calibrations demonstrated better performance, improving the e W r , Z ef , and SPR computations in comparison with the SECT approaches (Taasti et al 2018, McCollough et al 2020).These calibrations mostly are optimization algorithms that convert HU High and HU Low to e W r /Z ef (Saito 2012, Bourque et al 2014).However, the aforementioned articles focus on the calibration of DS-DECT scanners and cannot be directly applied to KVS-DECT scanners, which provide no HU pairs, as HU pairs are not typically available from these scanning systems.
To the best of the authors' knowledge, there is only one previous scientific reference aimed at calibrating KVS-DECT (Näsmark and Andersson 2021).This unique KVS-DECT method derives SPR maps by directly applying the formalism of Jackson and Hawkes (Hawkes and Jackson 1980) to the optimal VMI pair in the reconstruction interval.However, VMIs do not present one-to-one correspondence with the SECT HUs since they are a linear combination extracted from the MD maps.Additionally, this method is presented as an initial proof of concept that needs further fine-tuning for some tissue substitutes.
As MD maps offer information with physical meaning, we hypothesize that these maps can be employed for direct DECT calibration.Thus, in this work, we present and evaluate a novel calibration method based on the MD maps of iodine and water to obtain e W r , Z ef , I, and SPR maps.

Materials
Dual and single energy acquisitions have been performed using a general electric (GE) Revolution CT scanner (General Electric Healthcare, Waukesha, WI, USA).This scanner includes the Gemstone Spectral Imaging system by GE, which employs ultra-fast kV switching during DECT acquisition, producing nearly simultaneous projections using a unique pair source-detector.
The acquisition pipeline consists of one KVS-DECT scan from which iodine and water MD maps are reconstructed.In addition, two SECT acquisitions are performed at 80 and 140 kVp for evaluation purposes, although they are not necessary for the proposed calibration method.All image volumes have been acquired with a dynamic range of 16 bits, applying a metal artifact reduction algorithm and adaptative statistical iterative reconstruction (ASIR) with a parameter value of 60.The rest of the acquisition parameters are shown in table 1.
The density phantom by GAMMEX (Sun Nuclear, Florida, USA) has been used for calibration and evaluation purposes.The phantom chassis is composed of two separable parts with 16 available places for inserts, as shown in figure 1(a).Table 2 summarizes all the information about the specific inserts used in this work.The phantom's user certification has been supplied by the manufacturer in a private communication and provides the composition, mass density, and ρ e of each tissue surrogate.
The proton therapy dosimetry head phantom by CIRS (Sun Nuclear, Florida, USA) has been used for evaluation purposes.It approximates the average male human head in both size and structure, as shown in figure 1(b).The phantom is constructed of tissue-equivalent materials, which mimic reference tissues within 1.5% for protons, according to manufacturer specifications.It should be also mentioned that the composition of the different materials has been kindly provided by the manufacturer in a private communication 2.2.Basic definitions of the physical magnitudes.The HU images provided by SECT studies are related to the linear attenuation coefficient of the material, μ, in every voxel of the map.Given that for a particular energy E, the linear attenuation coefficient of the water μ W can be known, the value of μ for every voxel can be obtained through the following relationship: According to the work of Möhler et al (2018), the attenuation linear coefficient is related to Z ef and ρ e , through equation (2): where coefficients S and T are dependent on the particular E used for the image acquisition.Also, Z ef and ρ e are defined based on the composition of the material.In particular, for a material with a density ρ, which is composed of N elements with atomic numbers Z i , mass numbers A i , and percent weight composition w i , by using the Avogrado constant, N A , it can be stated that: In this work, we have decided to use equations (3) and (4) in order to define ρ e and Z ef based exclusively on the materials composition, as is the case of other works (Xie et al 2018).It is worth noting that there are other definitions for Z ef that take into account the particular energy spectra of the imaging system (Bonnin et al 2014, Bourque et al 2014).In this case, the particular x-ray spectrum should be accounted for in order to properly compute Z ef .
The SPR of a particular material, relative to water, can be calculated with the Bethe equation using the rest of electron mass, m e ; the electron density of the material, ρ e ; the electron density of the water, ρ e,W ; the speed of light in vacuum, c; the relative speed of the proton to light, β; and the mean ionization values of the material and water, I and I W , according to equation ( 5 Finally, it should be noted that the mean excitation energy of the material, I, is linked to its chemical composition, i.e. its Z ef .Among all options proposed in the literature, we use the relationship derived by Bourque et al (2014): In equation (6), the coefficients take the following values in the original work: a1 = 14.007762, a2 = -21.414214,a3 = -0.005342,a4 = 0.207079, a5 = -2.589844,a6 = 8.339473, a7 = 51.895887,a8 = -219.722173,a9 = 11.794847, and a10 = -47.707141.

Design of the study
This study is aimed at providing a novel calibration method for DECT studies based on MD maps and evaluating its performance versus a known DECT calibration method based on low and high HU images.SECT studies are acquired only to simulate the HU low and high acquisitions (80 and 140 KVp) and to allow applying the method proposed in Bourque et al (2014).The main hypothesis here is that SECT images of different energies may be employed as representative HU low and high images of a DECT study, but this methodology has been previously adopted in the work of Garcia et al (2015) to simulate DECT acquisitions.
In order to accomplish the aim of the study, some intermediate steps need to be previously completed and their potential impact should be properly evaluated to ensure the consistency of the MD calibration for DECT studies.Therefore, results from the intermediate steps will be presented and discussed to properly address all the inputs and procedures needed for the calibration.
First, we test the ability of the MD values to reproduce the HU values of SECT acquisitions.According to previous results from Lehmann et al (1981), the choice of the basis function employed as the root of the MD leads to errors lower than 1%, as detailed in section 2.4.However, according to McCollough et al (2020), the errors of the MD values when reproducing mass attenuation coefficient values may be as high as 3.5% depending on the energy spectrum.Although these errors may be partially overcome in the calibration process, it is important to ensure that the MD values provided by the imaging system are in good agreement with the HU values for different energy spectra frequently used in computed tomography such as 80, 120 and 140 KVp.
Second, we test the validity of the computed Z ef values of the inserts used in this work.As mentioned above, several definitions of Z ef can be found in the literature (Liu et al 2021) that may give rise to differences between the estimates of up to 10% depending on the material.Although providing a definition for Z ef values is beyond the scope of this work, we want to be consistent with published protocols and to show that this point introduces no bias in this study.Thus we compare in section 2.5 the definition used in this study based on the known composition of the materials versus the one introduced in the work of Bourque et al (2014), which requires knowledge of the x-ray spectrum and differential Z ef values.To this purpose, we employ the spectrum from the work of Terini et al (2020) for a GE DECT scanner, kindly provided by the authors in private communication, and differential μ values obtained from the XCOM database (National Institute of Standards and Technology, Gaithersburg, Maryland, USA).
Third, we evaluate the performance of the new DECT calibration based on water-iodine MD maps and compare it to the method of Bourque et al 2.6.In this evaluation, we present the fits of all the steps of the new method and discuss the numerical errors produced by both calibrations, versus the known values.As previously pointed out in the bibliography (Hünemohr et al 2013), the lung tissue surrogates may not work properly for calibration procedures and they can be removed from the fitting algorithm although they may be included in the error evaluation.It is beyond the scope of our work to give the final solution to this problem.However, we evaluate different calibration options and propose a way to account for this extra uncertainty source in the calibration protocol for both DECT calibration protocols.
Finally, the ρ e , Z ef , I, and SPR maps are presented for both methods and compared.Also, the maps are qualitatively and quantitatively analyzed to show differences in the final outcomes of the calibration procedure.

The material decomposition (MD) theory
The mass attenuation coefficient of a given material for an x-ray beam of a particular energy E may be written as a sum of three different terms related to the main photon interactions: photoelectric (P), compton (C), and coherent (Co) (McCollough et al 2020) As shown in the bibliography (Alvarez andMacovski 1976, Lehmann et al 1981), when k-edge effects are negligible and the effective energy is below 30 KeV, the coherent interaction effect may be discarded, and with a particular basis functions choice, equation (7) may be approximated with an error lower than 1% by: Regarding equation (8), it is worth mentioning that the photon energy dependency is fully contained in the well-known energy functions f P (E) and f C (E), while α P and α C are only related to the material composition.In particular, if we account for the atomic mass of the material, A, and the atomic number, Z, we obtain the relationship shown in equations ( 9) and (10), where K P and K C are the constants associated to the compton and photoelectric effect equations Now, for a given energy E, we can choose two materials as a basis, M and N (water and iodine in our case), and we can rewrite equation (8) for these materials as follows By working with equations (11) and (12), f P (E) and f C (E) may be expressed in terms of the particular characteristics of the materials M and N, i.e.
It should be noted that these characteristics of the materials M and N are constants for the given energy E, so equation (8) can be rewritten for every material with characteristic values, α P and α C , using the two material bases as: Finally, the linear coefficient of a material for a particular energy, μ(E), may be written in terms of the concentrations of the bases materials, ρ M and ρ N , as follows: It should be noted that ρ M and ρ N have dimensions of concentration, but they are not real concentrations.Indeed, they can even present negative values as they depend on the constants associated with the materials of the basis, M and N, and on the particular material in the voxel with ρ, Z, and A. In fact, by combining equations ( 9) and (10), ( 13), ( 14) it can be shown that: It should be noted that constants C 1 to C 5 have been introduced for the sake of clarity.These constants depend on the basis material characteristics α M,P , α M,C , α N,P , α N,C , and K C and K P .

The DECT calibration method based on MD maps
The starting point of our calibration are equations (15) and (16).By accounting for the definition of the electron density of a material, shown in equation (4), the output of the MD maps, ρ M and ρ N , can be written as a function of the characteristics of the material in each voxel, ρ e and Z ef , as: From equations (17) and (18), we can define new quantities related to the outputs of the MD maps, ρ M and ρ N , that allow us to obtain ρ e and Z ef .In particular, we can define X that only depends on the Z ef as: Also, we can define the quantity Y that can be written as a function of ρ e and Z ef as: Subsequently, the input calibration data, ρ M and ρ N , are obtained from the MD maps by computing mean and standard deviation in the ROIs shown in 1 for every insert.Also, the insert compositions provided by the GAMMEX phantom manufacturer are considered as the other input data needed for the calibration.Finally, with the whole input data, fits can be carried out and f X and f Y , which are the outputs of the calibration procedure, can be properly determined.
Once these functions are known from the calibration data, they can be used to estimate ρ e (or e W r , as ρ e,W may be considered as a constant), and Z ef maps from a particular study.At this point, it should be mentioned that according to equations (15) and ( 16), f X and f Y may be well represented by monotonic functions, as polynomial functions of low order, or a sum of exponential functions.In our case, a three-degree polynomial was found to properly fit the experimental data for both functions.

The Bourque et al method: brief summary
The Bourque et al method for DECT calibration (Bourque et al 2014) works with two HU images acquired with two different energy spectra, low and high, and begins with equations (1) and (2) to define the reduced HUs, u L and u H , and their dependence on Z ef and ρ e as follows: Based on equations ( 21) and ( 22), the dual energy ratio (DER) may be defined as the relationship between u L and u H , and it is clear that this magnitude only depends on the Z ef .
Then, the experimental data u L and u H can be employed to properly determine f DER , f L , and f H through fittings to polynomial functions.According to the original work of Bourque et al (2014), the degree of the polynomial functions may vary depending on the experimental data.In our case, three-degree polynomials were found to meet the accuracy requirements.Finally, the knowledge of these functions is the output of the calibration procedure and makes it possible to estimate Z ef and ρ e (or e W r ) on a voxel basis from a DECT study.

Uncertainty analysis
In order to estimate the uncertainty of the calibration procedures, we consider the standard deviation in the ROIs of the ρ e , Z ef , I, and SPR maps for every insert.In this way, we can evaluate the real output of the calibration procedure without propagating uncertainties and can estimate the behavior of the calibration procedures in a clinical scenario to compare the performance of both methods in terms of uncertainty.

The MD theory
Figure 2 displays in section (a) a front view of the GAMMEX phantom showing the actual configuration of the tissue substitute inserts used during our experiments, and a cross-section of (b) the iodine MD map depicting the segmented ROIs, and (c) the water MD map.We estimate the values for each insert as the mean and standard deviation measured in the selected ROIs in each of the MD maps and each of the three HU datasets acquired at different voltages, 80, 120, and 140 KVp.
Using these experimental data, a 2D linear fit was performed for every SECT energy using equation (14).In this way, the ROI values coming from the MD maps provided the ρ M and ρN input data for the basis materials, i.e. water and iodine, while the ROI values coming from the SECT images, provided the μ(E) data for every SECT value.In this way, ( ) ( ) are obtained for 80, 120, and 140 KVp, as the result of the fitting process.The residual errors of these fits are plotted in figure 3 for every SECT energy considered.Finally, as described in the bibliography (McCollough et al 2020), equation ( 14) should be valid for the lower and higher KVp values employed in the DECT acquisition, i.e. 80 and 140 KVp in our case.Thus, by evaluating the results in figure 3, we see that the residual error of the MD maps of our machine is lower than 1.2% when using energies higher than 80 KVp.

Z ef estimates
As stated in section 2.2, we decided to use equation (3) to define the values of Z ef based exclusively on the materials composition extracted from the GAMMEX user's certificate.On the contrary, Bourque selected a different definition for Z ef that accounts for the energy spectra of the imaging system.
In figure 4, the differences between our Z ef estimates and Bourque ones are shown for the lung, liver, bone, adipose tissue, breast, brain, and calcification inserts.As may be seen, the differences between the two estimates may differ up to 0.8, which is roughly a 10% of a percent difference.

New DECT calibration procedure
Following the same procedure stated before in 3.1, mean and standard deviation for each tissue substitute insert are extracted from the MD maps by segmenting ROIs (shown in figure 2(b)).Subsequently, this data is fit into equations ( 19) and (20) to obtain SPR maps.The coefficients for the exponential and the polynomial fits are presented in table 3, where the x and y data have been normalized using mean and standard deviation values of 0.9395  0.1052 and 8.743  2.225 respectively.

DECT calibration outputs
The comparison between the two DECT calibration methods (Bourque et al and the proposed method, hereafter referred to as URJC-QS) is shown in figure 6, where the errors produced by both methods are shown for calibrations without excluding any insert.It can be observed how both methods give rise to very similar residual errors.Also, the absolute maximum deviations are found for the lung inserts in both cases, Z ef and e W r , and a deviation close to 5% is found for the blood insert with both methods in the e W r evaluation.

Use of lung inserts
In table 4, the absolute errors through the whole calibration set of inserts are studied when the lung inserts are included or excluded from the calibration process but considered for the evaluation of the method's goodness.As shown, the method of Bourque is pretty insensitive to this consideration.However, for our method, the maximum absolute errors are incremented especially for the Z ef estimates.In particular, these maximum deviations are always related to these lung inserts.Furthermore, the mean absolute error is quite independent of the inclusion or exclusion of the lung inserts for both methods.Finally, it may be regarded that the mean absolute error of the e W r estimates of both methods are very similar, while the mean absolute errors of the Z estimates are about 1% higher for our method.
3.6.Z ef , e W r , I and SPR maps Two different analyzes have been conducted on the Z ef , e W r , I and SPR maps.First, the GAMMEX phantom maps have been created and a ROI analysis has been conducted for every map by following the ROI definition shown in figure 2. So, the output of the two calibration methods for the magnitudes of interest is shown in figure 7. Unlike the results in figure 6, where only the calibration data is considered for the residual errors computation, the results in figure 7 are obtained from the maps.In both cases results are compared to the reference values computed using phantom certificate compositions shown in table 2. The error bars account for the uncertainty computed from the maps as the ROI standard deviations for every insert.As may be seen, the maximum deviations and the maximum uncertainties are found for the lung inserts, as was the case with the residual errors found in figure 6.It should be noted that the blood insert also produces bad accuracy results for both calibrations with deviations close to 5% in e W r and SPR evaluations.It also should be pointed out that the Bourque method gives rise to higher experimental uncertainties, especially for lung inserts.In particular, when accounting for all the inserts in the study, the mean uncertainties of all the inserts were 9.0%, 2.4%, 8.1% and 3.5% in the case of the Bourque et al protocol for the Z ef , e W r , I and SPR maps; while in the case of our protocol, the mean uncertainties were 2.6%, 1.3, 1.8%, and 1.1% respectively.
In second place, the Proton Therapy Dosimetry Head phantom is also employed to compute the maps.In particular, the Z ef , e W r , I, and SPR maps and their relative differences are shown in figure 8.In accordance with the ROI analysis, the Bourque method gives rise to maps with lower definition and higher noise levels.This amount of noise is translated to the percent difference maps, although the absolute mean percent differences are Table 4. Absolute errors (in %) when the lung inserts are considered as part of the calibration procedure or excluded from it.These absolute errors are computed as the mean and the maximum deviations obtained when evaluating the calibration in all the inserts in table 2 regardless of the inserts excluded for the calibration process.

No exclusion
Lungs lower than 3%.In addition to the map evaluation, four detailed ROIs matching specific tissues have been introduced as shown in figure 9.The results of the ROI evaluation are presented in table 5 where it can be seen how the Bourque method gives rise to higher uncertainty levels, while the agreement between both methods is better than 3%.

Discussion
The MD theory is able to reproduce HU values for the tissue substitute inserts employed in this work with an accuracy within 1.5%, as shown in section 3.1.It should be stressed that this accuracy level is achieved in an imaging study with an insert configuration including soft and hard tissue in the same scan, i.e. mimicking the real patient scanners to be performed for the clinical study.Different insert configurations were studied (results not shown), and no significant differences were obtained as the stability of both HU and MD values were highly reproducible across the whole FOV independently of the insert configurations.In this sense, our results are similar to those from Ainsley and Yeager (2014), stating that a simplified setup with only one image acquisition may be employed for CT scanner calibration, even in the case of DECT studies.Also, our accuracy results when employing the MD are 0.5% higher than those from Lehmann et al (1981), but this may be attributed to different factors such as the insert configuration or the inclusion of lung inserts.It should be noted, however, that this initial deviation is not fully transferred to the calibration and it can be partially mitigated by the final procedure.
The Z ef is a magnitude with different definitions in the bibliography, as shown in section 2.2.The definition employed for computation may have a non-negligible influence on the final computation of the SPR values, so attention should be paid to this matter.Also, as may be seen in section 3.2, the two definitions evaluated in this work show differences that may reach up to 10%.Regarding clinical practice, we think that it would be desirable to have a working definition for Z ef based solely on material compositions in order to avoid potential discrepancies due to tube spectrum variations.In particular, we have worked with a general spectrum for our scanner kindly provided by colleagues from another research group, but it is a well-known fact that spectrum variations may be found in the clinical practice after some interventions of the technical service, and detailed spectrum characterization is a complex measurement that is not usually performed routinely.Thus, these expected spectrum variations during the clinical operations potentially would yield a new characterization of the real spectrum and a new computation of the Z ef values for the calibration inserts, even when their composition is known and remains unchanged.For this reason, we decided to work with the general definition shown in equation (4) based on the known compositions to simplify the clinical workflow related to CT calibration.Nevertheless, we agree with previous efforts found in the bibliography devoted to providing a unified method to compute Z ef .
With regard to the new calibration method proposed for MD maps, it is noteworthy that two new magnitudes, X and Y (i.e.: f Z X ef 1 -( ) and f Y (Z ef )), are computed from the water and iodine densities (i.e.: ρ M and ρ N ), employing always both of them and treating them symmetrically.It was decided to proceed in this way to take profit of all the available information in the two maps.Also, for both methods, attention should be paid to the choice of fitting functions that may impact the whole method due to overfitting.In our case, we decided to not use exponents higher than three when using polynomials, as some undesired behavior was observed with fourth-degree polynomials for both methods.According to the work of Bourque et al (2014), polynomials up to the sixth degree can be employed when dealing with other CT scanners, but it should be considered that these fits may be influenced by the experimental data and the unknown particular spectrum itself, so we recommend to keep a good balance between the polynomial order and the adjustment of the extreme data in the working interval.In this sense, it is always a good starting point to evaluate the theory discussed in sections 2.4, 2.5, and 2.6 that explicitly states functional dependencies.Although these functional dependencies may show errors up to 1% according to previously published results (Alvarez and Macovski 1976), this knowledge may help choose appropriate functions in terms of monotonic functional behavior and form.Thus we agree on the approach to the problem of previously published efforts (Bourque et al 2014), but we recommend to evaluate with the particular experimental results obtained with every single CT scanner, since the lack of formal definition in equations (8)-(10) suggests some freedom to choose the experimental fitting functions and improve the experimental results.
Another point related to the new calibration method is the accuracy deviation found for some inserts such as those corresponding to lung and blood.However, it is shown in section 3.4 that the two calibration methods give rise to a similar level of accuracy for all the inserts.Thus, we consider that those deviations are mainly related to the particularities of these inserts instead of the calibration procedures themselves.It should be also noted that the Bourque method gives rise to more accurate Z ef estimates (roughly 1% of the mean value); while the e W r estimates show the same accuracy level.However, the uncertainty outputs of the new method are lower than those from Bourque et al (2014).To properly analyze these results, the methodology followed in our study (see section 2.1) should be accounted for as we have tried to replicate the acquisitions of DS-DECT equipment, but some differences may still be found due to the different spectra employed, although we conjecture that these spectra differences have a minor impact.However, we hypothesize that the higher uncertainties in DS-DECT  equipment may be due to the lack of correlation between the two acquisitions while, with the fast-switching technique, the correlations in the iodine and water maps may be employed to reduce the amount of noise in images and uncertainty (McCollough et al 2020).On the other hand, we consider that the difference in Z ef accuracy is due to the MD theory itself, but it has a negligible impact on the SPR estimates.It is also worth noting that both calibration procedures are quite independent of the inclusion or exclusion of the lung inserts in terms of mean AE.However, the potential impact on the maximum AEs for our method makes us recommend a similar study when dealing with one of these two calibration methods in clinical practice.
The two methods evaluated in this work have shown similar performance in terms of clinical applicability to provide better tissue characterization and SPR determination.They have been applied to the experimental data acquired in the same conditions to provide a comparison between both of them.However, it should be mentioned that our method should be the preferred one for CT scanners that provide MD maps, while the Bourque et al method should be employed for CT scanners that provide HU low and high images.The slight difference found in accuracy between the two methods may be explained by the MD theory itself, since deviations found in figure 3 are mitigated but not fully compensated by the DECT calibration procedure.The difference in precision is found mainly in figure 7 and may be attributed to different noise levels in the calibration data.On the one hand, the spatial and temporal correlation in the MD maps results in reduced noise in the VMIs, compared to the HU maps employed by the Bourque method (McCollough et al 2020).It should be also considered that the HU map acquisitions employed for the other protocol were performed by simulating an SS-DECT with two independent acquisitions that may be affected by different noise components that are propagated to the final results.In this sense, we think that some efforts should be performed in the medical physics community to provide an independent analysis showing the main differences found with the different technologies providing DECT images and the potential impact on final results.
Finally, we see that a global uncertainty analysis together with a biological tissue validation to determine the capabilities of these DECT calibrations in range uncertainty reduction is missing.Thus, the next steps in our research plan are to carry out these final end-to-end tests, perform full uncertainty analyzes, and implement the use of the maps in clinical practice.Also, some further investigations will be conducted to optimize the approach to DECT calibrations with these CT scanners.

Conclusion
In this paper, a new method to calibrate DECT scanners based on MD maps has been introduced and compared to a method based on HU low and high images.The equations of the method are straightforwardly derived from the MD theory, without additional hypothesis.Also, it has been shown how this new method yields similar accuracy with reduced experimental uncertainty results to previously published methods.Also, the performance of the new method is similar to the previous one by providing ρ e , Z ef , I, and SPR maps.Therefore, SPR-calibrated maps from DECT studies may be employed to reduce uncertainty in clinical proton and hadron therapy treatments after clinical validations.

Figure 1 .
Figure 1.(a) Density phantom by GAMMEX, used for calibration and evaluation; (b) Head phantom by CIRS, used for evaluation.

1-
(figure 5(a)), a sum of two exponential functions is employed, while in the case of f Y (figure 5(b)), a third-order polynomial function is employed for the fit.In this case, all the inserts have been included to perform the calibration.It is worth noting that the two lung inserts are the ones whose behavior deviates the most in f X and f Y characterization, together with a blood insert in the case of f Y .

Figure 2 .
Figure 2. (a) Tissue substitute inserts configuration in the GAMMEX phantom (b) iodine and (b) water maps.All different tissue substitutes contained in the phantom GAMMEX are highlighted with a blue circle in (a).

Figure 3 .
Figure 3. Residual errors after fitting experimental data to equation (14).The inserts are listed by their e W r , presented in table2.The residual errors are presented as the percent difference between the linear attenuation coefficients obtained from the SECT images and the fit using the equation coming from the MD theory.

Figure 4 .
Figure4.Differences between Z ef estimates obtained using the proposed method and those obtained using Bourque's.These differences are directly computed as the subtraction between both quantities.

Figure 5 .
Figure 5. Output of the calibration method, according to section 2.5 (a) f X , (b) f Y .

Figure 6 .
Figure 6.Residual percentage errors for the two DECT calibration methods relative to reference values from phantom certificate compositions (a) Z ef , (b) e W r .

Figure 7 .
Figure 7.Comparison of the output of the two methods for the map evaluations: (a) Z ef , (b) e W r , (c) I and (d) SPR.For the sake of clarity, we have added a small displacement in the horizontal coordinates of the URJC-QS values in the graphs to avoid a complete overlap.

Figure 8 .
Figure 8. Maps obtained: (a) Z ef bourque, (b) Z ef URJC-QS, (c) percent difference between a and b, (d) e W r bourque, (e) e W r URJC-QS, (f) percent difference between e and f, (g) I Bourque, (h) I URJC-QS, (i) percent difference between g and h, (j) SPR Bourque, (k) SPR URJC-QS and (l) percent difference between j and k.

Figure 9 .
Figure 9. ROIs employed to evaluate the performance of both methods shown on a HU map image.The yellow ROI stands for compact bone, blue for brain, red for soft tissue, and green for spongy bone.

Table 2 .
Elemental composition (in %), mass density ρ (in g/cm 3 ), and relative electron density to water e W r of the inserts that constitute the Gammex.Both blood inserts contain 0.10% of Fe.

Table 5 .
Results of the ROI analysis for the different tissues in the Proton Therapy Dosimetry Head Phantom.