Calculation of Stopping-Power Ratio from Multiple CT Numbers Using Photon-Counting CT System: Two- and Three-Parameter-Fitting Method

The two-parameter-fitting method (PFM) is commonly used to calculate the stopping-power ratio (SPR). This study proposes a new formalism: a three-PFM, which can be used in multiple spectral computed tomography (CT). Using a photon-counting CT system, seven rod-shaped samples of aluminium, graphite, and poly(methyl methacrylate) (PMMA), and four types of biological phantom materials were placed in a water-filled sample holder. The X-ray tube voltage and current were set at 150 kV and 40 μA respectively, and four CT images were obtained at four threshold settings. A semi-empirical correction method that corrects the difference between the CT values from the photon-counting CT images and theoretical values in each spectral region was also introduced. Both the two- and three-PFMs were used to calculate the effective atomic number and electron density from multiple CT numbers. The mean excitation energy was calculated via parameterisation with the effective atomic number, and the SPR was then calculated from the calculated electron density and mean excitation energy. Then, the SPRs from both methods were compared with the theoretical values. To estimate the noise level of the CT numbers obtained from the photon-counting CT, CT numbers, including noise, were simulated to evaluate the robustness of the aforementioned PFMs. For the aluminium and graphite, the maximum relative errors for the SPRs calculated using the two-PFM and three-PFM were 17.1% and 7.1%, respectively. For the PMMA and biological phantom materials, the maximum relative errors for the SPRs calculated using the two-PFM and three-PFM were 5.5% and 2.0%, respectively. It was concluded that the three-PFM, compared with the two-PFM, can yield SPRs that are closer to the theoretical values and is less affected by noise.


Introduction
An X-ray computed tomography (CT) image is a map of the photon linear attenuation coefficient, which is usually calculated as the sum of the cross-sectional contributions of the photoelectric effect, Compton scattering, and Rayleigh scattering. This linear attenuation multiple CT images. On the other hand, the two-PFM proposed by Torikoshi et al. [19] has been used more widely because of its mathematical simplicity. Two-PFM is a method of approximating the atomic cross section to calculate the linear attenuation coefficient with two terms, namely photoelectric absorption and combined Compton and coherent scattering, and calculating the parameters of each term to fit the theoretical value. For example, Ohno et al. [20] used two-PFM to calculate the effective atomic number and electron density from multiple CT images. However, because the current two-PFM ignores energy dependence while approximating the fitting parameter from the theoretical value, the method itself has a large error [3].
In this study, we present the inherent error of the aforementioned PFMs that fundamentally affect the results of the calculation of the effective atomic number and electron density. This study proposes a new formalism: a three-PFM, which can be used in multiple spectral CT. Both two-and three-PFMs were used to calculate the effective atomic number and electron density from multiple CT numbers using a photon-counting CT system. A semi-empirical correction method that corrects the differences between the CT values from the photon-counting CT images and the theoretical values in each spectral region was also introduced. The I-value was calculated via parameterisation with the effective atomic number, after which the SPR was calculated from the calculated electron density and the I-value. Then, the SPRs from both methods were compared with their corresponding theoretical values. To estimate the noise level of the CT numbers obtained from the photon-counting CT, CT numbers, including noise, were simulated to evaluate the robustness of the twoand three-PFMs.

Effective Atomic Number and Electron Density from the Linear Attenuation Coefficient
The X-ray linear attenuation coefficient µ can be related to the atomic cross-section as a function of energy E and a single element Z [1, 2] µ(E, where σ PE and σ KN are the photoelectric and Klein-Nishina cross sections, respectively; σ SCA is the correction term for the coherence and binding energy effects; and ρ, N A , and A are the mass density, Avogadro's number, and atomic mass, respectively. Each cross-section term can then be approximated via the PFM [1] µ(E, Z) = ρ e K PE (E)Z m + K KN (E) + K SCA (E)Z n , where ρ e is the electron density, defined by ρN A Z A , and K PE , K KN , and K SCA are coefficients for the photoelectric and Klein-Nishina terms and for the coherence and binding energy effects term, respectively. These terms can be determined via least-square fits to NIST attenuation coefficients from the NIST XCOM photon cross sections database [21]. Rutherford et al. [1] proposed values for the exponents m and n at 3.62 and 1.86, respectively. For mixtures, the electron density is calculated using Equation (3), and Z is represented by the effective atomic number (Z) [22] where ω i is the fractional weight (i.e., ∑ i ω i = 1). In general, α in Equation (4) is assigned values between 2.94 and 4 because it varies depending on the attenuation process and the materials being attenuated [22,23]. In our study, we adopted an α value of 3.3 for the theoretical calculations in Section 2.4.5, referring to Landry et al.'s [24] results [3,9,24]. Torikoshi et al. [19] proposed an expression of the linear attenuation coefficient as a photoelectric term and a scattering term based on Jackson and Hawkes' formulation [2] µ(E, Z) = ρ e Z 4 F(E, Z) + G(E, Z) , where ρ e Z 4 F(E, Z) is the photoelectric attenuation, and ρ e G(E, Z) is the combined Compton and coherent scattering term. Yang et al. [3] derived the F(E, Z) and G(E, Z) terms from fitting to attenuation cross-sections, referred from the NIST XCOM database, for each Z in a quadratic equation and interpolating over the entire Z range for several energies.
In this study, we employed the method expressed in Equation (5) and referred to it as two-PFM. Figure 1a shows the µ(E, Z)/ρ e calculated using Equation (5) from 30 keV to 150 keV obtained using the NIST XCOM database, and their corresponding NIST values and Figure 1b shows the relative errors. The functions F(E, Z) and G(E, Z) are obtained via quadratic fits of the photoelectric and scattering terms of the NIST attenuation coefficients.
where is the fractional weight (i.e., ∑ = 1). In general, in Equation (4) is assigned values between 2.94 and 4 because it varies depending on the attenuation process and the materials being attenuated [22,23]. In our study, we adopted an value of 3.3 for the theoretical calculations in Section 2.4.5, referring to Landry et al.'s [24] results [3,9,24].

Two-Parameter-Fitting Method
Torikoshi et al. [19] proposed an expression of the linear attenuation coefficient as a photoelectric term and a scattering term based on Jackson and Hawkes' formulation [2] ( , ) = ( , ) + ( , ) , where ( , ) is the photoelectric attenuation, and ( , ) is the combined Compton and coherent scattering term. Yang et al. [3] derived the ( , ) and ( , ) terms from fitting to attenuation cross-sections, referred from the NIST XCOM database, for each in a quadratic equation and interpolating over the entire range for several energies. In this study, we employed the method expressed in Equation (5) and referred to it as two-PFM. Figure 1a shows the ( , )/ calculated using Equation (5) from 30 keV to 150 keV obtained using the NIST XCOM database, and their corresponding NIST values and Figure 1b shows the relative errors. The functions ( , ) and ( , ) are obtained via quadratic fits of the photoelectric and scattering terms of the NIST attenuation coefficients. Equation (5) is a valid equation for a monochromatic X-ray beam. Using continuum energy spectrum X-rays, we adopted an average attenuation coefficient, which is approximated via integration of the photoelectric and scattering terms weighted with a normalised energy spectrum function [9]: The subscript indicates individual energy threshold levels in the spectral CT: = 1 to 4, and is the weighting function of the spectrum in each threshold Equation (5) is a valid equation for a monochromatic X-ray beam. Using continuum energy spectrum X-rays, we adopted an average attenuation coefficient, which is approximated via integration of the photoelectric and scattering terms weighted with a normalised energy spectrum function [9]: The subscript k indicates individual energy threshold levels in the spectral CT: k = 1 to 4, and ω E k is the weighting function of the spectrum in each threshold where Φ(E, E max ) represents the energy spectrum generated in an X-ray tube with a maximum energy E max . From Equation (6), the effective atomic number is calculated from four CT numbers obtained at four energy thresholds using Equation (8) where k represents each of the four energy threshold levels in the spectral CT; and F(E k , Z) and G(E k , Z) are the weighted sums of the spectra for F(E, Z) and G(E, Z), respectively. The electron density is then calculated using Equation (9), using the effective atomic number calculated in Equation (8): The Appendix A details the derivation process for Equations (8) and (9).

Three-Parameter-Fitting Method
In our study, we improved the accuracy of µ(E, Z)/ρ e by adding energy and atomic number dependency to each index of Equation (2) µ (10) where K PE (E, Z)Z PE(E,Z) , K I NC (E, Z)Z I NC(E,Z) , and K COH (E, Z)Z COH(E,Z) are terms for photoelectric absorption and incoherent and coherent scattering, respectively. In this model, each term has a new dependence on E and Z to improve the approximation accuracy more than that in Equation (2). Hereafter, we refer to this method as three-PFM. Figure 2a shows the µ(E, Z)/ρ e calculated using Equation (10) from 30 keV to 150 keV obtained using the NIST XCOM database, and their corresponding NIST values and Figure 2b shows the relative error values. We can see that, compared to two-PFM, as shown in Figure 1b, three-PFM can further reduce the relative errors with respect to the theoretical values.
where Φ( , ) represents the energy spectrum generated in an X-ray tube with a maximum energy . From Equation (6), the effective atomic number is calculated from four CT numbers obtained at four energy thresholds using Equation (8) where represents each of the four energy threshold levels in the spectral CT; and 〈 ( , )〉 and 〈 ( , )〉 are the weighted sums of the spectra for ( , ) and ( , ), respectively. The electron density is then calculated using Equation (9), using the effective atomic number calculated in Equation (8): The Appendix details the derivation process for Equations (8) and (9).

Three-Parameter-Fitting Method
In our study, we improved the accuracy of ( , )/ by adding energy and atomic number dependency to each index of Equation (2) where ( , ) ( , ) , ( , ) ( , ) , and ( , ) ( , ) are terms for photoelectric absorption and incoherent and coherent scattering, respectively. In this model, each term has a new dependence on and to improve the approximation accuracy more than that in Equation (2). Hereafter, we refer to this method as three-PFM. Figure 2a shows the ( , )/ calculated using Equation (10) from 30 keV to 150 keV obtained using the NIST XCOM database, and their corresponding NIST values and Figure 2b shows the relative error values. We can see that, compared to two-PFM, as shown in Figure Fig  Equation (10) can be approximated for polychromatic X-rays in the same manner as that used for Equation (6): Equation (10) can be approximated for polychromatic X-rays in the same manner as that used for Equation (6): From Equation (11), the effective atomic number can be numerically calculated from four CT numbers obtained at four energy thresholds, as follows: Sensors 2021, 21, 1215 6 of 20 The Appendix A details the derivation process for Equation (12). σ Weight E k=1 to 4 in Equation (12) represents the sum of the three terms in parentheses in Equation (11) at each threshold k, where each coefficient and exponent are weighted by the spectrum. The electron density was then calculated as follows: 2.1.3. Errors of the Two-and Three-PFM Figure 3a,b show the relative errors of the atomic number and electron density calculated using Equations (8) and (9) for the two-PFM and calculated using Equations (12) and (13) for the three-PFM, respectively, compared to the theoretical values. The calculations were performed under the assumption of monochromatic X-rays, for example, using values of 40, 60, 80, and 100 keV, which are the energies used in the photon-counting CT in Section 2.3. For the two-PFM, the calculated effective atomic number and electron density had relative errors within 1% and 0.35%, respectively, whereas, for the three-PFM, the calculated effective atomic number and electron density have relative errors, based on the theoretical values, within 0.5%, except for atomic numbers 2 and 3 (further discussed in Section 4), and 0.05%, respectively. However, because the effective atomic number of human tissue is between 6 and 15 [4], the errors for atomic numbers 2 and 3 will have a negligible effect on the calculation.
. (11) From Equation (11), the effective atomic number can be numerically calculated from four CT numbers obtained at four energy thresholds, as follows: The Appendix details the derivation process for Equation (12). in Equation (12) represents the sum of the three terms in parentheses in Equation (11) at each threshold , where each coefficient and exponent are weighted by the spectrum. The electron density was then calculated as follows: 2.1.3. Errors of the Two-and Three-PFM Figure 3a,b show the relative errors of the atomic number and electron density calculated using Equations (8) and (9) for the two-PFM and calculated using Equations (12) and (13) for the three-PFM, respectively, compared to the theoretical values. The calculations were performed under the assumption of monochromatic X-rays, for example, using values of 40, 60, 80, and 100 keV, which are the energies used in the photon-counting CT in Section 2.3. For the two-PFM, the calculated effective atomic number and electron density had relative errors within 1% and 0.35%, respectively, whereas, for the three-PFM, the calculated effective atomic number and electron density have relative errors, based on the theoretical values, within 0.5%, except for atomic numbers 2 and 3 (further discussed in Section 4), and 0.05%, respectively. However, because the effective atomic number of human tissue is between 6 and 15 [4], the errors for atomic numbers 2 and 3 will have a negligible effect on the calculation.
(a) (b) Figure 3. (a) shows the relative errors of the atomic number calculated via two-PFM using Equation (8) and three-PFM using Equation (12), respectively, compared to the theoretical values, and (b) shows the relative errors of the electron density calculated via two-PFM using Equation (9) and three-PFM using Equation (13), respectively, compared to the theoretical values.

Stopping-Power Ratio from I-Value Parameterisation
The SPR is a parameter used to determine the Bragg curve's position in ion beam dose calculation and can be calculated using the Bethe-Bloch formula [25], which can be approximated as follows [3,26] (8) and three-PFM using Equation (12), respectively, compared to the theoretical values, and (b) shows the relative errors of the electron density calculated via two-PFM using Equation (9) and three-PFM using Equation (13), respectively, compared to the theoretical values.

Stopping-Power Ratio from I-Value Parameterisation
The SPR is a parameter used to determine the Bragg curve's position in ion beam dose calculation and can be calculated using the Bethe-Bloch formula [25], which can be approximated as follows [3,26] where ρ e,m and ρ e,w are the electron densities of the medium and water, respectively; m e c 2 is the equivalent energy of the electron rest mass; I m and I w are the I-values of the medium and water, respectively, and β is the ratio of the ion velocity to the speed of light. The I-value of the mixture can be calculated using the Bragg additivity rule [27], as follows: In this study, the theoretical I-values of the elements were obtained from the Particle Data Group (PDG) database [28], which refers to ICRU Report 37 [29], and those of the mixtures were calculated using Equation (15).
Yang et al. [3] and Bourque et al. [4] have demonstrated that I-values can be parameterised as linear or polynomial functions of effective atomic numbers. Figure 4 shows the I-value for each element from the PDG and the I-values for human tissues, for which we refer to the effective atomic numbers and I-values in Table 1 of Bourque et al. [4], calculated using the Bragg additivity rule. The I-value for atomic number 1 was added to the mixture value in Figure 4, and least-squares fitting was applied to the fourth-order polynomial of the atomic number, denoted with a dashed line in Figure 4. This fourth-order polynomial fitting result was adopted in the SPR calculation. The I-values of the elements, obtained from PDG, were also fitted to the fourth-order polynomial, shown as a solid line in Figure 4.
where , and , are the electron densities of the medium and water, respectively; is the equivalent energy of the electron rest mass; and are the I-values of the medium and water, respectively, and is the ratio of the ion velocity to the speed of light. The I-value of the mixture can be calculated using the Bragg additivity rule [27], as follows: In this study, the theoretical I-values of the elements were obtained from the Particle Data Group (PDG) database [28], which refers to ICRU Report 37 [29], and those of the mixtures were calculated using Equation (15).
Yang et al. [3] and Bourque et al. [4] have demonstrated that I-values can be parameterised as linear or polynomial functions of effective atomic numbers. Figure 4 shows the I-value for each element from the PDG and the I-values for human tissues, for which we refer to the effective atomic numbers and I-values in Table 1 of Bourque et al. [4], calculated using the Bragg additivity rule. The I-value for atomic number 1 was added to the mixture value in Figure 4, and least-squares fitting was applied to the fourth-order polynomial of the atomic number, denoted with a dashed line in Figure 4. This fourth-order polynomial fitting result was adopted in the SPR calculation. The I-values of the elements, obtained from PDG, were also fitted to the fourth-order polynomial, shown as a solid line in Figure 4.    Figure 5 shows the experimental setup. The photon-counting CT system was housed in a lead-shielded box, which consisted of an X-ray tube with a tungsten target (L12161-07, Hamamatsu, Japan), and included a sample holder, collimator, detector module, and computer. The collimator was made of tungsten with a thickness of 7 mm and 64 circular converging holes, each with a diameter of 0.6 mm and spacings of 1 mm each. This design lowers the transmission of the scattered radiation of 150 keV energy to less than 1% [30]. The sample holder was made of poly(methyl methacrylate) (PMMA) and was cylindrically shaped, 30 mm in diameter, and 85 mm in height, with the capability to rotate 360 • . It can hold up to seven samples, each with a 5 mm diameter and 65 mm height. Water fills the space surrounding the sample rods [31]. An energy-differentiation-type 64-channel cadmium telluride (CdTe) radiation line sensor module (C10413, Hamamatsu, Japan) was used as the imaging system. The CdTe detector consists of 64 pixels of CdTe elements, and the size of each element is 0.8 mm (W) × 0.5 mm (H) × 5 mm (T). The elements are aligned with a pitch of 1 mm and a space of 0.2 mm between them. Each CdTe detector was connected to a 64-channel application-specific integrated circuit (ASIC) with 64 low-noise amplifiers. Signal pulses from each detector pixel were differentiated using five comparators according to their energy levels, where each threshold voltage could be controlled by user preference. The measured energy resolution at the full width at half maximum was 10% at 122 keV, which is the main peak of 57 Co. [32].  Figure 5 shows the experimental setup. The photon-counting CT system was housed in a lead-shielded box, which consisted of an X-ray tube with a tungsten target (L12161-07, Hamamatsu, Japan), and included a sample holder, collimator, detector module, and computer. The collimator was made of tungsten with a thickness of 7 mm and 64 circular converging holes, each with a diameter of 0.6 mm and spacings of 1 mm each. This design lowers the transmission of the scattered radiation of 150 keV energy to less than 1% [30]. The sample holder was made of poly(methyl methacrylate) (PMMA) and was cylindrically shaped, 30 mm in diameter, and 85 mm in height, with the capability to rotate 360°. It can hold up to seven samples, each with a 5 mm diameter and 65 mm height. Water fills the space surrounding the sample rods [31]. An energy-differentiation-type 64-channel cadmium telluride (CdTe) radiation line sensor module (C10413, Hamamatsu, Japan) was used as the imaging system. The CdTe detector consists of 64 pixels of CdTe elements, and the size of each element is 0.8 mm (W) × 0.5 mm (H) × 5 mm (T). The elements are aligned with a pitch of 1 mm and a space of 0.2 mm between them. Each CdTe detector was connected to a 64-channel application-specific integrated circuit (ASIC) with 64 lownoise amplifiers. Signal pulses from each detector pixel were differentiated using five comparators according to their energy levels, where each threshold voltage could be controlled by user preference. The measured energy resolution at the full width at half maximum was 10% at 122 keV, which is the main peak of 57 Co. [32].

Experimental Procedures 2.4.1. Photon-Counting CT Measurement
Seven rod-shaped samples were placed in a water-filled sample holder. The seven materials used for the samples were aluminium, graphite, PMMA, and four types of biological phantom materials (Kyoto Kagaku, Japan), that is, BE-T-10, BE-H-10, BE-N-10, and WD-3010, mimicking the properties of compact bone, cortical bone, internal bone, and water, respectively [31]. For the experimental conditions, the X-ray tube voltage and current were set at 150 kV and 40 µA, respectively, with focal spot sizes of 7 µm. All measurements were performed at room temperature, which was controlled between 20 ± 1 • C. The sample holder rotated 360 • at intervals of 2 • . The count source period of the detector for each projection was 100 ms and was repeated 100 times to reduce statistical uncertainty, which required 10 s per projection. Each CT number within each energy window was derived by subtracting the data taken above the threshold energy, k + 1, from those of k: 40 to 60 keV for energy window one, 60 to 80 keV for energy window two, 80 to 100 keV for energy window three, and 100 to 120 keV for energy window four. Each CT image was reconstructed using the fan-beam-based filtered-backprojection method with a Shepp-Logan filter [33], and ring artefacts were removed via a ring correction method in polar coordinates [34] and a three-step correction method introduced by Prell et al. [35].

Spectrum Measurement
To obtain the weighted sum from the spectrum to the coefficients and exponent values of each term in Equations (6) and (11), information on the X-ray energy spectrum ω E k at each threshold energy is required. After determining the weighting sum of each term, we could calculate the atomic number and electron density for the two-PFM and three-PFM. Because the diameter of the sample holder used in the CT measurement was 30 mm, a 30 mm-thick slab phantom was prepared [9] from 15 mm-thick PMMA and 15 mm-thick polypropylene and was attached to the X-ray tube. Fifty-four points (at approximately 2.1 keV intervals) were measured between 40 keV and 150 keV, with the X-ray tube voltage set to 150 kV. Figure 6 shows the spectra, measured using the CdTe line-sensor array, without and with the 30 mm-thick slab phantom. To obtain the spectral weights in Equations (6) and (11), the spectra measured using the 30 mm-thick slab phantom were interpolated using the cubic spline interpolation method to increase the energy resolution. The mean energy at each threshold can be calculated using Equation (16), where Φ(E) is the spectrum obtained from the experiment (Figure 6), and each sum of the spectra were obtained within each energy window.
The sample holder rotated 360° at intervals of 2°. The count source period of the detector for each projection was 100 ms and was repeated 100 times to reduce statistical uncertainty, which required 10 s per projection. Each CT number within each energy window was derived by subtracting the data taken above the threshold energy, k + 1, from those of k: 40 to 60 keV for energy window one, 60 to 80 keV for energy window two, 80 to 100 keV for energy window three, and 100 to 120 keV for energy window four. Each CT image was reconstructed using the fan-beam-based filtered-backprojection method with a Shepp-Logan filter [33], and ring artefacts were removed via a ring correction method in polar coordinates [34] and a three-step correction method introduced by Prell et al. [35].

Spectrum Measurement
To obtain the weighted sum from the spectrum to the coefficients and exponent values of each term in Equations (6) and (11), information on the X-ray energy spectrum at each threshold energy is required. After determining the weighting sum of each term, we could calculate the atomic number and electron density for the two-PFM and three-PFM. Because the diameter of the sample holder used in the CT measurement was 30 mm, a 30 mm-thick slab phantom was prepared [9] from 15 mm-thick PMMA and 15 mm-thick polypropylene and was attached to the X-ray tube. Fifty-four points (at approximately 2.1 keV intervals) were measured between 40 keV and 150 keV, with the X-ray tube voltage set to 150 kV. Figure 6 shows the spectra, measured using the CdTe line-sensor array, without and with the 30 mm-thick slab phantom. To obtain the spectral weights in Equations (6) and (11), the spectra measured using the 30 mm-thick slab phantom were interpolated using the cubic spline interpolation method to increase the energy resolution. The mean energy at each threshold can be calculated using Equation (16), where Φ( ) is the spectrum obtained from the experiment (Figure 6), and each sum of the spectra were obtained within each energy window. Figure 6. Spectra measured using CdTe line-sensor array, without and with 30 mm-thick slab phantom. Solid lines were calculated using the cubic spline interpolation method.

Semi-Empirical Correction Method for CT Values
The linear attenuation coefficients derived from the CT numbers using the photoncounting CT system were lower than those of the NIST values. This result can be attributed to the scattered X-rays and inter-pixel charge sharing of the CdTe detector [5,14,36]. Because the coefficients and exponent values of each term in Equations (6) and (11) were determined using the NIST values (by fitting), accurate values cannot be obtained in the effective atomic number and electron density calculations if the linear attenuation coefficients obtained from the CT number are different from the NIST values. Therefore, we applied a correction method that directly compares the linear attenuation coefficient from the CT number with the NIST value. To perform such a correction, each rod sample (5 mm in diameter) of PMMA, graphite, magnesium, aluminium, and titanium was separately placed at the centre of the water-filled sample holder to obtain CT images of each sample with respect to the four thresholds. Each linear attenuation coefficient derived from the CT image was then compared with the NIST value and fitted with a power function. Figure 7 compares the linear attenuation coefficients obtained from photon-counting CT and those calculated using Equation (17) with NIST values for each of the four energy windows. (17) of each sample with respect to the four thresholds. Each linear attenuation coefficient derived from the CT image was then compared with the NIST value and fitted with a power function. Figure 7 compares the linear attenuation coefficients obtained from photoncounting CT and those calculated using Equation (17) with NIST values for each of the four energy windows.
The CT numbers measured in each window were corrected to the expected NIST theoretical value via the plotted power function for each energy window in Figure 7. The corrected CT numbers were then used to calculate the effective atomic number and electron density using the two-and three-PFMs.

SPR Calculation Using Two-PFM and Three-PFM
With the four reconstructed images corrected using the method discussed in Section 2.4.3, we calculated the effective atomic number and electron density of each sample using the two-PFM (Equations (8) and (9)) and three-PFM (Equations (12) and (13)). The I-value was then calculated using the calculated effective atomic number, using the parameterisation method shown in Figure 4. Because human tissues, which are mixtures, are actually needed in a clinic, a fitting curve for the fourth-order polynomial fitting of the mixture, calculated using the Bragg additivity rule, was used for the conversion to the Ivalue. For the 290 MeV/u carbon-ion beam, SPR was calculated via the application of the electron density, and the I-value derived from the photon-counting CT images to Equation (14). The SPR depends on the energy of ions and has an error within 1% for the energy of the carbon ions used clinically, i.e., from 55 MeV/u to 430 MeV/u. In this study, the energy of 290 MeV/u carbon ions, which will be used to measure the Bragg curve in Section 2.4.5, was utilised for SPR calculation.  The CT numbers measured in each window were corrected to the expected NIST theoretical value via the plotted power function for each energy window in Figure 7. The corrected CT numbers were then used to calculate the effective atomic number and electron density using the two-and three-PFMs.

SPR Calculation Using Two-PFM and Three-PFM
With the four reconstructed images corrected using the method discussed in Section 2.4.3, we calculated the effective atomic number and electron density of each sample using the two-PFM (Equations (8) and (9)) and three-PFM (Equations (12) and (13)). The I-value was then calculated using the calculated effective atomic number, using the parameterisation method shown in Figure 4. Because human tissues, which are mixtures, are actually needed in a clinic, a fitting curve for the fourth-order polynomial fitting of the mixture, calculated using the Bragg additivity rule, was used for the conversion to the I-value. For the 290 MeV/u carbon-ion beam, SPR was calculated via the application of the electron density, and the I-value derived from the photon-counting CT images to Equation (14). The SPR depends on the energy of ions and has an error within 1% for the energy of the carbon ions used clinically, i.e., from 55 MeV/u to 430 MeV/u. In this study, the energy of 290 MeV/u carbon ions, which will be used to measure the Bragg curve in Section 2.4.5, was utilised for SPR calculation.

Theoretical Values of Effective Atomic Number, Electron Density, I-Value, and SPR
To compare the experimental values for the rod samples used in Section 2.4.1 with their corresponding theoretical values, the theoretical electron density and effective atomic number of the rod samples were calculated using Equations (3) and (4), respectively, and the I-values for the elements were obtained from the PDG database, and those for mixtures were calculated using Equation (15). The elemental compositions and mass densities of the biological phantom materials mentioned in Section 2.4.1 were obtained from the catalogue of Kyoto Kagaku [37] and are summarised in Table 1. The method of calculating the SPR by measuring the depth difference of the Bragg curve is conventionally referred to as water-equivalent length (WEL) or water-equivalent path length (WEPL) and has been adopted because of its reliability [38,39]. Here, we adopt the WEL calculation method for the theoretical value of SPR. This method is expressed as Equation (18), where R H,w is the depth of the Bragg curve from the distal H% in only water. R H,m is the depth of the Bragg curve from the distal H% in the water column when a slab material is placed between the water and the ionisation chamber; H is the height from the Bragg peak at the distal portion of the Bragg curve, which is generally represented as a percentage; and t m is the physical thickness of the slab material. In this study, the value of H was 80.
A height-adjustable water column was used to measure the SPR. Depth profiles of the Bragg curves from the 290 MeV/u carbon-ion beam were obtained in water only and for each slab absorber using a plane-parallel chamber (PTW, advanced Markus model 34045) [31]. A slab phantom of each material (aluminium, graphite, PMMA, and four biological phantom materials) was placed between water and an ionisation chamber in a water column. After each Bragg curve was obtained, each SPR was calculated using Equation (18). The values of effective atomic number, electron density, I-value, and SPR obtained from photon-counting CT were compared with their theoretical values to calculate their relative errors using Equation (19): For the values of the photon-counting CT images, the mean value and standard deviation were calculated by setting the region of interest (ROI) corresponding to the size of each rod sample (5 mm in diameter).

Random Noise Effect for Theoretical Linear Attenuation Coefficient
Compared to single-energy CT or dual-energy CT, photon-counting CT is more susceptible to image noise because only a portion of the photon is counted for each energy bin. Thus, it is necessary to test how the noise in the image affects the results [40,41]. CT numbers, including noise, were simulated to evaluate the robustness of the calculation method described in Sections 2.1 and 2.2. For noise, a Gaussian distribution was assumed, and a 6% noise level was considered (σ = 6%). Each noise was added to a total of 4096 theoretical linear attenuation coefficients obtained using Equation (17) for each of the seven samples mentioned in Section 2.4.1 with regard to four different thresholds. For 4096 noise-implemented linear attenuation coefficients, the effective atomic number and electron density were calculated for each sample using the two-PFM and three-PFM, and the I-value and SPR were calculated using the method described in Section 2.2. Figure 8 shows each reconstructed CT image of the seven samples, obtained via photon-counting CT system for the four different thresholds, corrected using the method discussed in Section 2.4.3. Ring artefacts are seen at thresholds one and four because of the high absorption rate at low energy and the low X-ray flux at the detector at high energy. In the case of BE-T-10, the sample itself cracked, and a dent portion occurred in the image.

Results
bin. Thus, it is necessary to test how the noise in the image affects the results [40,41]. CT numbers, including noise, were simulated to evaluate the robustness of the calculation method described in Sections 2.1 and 2.2. For noise, a Gaussian distribution was assumed, and a 6% noise level was considered ( = 6%). Each noise was added to a total of 4096 theoretical linear attenuation coefficients obtained using Equation (17) for each of the seven samples mentioned in Section 2.4.1 with regard to four different thresholds. For 4096 noise-implemented linear attenuation coefficients, the effective atomic number and electron density were calculated for each sample using the two-PFM and three-PFM, and the I-value and SPR were calculated using the method described in Section 2.2. Figure 8 shows each reconstructed CT image of the seven samples, obtained via photon-counting CT system for the four different thresholds, corrected using the method discussed in Section 2.4.3. Ring artefacts are seen at thresholds one and four because of the high absorption rate at low energy and the low X-ray flux at the detector at high energy. In the case of BE-T-10, the sample itself cracked, and a dent portion occurred in the image.  The I-values were parameterised using the fitting curve of the mixture in Figure 4, calculated using the Bragg additivity The values for each of the seven samples are plotted in Figure 11. The error of the SPR in Figure 11d tends to increase as the density of the sample increases. In the case of three-PFM, the overall relative error and standard deviation were less than those of two-PFM. Table 2 summarises the theoretical and experimental values and the relative errors for atomic number, electron density, I-value, and SPR. Figure 12 describes the theoretical values of the electron density and SPR in which 6% of noise level were added. As the electron density and SPR obtained from the experiment are within the error range of the theoretical value, it can be estimated that a noise level of about 6% may be included to the CT numbers obtained from the photon-counting CT. For three-PFM, the overall standard deviation was less than those of two-PFM. The values for each of the seven samples are plotted in Figure 11. The error of the SPR in Figure 11d tends to increase as the density of the sample increases. In the case of three-PFM, the overall relative error and standard deviation were less than those of two-PFM. Table 2 summarises the theoretical and experimental values and the relative errors for atomic number, electron density, I-value, and SPR. Figure 12 describes the theoretical values of the electron density and SPR in which 6% of noise level were added. As the electron density and SPR obtained from the experiment are within the error range of the theoretical value, it can be estimated that a noise level of about 6% may be included to the CT numbers obtained from the photon-counting CT. For three-PFM, the overall standard deviation was less than those of two-PFM.

Discussion
Most studies have adopted the two-PFM of Equation (5) for dual-energy CT calculations because it is possible to solve the effective atomic number and electron density clearly when there are two linear attenuation coefficients using two photon energies, assuring a certain degree of accuracy [3,9]. On the other hand, few publications have discussed the accuracy of the fitting method. In contrast, Yang et al. [3] reported that two-PFMs had errors within 1% of the effective atomic number and within 0.35% of the electron density. This result is in good agreement with the calculations shown in Figure 3a,b. Meanwhile, in the conventional approach for the three-PFM using Equation (2) because the exponents m and n of each term are fixed to constants, the calculated µ(E, Z)/ρ e has a somewhat large error, depending on the atomic number. We modified Equation (2) for each term to have a dependence on the atomic number, as suggested by Torikoshi et al. [19] for two-PFM, and significantly reduced errors in µ(E, Z)/ρ e (Figure 2). Therefore, the errors of the atomic number and electron density calculated using the three-PFM were within 0.5% (except for atomic numbers 2 and 3) and 0.05%, respectively, which were reduced by more than half of those of two-PFM. In Figure 3a, the relative errors of atomic numbers 2 and 3 are fairly high, which seems to be due to the polynomial fitting not fitting well for atomic numbers 2 and 3. For the range of atomic number 1 to 29, compared to the NIST value, atomic number 2 exhibited an error of 8%-10% for the coherent scattering term, and atomic number 3 exhibited an error of 0.6%-0.8% for the photoelectric attenuation. Thus, the relative errors of µ(E, Z)/ρ e for atomic numbers 2 and 3 appear to be greater than those for other atomic numbers (Figure 2b). However, because the atomic number of the sample used in this study was between 6 and 15, this error may not have had a significant effect on our calculations. Although both the two-and three-PFMs calculate the electron density after the effective atomic number is calculated, the calculated electron density does not seem to be significantly affected by the calculated effective atomic number (Figure 3).
To calculate the effective atomic number and electron density accurately using both two-PFM and three-PFM, the prerequisite is that the linear attenuation coefficient obtained from the CT number for each threshold should correspond to the NIST value. The linear attenuation coefficient measured using the CdTe detector; on the other hand, did not exhibit the theoretical value at low-energy (50-70 keV) and high-energy regions (120-150 keV). Therefore, the effective atomic number and electron density calculated using the parameters fitted by the NIST values would be lower or higher than the theoretical values. Nakashima et al. [36] and Matsumoto et al. [42] experimented with spectral CT using a CdTe detector and observed that the linear attenuation coefficients measured in low-energy regions were lower than the theoretical values because of scattered X-rays. They installed a molybdenum collimator and reduced the relative error of the linear attenuation coefficient by 16%. We thus fabricated a collimator from tungsten to compare the X-ray spectra and confirmed that the collimator is effective in low-energy regions, as reported by Nakashima et al. [36]. In high-energy regions, on the other hand, the linear attenuation coefficients still differ from the theoretical values. Miyajima [43] reported that the impurities of the detector disturb the drift of charge carriers and cause carrier trapping, which reduces the output pulse height and results in spectrum distortion. They corrected the spectral distortion due to transmission of primary X-rays, escape of secondary X-rays, and tailing, using the stripping method with the detector response function. However, it is considerably difficult to perform all corrections for each 64-detector element for different thresholds and for each material in a single sample holder at the same time. To address these problems, instead of solving the spectrum itself, we chose a semi-empirical method that theoretically corrects the experimental values by directly comparing the CT numbers with the NIST values weighted by the experimental X-ray spectrum for each threshold. In Figure 7, the difference between the theoretical and experimental linear attenuation coefficients increased with increasing material density and was the largest at threshold 1. We experimentally confirmed that the difference between the two linear attenuation coefficients was reduced when the energy width of the threshold was reduced, and we suspect that photon-counting CT was still affected by beam hardening.
In Figure 8, the photon-counting CT images show ring artefacts at threshold 1 (40-60 keV) and threshold 4 (100-120 keV). The artefacts on these images had a direct impact on the calculation of the effective atomic number and electron density. The artefact at threshold 1 seems to be related to the low number of photons because of the large attenuation coefficient, and the artefact at threshold 4 seems to be related to the detector response capable of absorbing high-energy photons. Using Equation (1), proposed by Tsutsui et al. [44], we calculated the X-ray absorbed fraction of our 5 mm-thick CdTe line-sensor array. CdTe absorbs 100% of low-energy photons, but the absorption starts to decrease from 100 keV to 85% at 150 keV. We set the threshold at four, that is, 120 keV, because severe artefacts started appearing on the image from 130 keV.
The effective atomic number calculations for graphite and aluminium demonstrated better results than the relative errors of 35% for graphite and 12% for aluminium in Nakashima et al. [36] and Zou et al. [45], who used two-PFM ( Table 2). As shown in Figure 11b, most values of electron density calculated using three-PFM, compared to those of two-PFM, are closer to the theoretical values. The calculation for the SPR image can be considered to be a deformation of the electron density image as much as the influence of the I-value map. However, the I-value estimation method was different for elements and mixtures ( Figure 4). Because we used the I-value from the mixture, the relative errors of two-PFM and three-PFM for aluminium were −36.7% and −38.6%, respectively. However, through the use of the I-value estimated from the element (fourth-order polynomial fitting), the errors could be reduced to −1.0% and −2.9%, respectively (Table 2). When the SPR is calculated from the element of fourth-order polynomial fitting, the relative errors of two-PFM and three-PFM for aluminium could be reduced by 4% to 12.7% and 2.9%, respectively. The Particle Data Group [28] mentioned that the I-value calculation using the Bragg additivity rule could be underestimated because electrons in a mixture are more tightly bound than in free elements. Hiraoka et al. [46] also showed that I-values calculated using the Bragg additivity rule tended to be lower than the I-values obtained via experimentation, with an uncertainty of ±5%.
A potential disadvantage of multiple-energy imaging is that the signal-to-noise ratio is degraded because the noise increases in proportion to the sum of the squares of the number of images [47]. However, the current study did not consider the evaluation of image quality and noise to see how realistic the SPR calculation results would be with the noise included. Instead, we estimated the noise level included in the reconstructed CT numbers in Figure 8 by adding noise to the theoretically calculated linear attenuation coefficient. As shown in Figure 12, it is thought that the three-PFM will be less affected by the noise contained in the CT numbers than the two-PFM. As Taasti et al. [41] mentioned, this noise level would be expected to be twice that of single-energy CT. We calculated the SPR using a photoncounting CT system, but there seems to be room for improvement in the image-based calculations, such as the inclusion of noise and artefacts of the reconstructed images, to accurately calculate the effective atomic number and electron density.

Conclusions
The SPR was calculated using a photon-counting CT system after multiple CT numbers were obtained at each threshold. The values for the four CT numbers obtained for the four energy thresholds were calibrated to their corresponding NIST values using a semiempirical correction method. To calculate the effective atomic number and electron density, we developed a new method, namely, three-PFM, applicable to multiple spectral CT to improve upon the conventional method, namely, two-PFM. Because the three-PFM fits the NIST value more accurately than the two-PFM, the effective atomic number and electron density calculated using the three-PFM exhibited less error than those calculated using the two-PFM. The effective atomic numbers and electron densities of the target materials were derived using each method. The effective atomic numbers were converted into I-values via parameterisation. The I-values and electron densities were applied to the Bethe-Bloch formula to calculate the SPRs of the target materials. The results of this study demonstrate that three-PFM, compared to two-PFM, calculated SPRs that were closer to the theoretical values and less affected by noise. Therefore, the proposed three-PFM can be potentially used in photon-counting CT, which is more easily affected by noise than single-energy CT or dual-energy CT. Acknowledgments: This work was supported by the Program for Cultivating Global Leaders in Heavy Ion Therapeutics and Engineering of the Ministry of Education, Culture, Sports, Science and Technology of Japan, grant to Gunma University. We would like to express our gratitude to the members of the Accelerator Engineering Corporation for their support.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
We rewrite Equation (6) to Equation (A1) µ k = ρ e Z 4 F(E k , Z) + G(E k , Z) , (k = Th1 to Th4), where k represents each of the four energy threshold levels in the spectral CT, and F(E k , Z) and G(E k , Z) are weighted sums from the X-ray spectra of Equation (6). If the electron density at each energy threshold is the same, then the aforementioned equation can be written as: We can then write Equation (A3) by binning each of the two equations of (A2): µ k Z 4 F(E k+1 , Z) + G(E k+1 , Z) − µ k+1 Z 4 F(E k , Z) + G(E k , Z) = 0, (k = Th1 to Th3).
Finally, Equation (A3) is summarised with respect to the effective atomic number and becomes Equation (A4): To calculate the electron density, if the atomic number is the same for each energy threshold, Equation (A1) can be written as: In a similar manner, we can write Equation (A6) by binning each of the two equations of (A5): F(E k+1 , Z) ( µ k − ρ e G(E k , Z) ) − F(E k , Z) µ k+2 − ρ e G(E k+2 , Z) = 0, (k = Th1 to Th3). (A6) We can then summarise Equation (A6) with respect to electron density and write it as Equation (A7):