Home

Search

Collections

Journals

About

Contact us

My IOPscience

A stoichiometric calibration method for dual energy computed tomography

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 Phys. Med. Biol. 59 2059 (http://iopscience.iop.org/0031-9155/59/8/2059) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 137.207.120.173 This content was downloaded on 28/05/2014 at 08:59

Please note that terms and conditions apply.

Institute of Physics and Engineering in Medicine Phys. Med. Biol. 59 (2014) 2059–2088

Physics in Medicine and Biology

doi:10.1088/0031-9155/59/8/2059

A stoichiometric calibration method for dual energy computed tomography Alexandra E Bourque 1,2 , Jean-Franc¸ois Carrier 2,3 and Hugo Bouchard 2,3,4 1

Medical Physics Unit, Montreal General Hospital (L5-113), McGill University, 1650 Cedar Avenue, Montreal, Quebec H3G 1A4, Canada 2 Centre hospitalier de l’Universit´e de Montr´eal (CHUM), 1560 Sherbrooke est, Montr´eal, Qu´ebec H2 L 4M1, Canada 3 D´epartement de physique, Universit´e de Montr´eal, Pavillon Roger-Gaudry (D-428), ´ 2900 Boulevard Edouard-Montpetit, Montr´eal, Qu´ebec H3T 1J4, Canada E-mail: [email protected] Received 31 July 2013, revised 20 December 2013 Accepted for publication 23 January 2014 Published 2 April 2014 Abstract

The accuracy of radiotherapy dose calculation relies crucially on patient composition data. The computed tomography (CT) calibration methods based on the stoichiometric calibration of Schneider et al (1996 Phys. Med. Biol. 41 111–24) are the most reliable to determine electron density (ED) with commercial single energy CT scanners. Along with the recent developments in dual energy CT (DECT) commercial scanners, several methods were published to determine ED and the effective atomic number (EAN) for polyenergetic beams without the need for CT calibration curves. This paper intends to show that with a rigorous definition of the EAN, the stoichiometric calibration method can be successfully adapted to DECT with significant accuracy improvements with respect to the literature without the need for spectrum measurements or empirical beam hardening corrections. Using a theoretical framework of ICRP human tissue compositions and the XCOM photon cross sections database, the revised stoichiometric calibration method yields Hounsfield unit (HU) predictions within less than ±1.3 HU of the theoretical HU calculated from XCOM data averaged over the spectra used (e.g., 80 kVp, 100 kVp, 140 kVp and 140/Sn kVp). A fit of mean excitation energy (I-value) data as a function of EAN is provided in order to determine the ion stopping power of human tissues from ED–EAN measurements. Analysis of the calibration phantom measurements with the Siemens SOMATOM Definition Flash dual source CT scanner shows that the present formalism yields mean absolute errors of 4 Present address: Acoustics and Ionizing Radiation Team, National Physical Laboratory, Hampton Road, Teddington TW11 0LW, UK.

0031-9155/14/082059+30$33.00

© 2014 Institute of Physics and Engineering in Medicine Printed in the UK & the USA 2059

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

(0.3 ± 0.4)% and (1.6 ± 2.0)% on ED and EAN, respectively. For ion therapy, the mean absolute errors for calibrated I-values and proton stopping powers (216 MeV) are (4.1 ± 2.7)% and (0.5 ± 0.4)%, respectively. In all clinical situations studied, the uncertainties in ion ranges in water for therapeutic energies are found to be less than 1.3 mm, 0.7 mm and 0.5 mm for protons, helium and carbon ions respectively, using a generic reconstruction algorithm (filtered back projection). With a more advanced method (sinogram affirmed iterative technique), the values become 1.0 mm, 0.5 mm and 0.4 mm for protons, helium and carbon ions, respectively. These results allow one to conclude that the present adaptation of the stoichiometric calibration yields a highly accurate method for characterizing tissue with DECT for ion beam therapy and potentially for photon beam therapy. Keywords: tissue characterization, proton therapy, computed tomography, dual energy, electron density, stopping power, mean excitation energy 1. Introduction In radiotherapy, computed tomography (CT) imaging is the gold standard method to obtain patient-specific data for treatment planning. To account for tissue heterogeneities, dose calculation engines require CT calibration curves which link, in each voxel, the Hounsfield unit (HU) to other physical properties such as the electron density (ED) or a property like stopping power (SP). To determine ED or SP maps using conventional single energy CT (SECT) scanners, methods based on the stoichiometric calibration proposed by Schneider et al (1996) yield more reliable results than the use of direct phantom-based CT calibrations (Schneider et al 1996, 2000, Guan et al 2002, Vanderstraeten et al 2007). For analytical types of dose calculation algorithms, such as photon convolution/superposition methods (Ahnesj¨o 1989) or proton range calculations (Schaffner and Pedroni 1998), SECT calibration curves (i.e., HU–ED or HU–SP relations) have sufficient information to perform these calculations. However, when using transport-based algorithms such as Monte Carlo techniques, the elemental compositions need to be known to establish medium cross sections (Chetty et al 2007, Paganetti et al 2008). Common tissue characterization methods using SECT calibration curves assume that medium compositions are discretizable into specific ranges of HU where either a fixed composition is used (Du Plessis et al 1998, Jiang and Paganetti 2004) or a mixture of two tissues is assumed (Schneider et al 2000). Although such methods can yield reasonable accuracy in determining ED (Paganetti 2012), uncertainty on HU measurements can cause tissue misallocation and yield errors on absorbed dose calculation of the order of 10% or more (Bedwani et al 2013). Moreover, while the nature of HU–ED (or HU–SP) relations requires performing multi-linear fits, subjectivity arises when it comes to defining the number of segments and their general trend (Thomas 1999, Verhaegen and Devic 2005). Hence it is recognized that SECT calibration curves limit the accuracy of characterizing tissue for treatment planning. In a recent study, Yang et al (2010) theoretically showed the superiority of dual energy CT (DECT) against standard clinical practice for determining medium properties relevant to dose calculations. While early developments in CT imaging proposed the use of DECT to extract several physical properties simultaneously (Hounsfield 1973, Alvarez and Macovski 1976, Rutherford et al 1976), recent technological developments in spectral CT (Schlomka et al 2008), DECT (Chandra and Langan 2011) and dual source CT (DSCT) (Flohr et al 2006, Petersilka et al 2008) are expected to improve tissue characterization for radiotherapy. 2060

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Several studies dealing with the extraction of the physical properties of tissues from DECT images were recently published in the literature. An approach dealing with monoenergetic synchrotron radiation was proposed by Torikoshi et al (2003), allowing the determination of ED and the effective atomic number (EAN) and yielding an overall accuracy around 1% when determining ED experimentally (Tsunoo et al 2008). Another method was also developed for monoenergetic synchrotron radiation in the same year by Kirby et al (2003) allowing the prediction of ED measurements with an accuracy of about 5%. The formalism of Torikoshi et al (2003) was adapted further by Bazalova et al (2008b) for commercial scanners which use polyenergetic photon spectra yielding an accuracy of 1.8% on average for ED and 2.8% on average for EAN, this using phantom measurements. The same formalism was used by Landry et al (2011) to compare experiments with Monte Carlo simulations of a CT scanner, yielding an agreement within ±5% both for ED and EAN. A method was published by Mahnken et al (2009) and demonstrated the ability to distinguish body fluids with a precision of about 0.02 g cm−3 and 0.1 in determining the mass density and EAN, respectively. A method was also proposed by Midgley (2011) which parametrized the attenuation coefficient as a function of ED and a compositional ratio (i.e., being a function of the EAN) yielding a theoretical accuracy of about 1–2% for ED. Finally, the method of Saito (2012) to establish a relation between HU (i.e., being the HU difference between two energies) and ED theoretically demonstrated an agreement within 0.7% and 2.5% for simulated and measured ED respectively. The mathematical robustness of the methods proposed in the literature to be used with commercial CT scanners relies on the accuracy to which cross sections are parametrized as a function of physical properties as well as the knowledge of the imaging photon spectrum. At the basis of the stoichiometric calibration of Schneider et al (1996), the cross section parametrization introduced by Jackson and Hawkes (1981) relies on Mayneord’s power law (Mayneord 1937, Spiers 1946) and is limited when it comes to modeling the photoelectric effect in human tissues containing high-Z media, such as bone (Williamson et al 2006) or the thyroid (Yang et al 2010). Other problems also arise with the attempt to consistently define the EAN using a power law, as the model is approximate and therefore the choice of the exponent is not clearly defined (Bazalova et al 2008b, Landry et al 2011, Saito 2012). One advantage in using an improved cross section parametrization (Torikoshi et al 2003, Kirby et al 2003, Midgley 2004) is to increase the accuracy to which one can model high-Z components. As for the knowledge of the photon spectrum, the method of Bazalova et al (2008a, 2008b), Landry et al (2011) is based on experimental spectrum measurements and empirical beam hardening corrections, while the method of Mahnken et al (2009) assumes a predetermined spectrum and neglects beam hardening effects. In both methods, either experimental measurements or reliable information from the manufacturer is required. The main idea of this study is to implement the stoichiometric calibration method of Schneider et al (1996) in DECT with the goal of making it suitable for commercial spectral CT, DSCT or DECT scanners. The theory is constructed as follows. The XCOM database (Berger et al 1998), which contains photon interaction cross sections for energies between 1 keV and 100 GeV, is parametrized as a function of the atomic number using a polynomial expansion (Kirby et al 2003, Midgley 2004), yielding more degrees of freedom than the model of Jackson and Hawkes (1981). Providing a rigorous definition of the EAN, the relation between HU, ED and EAN is theoretically established and expressible using linear matrix systems. The equation system is generalized for DECT, providing solutions for ED and EAN based on HU pairs. Using a XCOM-based theoretical framework, the developed mathematical formalism is thoroughly validated. A solution for the determination of the mean excitation energy (the I-value) from EAN is provided, making the method applicable to proton and ion therapy. Experimental measurements using a Siemens SOMATOM Definition Flash DSCT 2061

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

are performed to validate the method for clinical use. A thorough uncertainty analysis allows determining the mathematical robustness of the method and compares it to other methods in the literature. Following the discussion and conclusion on this study, a summary of the method is also provided to ease the clinical implementation. 2. Definitions and derived quantities 2.1. Basic quantities

The electron density ne (in cm−3 ) of a given medium is given by   Z ne = ρNA , A med where ρ is the mass density (in g cm−3 ), NA is Avogadro’s constant (in mol−1 ) and the ratio of number of electrons per molecular weight given by    Zi Z = ωi , A med Ai i

(1) Z A med

is

(2)

Ai being the atomic mass (in g/mol) and wi where for element i, Zi being the atomic number,  being the fractional weight (i.e., such that i wi = 1). The ED relative to the one of water is calculated as follows:   ρ AZ med ne   . = (3) ρe ≡ ne,w ρw AZ w For the medium of interest, the ED and the mass density are noted ne and ρ, respectively. For water, these quantities are ne,w and ρw , respectively, and AZ w is the total number of electrons per molecular weight. The I-value of the medium of interest is obtained from the Bragg additivity rule (ICRU 1984)  λi ln Ii , (4) ln I = i

where Ii is the I-value of the ith element and λi is the fraction of electrons from the ith element in the medium given by ω i Zi A λi =  Z  i

.

(5)

A med

From the I-value and the relative ED, one can determine the ion SP (in MeV/cm) using Bethe’s formula (Bischel 1972):    2me c2 β 2 k0 z2 2 −β , (6) S = ρe 2 ln β I(1 − β 2 ) with k0 = 0.170 45 MeV cm−1 , z the charge of the particle (e.g., z = 1 for protons), me c2 the rest mass energy of the electron and β the proton velocity relative to the speed of light in vacuum. 2.2. Electronic cross section parametrization

For elements, the linear attenuation coefficient of photons is given by the following formula (Jackson and Hawkes 1981):

2.2.1. Elements.

μ(Z) = ne σe (Z),

(7) 2062

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Electron Cross Section (barn)

60 80 kVp 100 kVp 140 kVp 140/Sn kVp

50 40 30 20 10 0

0

20

10

30

40

50

60

Atomic Number Figure 1. Electronic cross section averaged over the spectra of the Siemens SOMATOM

Definition Flash DSCT as a function of the atomic number.

where σe (in cm2 ) is the electronic cross section defined as the atomic cross section σa divided by the atomic number: σa (Z)  σa, j (Z) = , Z Z j=1 3

σe (Z) ≡

(8)

with the index j representing the interaction type, i.e., Rayleigh scattering, photoelectric effect and Compton scattering for j = 1, 2, 3 respectively. The cross sections (atomic or electronic) are averaged over a given photon spectrum as follows:

hνmax σE (hν, Z)s(hν) dhν, (9) σ (Z) ≡ 0

where σE (hν, Z) represents the cross section (atomic or electronic) for a specific photon −1 energy hν, s(hν) is the photon spectrum  hνmax (in MeV ) equal to the relative number of photons of energy hν in the spectrum (with 0 s(hν)dhν = 1), and σ (Z) represents the cross section (atomic or electronic) averaged over the spectrum and is implicitly dependent of s(hν). For a monoenergetic spectrum of energy E, one can use s(hν) = δ(hν − E ), the Dirac delta function centered at E, in equation (9) without loss of generality. Let us define the parametric electronic cross section σˆ e for any Z (integers or non-integers) as follows: σˆ e (Z) ≡

M 

am Z m−1 ,

(10)

m=1

with Z ∈ R and Z  1. The coefficients am are obtained using a least square fit on the electronic cross section data obtained from the XCOM photon cross section database (Berger et al 1998) averaged over a given spectrum over a range Z ∈ [Zmin , Zmax ]. M is the order of the fit which establishes the level of accuracy to which σˆ e (Z) = σe (Z) for Z integers. Figure 1 shows the behavior of electronic cross sections as a function of Z for all four spectra of the Siemens SOMATOM Definition Flash DSCT: 80 kVp, 100 kVp, 140 kVp and 140/Sn kVp (i.e., 140 kVp with Sn filter). It is worth noting that, for all spectra, the behavior of σe with Z allows σˆ e (Z) to be a bijective function in the range Z ∈ [1, 52]. 2063

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

2.2.2. Mixtures. For a given photon energy or spectrum, the mass attenuation coefficient of a mixture is given by the following relation (Jackson and Hawkes 1981)    μ μ = wi . (11) ρ med ρ i i

By analogy, one can show that the electronic cross section of a mixture is given by μ  ρ med Z σe,med = = λi σe (Zi ). NA A med i

(12)

2.3. Effective atomic number

As described in the definition of the HU, CT data provides an indirect measure of the ED and the electronic cross section. Using equations (3) and (7) as well as the definition of the HU, one expresses the following for elements: HU = 1000(ρe f (Z) − 1),

(13)

σe (Z) σe,w

the electronic cross section relative to water. While for elements, HU is a with f (Z) = function of the ED and the electronic cross section, a definition of EAN for media based on electronic cross sections is suitable for tissue characterization and thus be representative of what is obtained with the CT data. This reasoning yields the following definition. Let us define the domain Z ∈ [Zmin , Zmax ] where σˆ e (Z) is a bijective function, i.e., there exists a function σˆ e−1 such that Z = σˆ e−1 (σe ). The EAN of a medium, noted Zmed , is defined such that for a given photon spectrum, the parametric electronic cross section evaluated at Zmed equals the electronic cross section of the medium averaged over the spectrum σe,med . Mathematically, this means that there exists a unique Zmed such that σˆ e (Zmed ) = σe,med ,

(14)

Zmed ≡ σˆ e−1 (σe,med ).

(15)

and thus Since this definition can only yield a unique Zmed if σˆ e is a bijective function, one must limit the domain of Zmed when performing the fit of σˆ e (Z) for a given spectrum in order to avoid discontinuities with Z due to the fact that the energy of the K-shell (or any other shells) becomes significant with respect to the spectrum at higher Z. It is worth noting that with this general definition, Zmed depends not only on the elemental composition but also on the photon spectrum s(E ). As shown in figure 1, the definition of Zmed is valid for Zmed ∈ [1, 52] over the spectra used in this study. Figure 2 shows examples of calculated Zmed for four different tissue compositions using equation (15) as a function of energy including seven spectra from two different scanners (i.e., scanner 1: Siemens SOMATOM Definition Flash, scanner 2: Philips Gemini GXL). In the range of CT imaging energies, Zmed varies weakly with the photon spectrum (or kVp). For further uncertainty analysis, it is convenient to fix the value of EAN to a central value and to consider the variation of Zmed over relevant energies as a non statistical uncertainty component caused by the assumption that the EAN is independent of the photon spectrum. For each tissue composition, the central value of Zmed over the different spectra is defined as max{Zmed } + min{Zmed } , (16) Z¯ med = 2 2064

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Monoenergetic 80 kVp - Scanner 1 100 kVp - Scanner 1 140 kVp - Scanner 1 140/Sn kVp - Scanner 1 90 kVp - Scanner 2 120 kVp - Scanner 2 140 kVp - Scanner 2

6.40

6.36

Monoenergetic 80 kVp - Scanner 1 100 kVp - Scanner 1 140 kVp - Scanner 1 140/Sn kVp - Scanner 1 90 kVp - Scanner 2 120 kVp - Scanner 2 140 kVp - Scanner 2

13.54

Effective Atomic Number

Effective Atomic Number

6.44

13.53

13.52

13.51 6.32 30

40

60

50

70

80

90

70

80

Monoenergetic 80 kVp - Scanner 1 100 kVp - Scanner 1 140 kVp - Scanner 1 140/Sn kVp - Scanner 1 90 kVp - Scanner 2 120 kVp - Scanner 2 140 kVp - Scanner 2

7.54

7.52

Effective Atomic Number

Effective Atomic Number

7.56

90

8.10

Monoenergetic 80 kVp - Scanner 1 100 kVp - Scanner 1 140 kVp - Scanner 1 140/Sn kVp - Scanner 1 90 kVp - Scanner 2 120 kVp - Scanner 2 140 kVp - Scanner 2

8.05

8.00

7.50 40

50

60

70

100

Photon Energy (keV)

Photon Energy (keV)

80

90

40

Photon Energy (keV)

50

60

Photon Energy (keV)

Figure 2. EAN values calculated as a function of energy from the XCOM database for human tissues: adipose tissue (top left), cortical bone (top right), muscle (bottom left) and thyroid (bottom right). Note that scanner 1 refers to the Siemens SOMATOM Definition Flash and scanner 2 refers to the Philips Gemini GXL.

where the minimum and maximum values of calculated Zmed are taken over the energy range of photon spectra relevant to CT imaging. The variation of the EAN is defined as max{Zmed } − min{Zmed } Zmed = . (17) 2 Tables 1 and 2 provide fixed Zmed values with their variation (i.e., Z¯ med ± Zmed ) for the Gammex 467 phantom and human tissues recommended by ICRP (Snyder et al 1975) respectively using the spectra of the Siemens SOMATOM Definition Flash DSCT. The values computed with the spectra of the Philips Gemini GXL are typically included in the range of the Siemens SOMATOM Definition Flash data, with the exception of bones where the differences in fixed Zmed for both scanners are less than 0.008. 2.4. Mean excitation energy parametrization

Although there exists a strong correlation between the mean excitation energy and the EAN, a clear trend cannot be established for all media (ICRU 1984). As proposed by Yang et al (2010), the I-value can be parametrized as a function of the EAN for a restrained set of media (i.e., 34 standard human tissues) with sufficient accuracy with respect to tabulated data uncertainties. Based on the present definition of Zmed , one can determine the following empirical relationship for human tissues (taking Z = Zmed ): ⎧ for Z < Zmin ⎨e1 Z + e2 ˆ (18) I(Z) = e3 Z 5 + e4 Z 4 + e5 Z 3 + e6 Z 2 + e7 Z + e8 for Zmin  Z  Zmax ⎩ e9 Z + e10 for Z > Zmax 2065

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Mean Excitation energy (eV)

150

Fitted Bragg

100

50 10

5

15

Effective Atomic Number Figure 3. Parametrization of the ICRP mean excitation energy as a function of the EAN

defined in this paper for human tissues. Table 1. Summary of the Gammex 467 phantom media properties computed from mass densities and elemental compositions provided by the manufacturer (for reference use only).

Medium

ρe

Zmed

I (eV)

Iˆ (eV)

Water LN300 lung RMI-455 LN450 lung RMI-485 AP6 adipose tissue RMI-453 BR12 breast RMI-454 CT solid water RMI-451 LV1 RMI-482 SR2 brain RMI-481 CB2—30% CaCO3 RMI-484 CB2—50% CaCO3 RMI-480 SB3 bone, cortical RMI-450 B200 bone mineral RMI-487 IB3 inner bone RMI-456

1.00 0.29 0.48 0.93 0.96 0.99 1.06 1.05 1.28 1.47 1.70 1.10 1.09

7.45 ± 0.02 7.55 ± 0.03 7.52 ± 0.03 6.17 ± 0.02 6.87 ± 0.03 7.66 ± 0.03 7.66 ± 0.03 6.04 ± 0.03 10.76 ± 0.02 12.40 ± 0.01 13.51 ± 0.01 10.29 ± 0.03 10.28 ± 0.03

75.3 73.9 73.8 66.6 68.2 70.4 70.3 63.5 80.7 93.2 104.5 80.2 80.1

74.2 74.7 74.5 61.9 70.2 75.2 75.2 59.9 83.1 98.0 111.7 80.9 80.9

where e1 = 14.007 762, e2 = −24.414 214, e3 = −0.005 342, e4 = 0.207 079, e5 = −2.589 844, e6 = 8.339 473, e7 = 51.895 887, e8 = −219.722 173,e9 = 11.794 847, e10 = −47.707 141, Zmin = 6.26 and Zmax = 13.52. The central portion of the function described by equation (18) is obtained using a least square fit on the data described in table 2. The two other sections of the function are extrapolated linearly to make the relation continuous over Z and assuming that the slope of the function remains constant beyond Zmin and Zmax . It is worth noting that the main differences between the present fit and the method of Yang et al (2010) are the definition of EAN and the fact that the present method is continuous over Zmed . A statistical analysis of the fit residuals yields a mean uncertainty below ±1 eV (1σ level), which is negligible as compared to reported uncertainties on ICRU data (ICRU 1984). Tables 1 and 2 provide I-values determined with the Bragg additivity rule (ICRU 1984) as well as equation (18) for the Gammex 467 phantom and human tissues, respectively. For human tissues, the data are presented in figure 3. The elemental I-values used in the additivity 2066

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Table 2. Summary of human tissue properties computed from mass densities and elemental compositions recommended by the ICRP (Snyder et al 1975). The dataset used is consistent with the stoichiometric method of Schneider et al (1996).

Medium

ρe

Zmed

Adipose tissue Blood Brain Breast Cell nucleus Eye lens GI tract Heart Kidney Liver Lung (inflated) Lung (deflated) Lymph Muscle Ovary Pancreas Cartilage Red marrow Yellow marrow Skin Spleen Testis Skeleton-cortical bone Skeleton-cranium Skeleton-femur Skeleton-humerus Skeleton-mandible Skeleton-ribs (2,6) Skeleton-ribs (10) Skeleton-sacrum Skeleton-spongiosa Skeleton-vertebral cl (C4) Skeleton-vertebral cl (D6,L3) Thyroid

0.95 1.06 1.04 1.02 1.00 1.07 1.03 1.06 1.05 1.06 0.26 1.05 1.03 1.05 1.05 1.04 1.10 1.03 0.98 1.09 1.06 1.04 1.92 1.61 1.33 1.46 1.68 1.41 1.52 1.29 1.18 1.42 1.33 1.05

6.37 7.62 7.54 6.98 7.81 7.27 7.41 7.59 7.52 7.54 7.56 7.56 7.53 7.53 7.53 7.38 7.98 7.07 6.26 7.33 7.55 7.50 13.52 12.59 11.47 12.03 12.81 11.64 12.20 10.84 10.10 11.72 11.18 8.04

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.03 0.02 0.02 0.02 0.03 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.03 0.03 0.02 0.02 0.02 0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.02 0.03 0.01 0.02 0.05

I (eV)

Iˆ (eV)

64.8 75.2 73.9 70.3 75.5 74.3 74.3 74.8 74.7 74.8 75.2 75.2 75.5 74.6 75.0 73.4 78.0 69.2 63.9 73.7 75.0 74.7 112.0 100.0 87.0 92.5 102.7 90.7 95.9 84.6 78.4 91.3 87.0 74.7

64.8 75.0 74.6 71.1 75.7 73.2 74.0 74.9 74.5 74.7 74.7 74.7 74.6 74.6 74.6 73.8 76.3 71.8 63.3 73.5 74.7 74.5 111.8 100.3 88.2 93.7 103.0 89.7 95.6 83.6 80.2 90.5 85.9 76.4

rule are the ones recommended by the ICRU report 37 (ICRU 1984) and are shown in table 3. The differences between fitted and ICRU-recommended I-values shown in table 1 for the Gammex 467 phantom are much higher than the uncertainty on the fit (i.e., as high as 7 eV) due to the fact that the phantom media are not sufficiently tissue equivalent. 3. Methodology 3.1. A summary of the SECT stoichiometric calibration of Schneider et al

This section summarizes the SECT stoichiometric calibration proposed by Schneider et al (1996). Note that the following description is made in the scope of this paper and therefore a few differences are present, although the main idea of Schneider et al (1996) is respected. For a given energy, the attenuation coefficient (in cm−1 ) is modeled with the following parametrization introduced by Jackson and Hawkes (1981):   (19) μ = ne k ph z˜3.62 + kcoh zˆ1.86 + kKN . 2067

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Table 3. Recommended values of mean excitation energies to determine the I-values of mixtures using Bragg’s additivity rule (ICRU 1984).

Element symbol

Atomic number

I (eV)

H C N O Na Mg Si P S Cl K Ca Fe I

1 6 7 8 11 12 14 15 16 17 19 20 26 53

19.2 81.0 82.0 106.0 168.4 176.3 195.5 195.5 203.4 180.0 214.7 215.8 323.2 554.8

The values of z˜ and zˆ are calculated with the known chemical composition of each medium (or mixture) with the following power law additivity rule (Mayneord 1937, Spiers 1946) 1/m  λi Zim , (20) zm = i

such that z˜ ≡ z3.62 and zˆ ≡ z1.86 . From equation (19), the linear attenuation coefficient of each medium relative to water can be expressed as u≡

  μ ∗ ∗ , = ρe k∗ph z˜3.62 + kcoh zˆ1.86 + kKN μw

(21)

∗ ∗ and kKN are equal with ρe the ED of the medium relative to water and the coefficients k∗ph , kcoh to the coefficients k ph , kcoh and kKN , respectively, divided by μw . It is convenient to define the relative attenuation coefficient u as reduced Hounsfield Unit (or reduced HU), as it is related to HU as follows:

u=

HU + 1000 . 1000

(22)

∗ ∗ In order to find the values of k∗ph , kcoh and kKN , one needs to use a least square fit on experimental data measured with a calibration phantom (and known densities and elemental compositions). Once these scanner-specific constants are known for each energy spectrum, a calibration with human tissues is performed by computing the theoretical values of HU with equation (21) using recommended elemental compositions of body tissues (Snyder et al 1975, ICRU 1989, 1992, ICRP 2002). One should consider that in practice, each distinctive scanning protocol use its own combination of energy spectrum and algorithms, which include additional corrections during reconstruction which could affect the effective energy of the scan. Thus, these constants could also differ for each scanning protocol. The CT calibration curves are finally obtained by making a fit through human tissue data points which describe the relationship between theoretical HU (i.e., equations (21) and (22)) and the physical property or quantity in question such as ED (i.e., equation (3)) or SP (i.e., equation (6)) in a clinical situation. 2068

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

3.2. Adaptation of the SECT stoichiometric calibration

In order to develop a rigorous formalism for DECT based on a stoichiometric calibration, the method of Schneider et al (1996), which uses a three-variable parametrization with ρe , z˜ and zˆ, needs to be adapted to a two-variable cross section parametrization. Combining the parametrization of equation (10) with (7), the relationship of Schneider et al (1996) described at equation (21) can be rewritten by following for a given energy spectrum u = ρe f (Zmed ),

(23)

where f is a function given by f (Z) ≡

M 

bm Z m−1 .

(24)

m=1

Here M is a parameter establishing the level of accuracy of the relation. Note that, by definition, equation (24) is linked to equation (10) by bm ≡ μamw . In practice, the photon spectrum throughout the scanned phantom is not known accurately and, therefore, the coefficients bm are obtained using a least square fit of the following system based on the set of measurements HU = {HUn } with N distinctive media: u = Fb,

(25)

with u the N-elements array with un = m−1 , that is elements Fn,m = ρe,n Zmed,n ⎛ ρe,1 ρe,1 Zmed,1 ⎜ ⎜ ρe,2 ρe,2 Zmed,2 ⎜ F≡⎜ . .. ⎜ .. . ⎝ ρe,N ρe,N Zmed,N

HUn +1000 1000

...

M−1 ρe,1 Zmed,1

⎞

M−1 ⎟ ⎟ ρe,2 Zmed,2 ⎟ ⎟. .. ⎟ . ⎠ M−1 ρe,N Zmed,N

... ..

and F a matrix of dimensions N × M with

.

...

(26)

The least square solution of the coefficients bm is given by b ≈ bˆ = (Fcal T Fcal )−1 Fcal T ucal ,

(27)

with Fcal obtained with medium data of the calibration phantom and ucal the corresponding measurements. Finally, the theoretical values of HU for a given set of media are estimated with uˆ = ρe

M 

m−1 . bˆ m Zmed

(28)

m=1

3.3. DECT stoichiometric calibration

Let there be two sets of CT data (L = low kVp and H = high kVp), noted HUL = {HUL,1 , . . . , HUL,N } and HUH = {HUH,1 , . . . , HUH,N }, taken simultaneously with two distinctive energy spectra (or kVp). Using equation (22), for each voxel one can write HUL + 1000 1000 HUH + 1000 . uH = 1000 uL =

(29) 2069

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

To define a variable independent of the relative ED and depending on the EAN, the two following definitions are appropriate for constructing a DECT formalism. The variable

encompasses both definitions and is written as follows: ⎧ uL − uH ⎪ ⎪ for = DEI ⎨ uL + uH (30)

≡ uL ⎪ ⎪ ⎩ for = DER. uH Here DEI stands for dual energy index (Johnson 2011) while DER stands for dual energy ratio. Since is independent of the ED and assuming that Zmed is independent of the photon spectrum over the range used, one can write ≡ (Z) over a domain Z ∈ [Zmin , Zmax ] where

and Z are bijective. Therefore, one can write Z=

K 

ck k−1 ,

(31)

k=1

with K a parameter establishing the level of accuracy of the relationship. The equation can be written using the following matrix system: Z = c,

(32)

with Z the N-elements array with Zn = Zmed,n and  a matrix of dimensions N × K with elements n,k = nk−1 , that is ⎛

⎞

1

···

2 .. .

··· .. .

2K−1 ⎟ ⎟ ⎟ .. ⎟. . ⎠

1 N

···

NK−1

1

⎜1 ⎜ ≡⎜ ⎜ .. ⎝.

1K−1

(33)

The least square solution of the coefficients ck is then given by c ≈ cˆ = (cal T cal )−1 cal T Zcal ,

(34)

with cal and Zcal obtained respectively from measurements and medium data using the calibration phantom. For a post-calibration measurement, the estimator of the EAN is given by Zˆ med =

K 

cˆk k−1 ,

(35)

k=1

where is measured experimentally and K is a parameter establishing the level of accuracy of the relationship. Using the relationship established at equations (23)–(24), the estimator of the ED is defined using the following maximum likelihood estimation ρˆe ≡ 12 [ρˆe,L + ρˆe,H ],

(36)

2070

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

with ρˆe,L/H ≡ M m=1

uL/H , bˆ L/H,m Zˆ m−1

(37)

med

where the subscript L/H represents L or H. Here, the coefficients bˆ L,m and bˆ H,m are estimated for the two distinctive sets of CT data. 3.4. Uncertainty analysis

This section aims to derive analytic expressions for evaluating uncertainties of the estimated physical properties during the stoichiometric calibration. The present approach is particularly relevant since, despite that HU values can be assumed statistically independent, the mathematical formalism yields values of ED and EAN that are statistically correlated. The general idea of the uncertainty estimation presented here is to obtain an overall uncertainty on the extracted quantities, based on the statistical analysis of experimental results with respect to the parametrization described by equations (28) and (35)–(37). It is therefore assumed that all sources of uncertainties present in the method are included in this analysis. There are two situations of interest for evaluating uncertainties: (1) the stoichiometric calibration and (2) clinical situations. During the stoichiometric calibration, uncertainties are smaller since HU values are averaged over a region of interest constituted uniformly (typically a few cm3 ). In clinical situations, HU values are reported to single voxels (i.e., less than 1 mm3 ) and their uncertainty is dominated by noise, typically about 4–10 times the uncertainty on the average HU value. As it is shown in this section, the uncertainty on all physical properties studied is directly proportional to the uncertainty on HU and therefore the uncertainty on given properties in clinical situations are about 4–10 times higher than for the calibration. The uncertainty on the estimation of the ED and the EAN depends entirely on the random behavior of the sets uL and uH with respect to their expected values during the measurement (i.e., uL,meas and uH,meas ). One simple way to evaluate the uncertainty on these datasets is to assume that the variables behave normally (i.e., a Gaussian distribution) with an uncertainty being constant over the HU domain. This yields the following relationship to estimate the uncertainty on the reduced HU (Press et al 2002)   N 1  (uL/H,meas,i − uˆL/H,i )2 , (38) σL/H,meas ≈  ν i=1 3.4.1. Statistical uncertainties on experimental measurements.

where ν is the number of degrees of freedom of the fit (i.e., ν = N − M). The estimated values of reduced HU for the calibration phantom are evaluated, by applying equation (28) to the media of interest. However, for clinical measurements one can simply assume the uncertainty to be dominated by noise, that is σL/H,clin ≈ σL/H,noise .

(39)

To derive an analytic expression of the uncertainty of Zˆ med and ρˆe , one can make the following approximation. While the uncertainty on the estimation depends on fits (i.e., uˆL and uˆH ) and measurements (i.e., uL,meas and uH,meas ), the uncertainty contribution from the fit is small as compared to the measurement uncertainty when a large number of media is used in the stoichiometric calibration (e.g. N > 12, see Bouchard et al (2009)). Therefore, one can assume that the coefficients bˆ L,m , bˆ H,m , cˆL,m and cˆH,m have a negligible contribution to the uncertainty and this simplifies the variance analysis considerably. 2071

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Applying the rule of error propagation to equations (35) and (36), the uncertainties are given by the following relations (see appendix A):  ⎧   K k−2  ⎪ 2 (k − 1) c ˆ

  u2H σL2 + u2L σH2 ⎪ k k=1 ⎪ ⎪ ⎪ for = DEI ⎨ (uL + uH )2  σZˆ ≈  (40)  K ⎪ k−2  2 σ 2 + u2 σ 2 ⎪ (k − 1) c ˆ

u   k ⎪ k=1 H L L H ⎪ ⎪ for = DER ⎩ u2H and ⎧  2 2  ⎪ 1 ρˆe,L 2AuH ρˆe,H 2AuL ⎪ 2 ⎪ ⎪ − σL + + σH2 for = DEI ⎨2 uL (uL + uH )2 uH (uL + uH )2  σρˆe ≈  ⎪ ⎪ ρˆe,L 1 A 2 2 ρˆe,H AuL 2 2 ⎪ ⎪ ⎩ − σL + + 2 σH for = DER, 2 uL uH uH uH (41)

3.4.2. Statistical uncertainties on estimated properties.

where A≡

M  m=1

! (m − 1)

2 ρˆe,L

uL

bˆ L,m +

2 ρˆe,H

uH

" bˆ H,m

# m−2 Zˆ med

K 

# (k − 1)cˆk

k−2

. (42)

k=1

Note here that uL and uH are assumed statistically independent. However, ρˆe and Zˆ med are correlated and their covariance is given as follows (see appendix A). For = DEI, #  K  uH σL2 ρˆe,L 2AuH k−2 σρˆe ,Zˆ ≈ (k − 1)cˆk

− (uL + uH )2 uL (uL + uH )2 k=1  ρˆe,H 2AuL uL σH2 , (43) + − (uL + uH )2 uH (uL + uH )2 and for = DER, #   K  σL2 ρˆe,L uL σH2 ρˆe,H A AuL k−2 − . (44) (k − 1)cˆk

− + 2 σρˆe ,Zˆ ≈ 2uH uL uH uH 2u2H uH k=1 The uncertainty on the estimated I-value is given by the following rule of error propagation    ∂ Iˆ    σIˆ ≈   σZˆ , (45)  ∂Z  with the derivative of Iˆ estimated from equation (18) as follows: ⎧ e1 for Z < Zmin ⎨ ˆ ∂I for Zmin  Z  Zmax . = 5e3 Z 4 + 4e4 Z 3 + 3e5 Z 2 + 2e6 Z + e7 ⎩ ∂Z e9 for Z > Zmax The uncertainty on ion SP is given by the following relation (see appendix B):  # # #2 #   Sˆ 2 Sˆ ρe k0 z2 ∂ Iˆ ρe k0 z2 ∂ Iˆ  2 2 σρˆe + − σˆ − 2 σSˆ ≈ σρˆe ,Zˆ . Z ρe ρe Iˆ β 2 ∂Z Iˆ β 2 ∂Z

(46)

(47)

Finally, the uncertainty on the range of ions is approximately equal to (see appendix C):   $ σ2 ˆ  S, σRˆ ≈ Rx (48) 2 S 2072

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

where Rˆ is the estimated range (in mm), x is the average size of a CT-reconstructed voxel in the direction of the beam (in mm), σS2 is the average variance of the SP and S is the average SP, both quantities averaged over the ion trajectory. Since Zmed is by definition energy-dependent, an additional uncertainty component, which cannot be determined using statistical analysis, must be considered when applying the DECT formalism. In tables 1 and 2, values of Zmed and Zmed are given for different phantom media. These values are calculated using equations (16) and (17) over the spectra of the Siemens SOMATOM Definition Flash and the variation Zmed can be treated as a non statistical uncertainty due to the theoretical flaw of defining a fixed EAN. Therefore, the overall uncertainty (statistical and non statistical) of Zˆ med and ρˆe can be written as follows:  σZ,overall ≈ σ 2ˆ + Zmed 2 (49) ˆ 3.4.3. Non statistical uncertainties.

Z

and σρˆe ,overall ≈



σρ2ˆe + ρˆe2 .

(50)

Finally, ρˆe is evaluated as follows  % & % & M  M m−2 m−2  ˆ ˆ ˆ ˆ  u u (m − 1) b (m − 1) b Z Z L,m med H H,m med m=1 m=1 1 L  + ρˆe =   Zmed , % & % &2 2    2 M M m−1 m−1 ˆ L,m Zˆ ˆ H,m Zˆ   b b m=1 m=1 med med     ρˆe assuming Z ≈  ∂∂Zˆρˆe . med

(51)

med

4. Application of the methodology 4.1. Theoretical calculations

In order to evaluate the mathematical robustness of the formalism without the presence of artifacts such as scatter, noise and beam hardening observed in CT images, the stoichiometric calibration method is simulated using theoretical values of HU. This allows one to perform a self-consistency test of the method and compare its performance with the method of Schneider et al (1996). The theoretical framework is based on the XCOM database to calculate the attenuation coefficients of the phantom media being used during the calibration process, which are listed in table 4, using the spectra of 80, 100, 140 and 140/Sn kVp. While the XCOM database provides mass attenuation coefficients for elements for Z  100 at energies from 1 keV to 100 GeV (Berger et al 1998), the attenuation coefficients for mixtures and spectra are computed with equations (7), (9), (11) and (12). The linear attenuation coefficients of the phantom media are then used to estimate the parameters of the model as required by equations (27) and (34) (i.e., {bˆ m } and {cˆm }). Finally, the stoichiometric calibration is simulated for SECT and DECT applying the model to human tissue compositions (Snyder et al 1975), this using equation (28) for SECT and equations (35) and (36) for DECT. Note that the body tissues compositions used are consistent with the ones used by Schneider et al (1996). 4.2. Experimental measurements

The calibration phantom used is the Gammex 467 tissue characterization phantom, a 33 cm diameter Solid Water disk containing interchangeable rods of 12 human tissue substitutes 2073

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Table 4. Electron densities and elemental compositions of the Gammex 467 phantom media (for reference use only).

Phantom medium

ρe

ρ (g cm−3 )

H

C

N

O

Water LN300 lung LN450 lung AP6 adipose tissue BR12 breast CT solid water LV1 liver SR2 brain CB2—30% CaCO3 CB2—50% CaCO3 SB3 cortical bone B200 mineral bone IB3 inner bone

1.000 0.292 0.477 0.930 0.957 0.988 1.061 1.047 1.280 1.471 1.696 1.101 1.089

1.000 0.300 0.490 0.948 0.980 1.017 1.092 1.051 1.334 1.561 1.823 1.149 1.136

11.19 8.46 8.47 9.06 8.59 8.00 8.06 10.83 6.68 4.77 3.41 6.65 6.67

59.37 59.56 72.29 70.10 67.29 67.01 72.54 53.47 41.61 31.41 55.51 55.65

1.96 1.97 2.25 2.33 2.39 2.47 1.69 2.12 1.52 1.84 1.98 1.96

88.81 18.14 11.19 0.78 0.10 18.11 11.21 0.58 0.10 16.27 0.13 17.90 0.13 0.95 19.87 0.14 2.31 20.01 0.14 2.31 14.86 0.08 25.61 0.11 12.01 32.00 0.08 20.02 36.49 0.04 26.81 23.64 3.24 0.11 8.87 23.52 3.23 0.11 8.86

Mg

Si

P

Cl

Ca

Table 5. Machine parameters used during experimental acquisition. A 0.8 mm Sn filter

is present in the 140 kVp source (i.e., denoted 140/Sn kVp) for the first and second energy couples. Energy couple

1

2

3

Source A Exposure of source A Source B Exposure of source B

80 kVp 500 mAs 140/Sn kVp 193 mAs

100 kVp 250 mAs 140/Sn kVp 193 mAs

140 kVp 75 mAs 80 kVp 413 mAs

and 1 vial for water measurements. The elemental compositions of the rods are summarized in table 4. The relative electron densities range from 0.292 to 1.692 while the EANs range from 6.04 to 13.51. The phantom is scanned with a Siemens SOMATOM Definition Flash DSCT at three energy couples shown on table 5 with a pitch of 0.6 and a reconstruction kernel D30f in order to optimize the dual energy analysis. The parameters used during the scans are summarized in table 5. For each rod, the analyzed volume of interest (VOI) consists of ten 496 voxels with dimensions of 0.77 × 0.77 × 0.66 mm3 (total volume ≈ 3.63 cm3 ). For the stoichiometric calibration, the reported HU values are taken as the average HU over each VOI. For the noise analysis, the HU uncertainty is taken as the standard deviation of HU over each VOI. 5. Results 5.1. Self-consistency of the method

With respect to the formalism, there must exist a range of Z where the electronic cross section and (DEI or DER) are bijective with Z. Using XCOM data over all spectra of interest (i.e., 80, 100, 140 and 140/Sn kVp), theoretical values of

of each energy pair are plotted as a function of the atomic number of elements. As shown in figure 4 for all energy pairs, is bijective with Z for Z ∈ [1, 38]. Since the electronic cross section is bijective with Z over Z ∈ [1, 52] (see figure 1), the domain where both functions are bijective is therefore Z ∈ [1, 38]. The Zmed value for human tissues ranges from 6.26 to 13.52 (see table 2) and from 6.04 to 13.51 for the calibration phantom (see table 1).

5.1.1. Parametrization domain of Z .

2074

A E Bourque et al

0.5 0.4 0.3 0.2 0.1

80 - 140/Sn kVp 100 - 140/Sn kVp 80 - 140 kVp

0

Dual Energy Ratio (DER)

Dual Energy Index (DEI)

Phys. Med. Biol. 59 (2014) 2059

3

80 - 140/Sn kVp 100 - 140/Sn kVp 80 - 140 kVp

2.5

2

1.5

1 10

20

30

40

10

50

Atomic Number

20

30

40

50

Atomic Number

Figure 4. Behavior of the DEI (left) and DER (right) as a function of the atomic number for elements. The data is obtained using XCOM cross sections averaged over the spectra.

40 Schneider et al. This paper

60

Schneider et al. This paper

Residual HU

Residual HU

30 40

20

20

10

0 0 -10 0.5

1.0

1.5

0.5

Relative Electron Density

1.0

1.5

Relative Electron Density

Figure 5. A theoretical comparison of residual HU error between the model used by

Schneider et al (1996) and the present model. A 80 kVp spectrum (left) and a 140 kVp spectrum (right) are used to calculate theoretical HU values.

ˆ − HU, is 5.1.2. SECT stoichiometric calibration. The residual HU error, defined as HU evaluated for different tissue compositions using both (Schneider et al 1996) and the present SECT formalism. Figure 5 shows the residual HU error of human tissues obtained from the stoichiometric calibration simulated with the Gammex 467 tissue compositions. Note that the order of the fit used in the present method is 5 (i.e., M = 6 in equation (24)). The present method is more accurate in predicting the HU of human tissues than the one used by Schneider et al (1996) , particularly for the thyroid where the difference goes from 65 to 9.3 HU for 80 kVp and from 37 to 1.7 HU for 140 kVp. In addition, a systematic trend is observed in the bone region with the model used by Schneider et al (1996) , leading to errors up to about 11 HU and 7 HU for 80 and 140 kVp respectively, while the present method yields residual errors less than ±1.3 HU and ±0.8 HU for the same energies (except for the thyroid). Using the DECT calibration method for all three energy pairs, mean absolute errors of 0.026 % and 0.057% for the ED and the EAN are obtained using the calibration phantom. The adaptation of the calibration to human tissues with the energy pair 100–140/Sn kVp is found to yield optimal results and is shown in figure 6. Note that the error bars in figure 6 are calculated from analyzing the non statistical component of the uncertainty on Zmed , which by definition varies with energy. Although the thyroid

5.1.3. DECT stoichiometric calibration.

2075

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

2

0.6

Relative Error (%)

Relative Error (%)

0.8

0.4 0.2 0

1.5 1 0.5 0

-0.2 -0.5 0.5

1.0

1.5

6

Relative Electron Density

8

10

12

14

Effective Atomic Number

Figure 6. Relative errors in ED (left) and EAN (right) using the DECT formalism of theoretical HU values for 34 human tissues. The error bars are non statistical uncertainties, i.e., ρˆe is calculated with equation (51) and Zmed is extracted from table 2.

presents higher relative errors because of its content in iodine, all other tissues stay under approximately ±0.2% for both extracted quantities, with mean absolute errors of 0.08% and 0.17% for the ED and the EAN, respectively. Comparatively, the study of Yang et al (2010) reported these mean absolute errors to be 0.16% and 0.79% for the theoretical evaluation of the DECT formalism of Bazalova et al (2008b) with the 100 and 140 kVp spectra. The results shown in figure 6 are obtained with = DEI and K = 6 (as in equation (31)). It is worth noting that a similar level of accuracy can be achieved choosing = DER. As shown in figure 6, the 1σ -level error bars accurately represent probable deviations of the estimated parameter (i.e., ED or EAN) from its actual value, with the exception of the thyroid. As reported in several studies (Williamson et al 2006, Yang et al 2010), the thyroid is difficult to model due to its content in iodine which makes it highly sensitive to the photon spectrum due to the photoelectric effect dominance for low energies. As observed during the SECT calibration method (see figure 5), residual HU errors are present for the thyroid during the DECT stoichiometric calibration and cannot be explained by the non statistical errors shown in figure 6. This behavior is different from other human tissues, where the agreement is acceptable with respect to predicted uncertainties. The explanation of these discrepancies resides in the relationship between and Z. Indeed, the assumption that Zmed is independent of energy allows the establishment of the equivalence

≡ (Z) ↔ Z ≡ Z( ). In reality, this yields inconsistencies in performing the inversion of the function (Z), as done to obtain equation (31), and the thyroid is found to be the worst case. Figure 7 shows the non compliance of the thyroid to the –Z relation. Therefore one can conclude from the results of figure 6 that the accuracy of the ED and EAN measurements cannot be better than 0.6% and 1.5% respectively. For a better accuracy requirement, the thyroid could be handled as an organ that is not subject to the presented formalism. 5.1.4. The problem of the thyroid.

5.2. Experimental measurements 5.2.1. SECT measurements. While the model of Schneider et al (1996) remains an excellent method of calibration, the performance of the present formalism with experimental measurements is slightly superior in terms of accuracy. As shown in table 6, the overall uncertainty on HU estimated with equation (38) is systematically smaller with the present 2076

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Dual Energy Index (DEI)

0.04

0.03

0.02

0.01

ICRP tissues ICRP thyroid General trend

0 7.8

7.6

8.0

8.2

8.4

8.6

Atomic Number Figure 7. A theoretical comparison of the DEI between human tissues and the thyroid

for the energy couple 100–140/Sn kVp.

Table 6. HU uncertainties for the SECT experimental calibration (1σ level).

Tube potential (kVp)

Schneider et al

This paper

80 100 140 140/Sn

12 10 7 2

9 5 5 4

Table 7. Impact on the standard deviation of the HU for filtered back projection (FBP) and SAFIRE reconstruction processes. The third level of iteration of SAFIRE is used. LN300 (lung) AP6 (adipose tissue) Tube potential (kVp) FBP SAFIRE FBP SAFIRE

FBP

SAFIRE

FBP

SAFIRE

FBP

SAFIRE

80 100 140 140/Sn

56 39 52 32

39 26 39 23

50 36 48 26

34 25 32 20

74 46 56 38

46 33 43 29

47 37 42 31

26 24 28 22

57 34 44 34

36 23 31 22

CT solid water

LV1 (liver)

SB3 (cortical bone)

formalism compared to Schneider et al (1996), except for 140/Sn kVp, yielding overall uncertainties ranging from 4 to 9 HU over all spectra studied with the present method. These results are consistent with the manufacturer’s specification of ± 4 HU for water. While these levels of uncertainty for experimental measurements are significantly higher than for the theoretical framework (i.e., using XCOM data), they are mostly affected by the impact of noise and artifacts. In a more clinical situation where voxel information is used (rather than VOI-averaged), the noise component dominates the uncertainty and the type of reconstruction algorithm is a determining factor. In order to demonstrate the dependence of the uncertainty on the reconstruction algorithm, table 7 shows standard deviations of HU in several VOI using two different reconstruction algorithms: (1) filtered back projection (FBP), and (2) sinogram affirmed iterative (SAFIRE). Note that the reconstruction kernel (D30f) is kept consistent 2077

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059 1.2

1.02

1.1

ZDECT / Zmed

ρe_DECT / ρe

1.01

1.00

0.99

1.0

0.9 0.98

0.5

1.0

0.8

1.5

10

12

14

Effective Atomic Number 1.04

1.10 1.05

1.02

SDECT / Smed

IDECT / Imed

8

6

Relative Electron Density

1.00 0.95

1.00

0.98

0.90 0.85

0.96 60

80

100

1

Mean Excitation Energy (eV)

2

3

4

6

5

7

Proton Stopping Power (MeV/cm)

Figure 8. The ratios of calculated to known relative electron densities, EANs, mean excitation energies and proton SPs for the 100–140/Sn kVp energy couple. Note that the proton kinetic energy is 216 MeV. The reported uncertainties include statistical and non statistical components. Table 8. Experimental RMS errors (in %) of the calibration phantom’s physical properties determined experimentally.

Theoretical

Experimental

Energy couple

ρe

Zmed

I

S

ρe

Zmed

I

S

80–140/Sn kVp 100–140/Sn kVp 80–140 kVp

0.13 0.13 0.20

0.34 0.30 0.42

1.17 1.17 1.16

0.17 0.17 0.21

0.60 0.46 0.63

3.58 2.50 4.26

5.01 4.86 5.56

0.78 0.67 0.87

for both algorithms. A comparison between tables 6 and 7 demonstrates that the uncertainty increases by a factor of about 6–10 for water-like media with the FBP algorithm, while a factor of about 4–8 is observed with SAFIRE. Table 8 summarizes the RMS errors on all estimated physical properties for the calibration phantom with the present formalism using the theoretical framework and experimental measurements. More specifically, figure 8 shows the performance of the formalism to extract different physical properties of the Gammex 467 media with the optimal energy couple (i.e., 100–140/Sn kVp). The reported uncertainties include statistical and non statistical uncertainty components (i.e., related to the definition of Zmed ). The mean absolute error of the measured ED and EAN are respectively (0.3 ± 0.4)% and (1.6 ± 2.0)%. The mean absolute error of the measured I-value and the proton SP (216 MeV) are respectively

5.2.2. DECT stoichiometric calibration.

2078

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Proton ( z = 1, m = 1 ) Helium ( z = 2, m = 4 ) Carbon ( z = 6, m = 12 )

σR (mm)

0.16

0.12

0.08

0.04

0

200

100

300

Range (mm) Figure 9. Ideal uncertainty on the range of protons, helium ions and carbon ions in water as a function of the range of the beam. Clinically relevant uncertainties can be estimated accounting for CT noise, which magnifies σR shown in this graph by a factor of about 6 to 8. Table 9. Experimental uncertainties on ED and EAN for the 100–140/Sn kVp energy

couple. Medium

σρˆe

σZˆ

LN300 (lung) AP6 (adipose tissue) CT solid water LV1 (liver) SB3 (cortical bone)

0.006 0.001 0.006 0.006 0.003

1.5 0.3 0.4 0.4 0.01

(4.1 ± 2.7)% and (0.5 ± 0.4)% . Comparatively, the method of Bazalova et al yielded mean absolute errors of (1.8 ± 1.6)% and (2.8 ± 2.6)% for ED and EAN measurements respectively under similar conditions. Note that for higher I-values, the errors in figure 8 are significantly larger than the uncertainties since the estimator of I as a function of Z (i.e., see equation (18)) is performed using human tissue compositions only, i.e., with the absence of the calibration phantom media. The uncertainty on the estimation of the ED and the EAN is summarized in table 9 for a few calibration phantom media. The chosen media cover the range of the classical HU–ED curve: lung, adipose tissue, water and cortical bone. Uncertainties in ED range from 0.001 to 0.006 (relative to water), while the uncertainty on the EAN range from 0.01 to 1.5. While no clear trend exists for the uncertainty on the ED, the uncertainty on the EAN is high for low Zmed and low for high Zmed . This trend is also observed in figure 8. 5.3. Proton and ion range uncertainties

Using Bethe’s formula (see equation (6)), one can estimate the range of ions from ED and EAN values, converting EAN into an estimation of the I-value using equation (18). The uncertainty analysis is performed based on the equations in section 3.4. Figure 9 shows the trend of the uncertainty of the range in water based on the results of the stoichiometric calibration performed on experimental measurements. Three particles are studied: protons, helium ions 2079

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

and carbon ions, with both ions fully-ionized. The energy of the ions is chosen such that the maximum range is 30 cm and the continuous slowing down approximation (CSDA) is used to estimate the range numerically. Using the uncertainty levels reported during the stoichiometric calibration (i.e., table 6), the ideal uncertainty on the range is less than 0.16 mm for protons, less than 0.09 mm for helium ions and less than 0.06 mm for carbon ions. Despite these uncertainty levels being clinically unrealistic since the presence of noise is neglected, they can be used as a basis to estimate uncertainties in clinical situations. Using the ratio between results in table 7 and table 6 for solid water at 100 kVp and 140/Sn kVp, the HU uncertainties in clinical situations are higher by a factor of about 8 and 6 times than during the calibration, for the FBP and SAFIRE algorithms, respectively. As discussed in section 3.4, one can get a good estimation of the uncertainty on physical properties by multiplying the ideal uncertainty obtained from stoichiometric calibration measurements with these factors. This yields to an estimation on the clinical range uncertainty with the FBP algorithm (using factor of 8) of less than about 1.3 mm, 0.7 mm and 0.5 mm for protons, helium ions and carbon ions, respectively. For the SAFIRE algorithm (using factor of 6), the clinical range uncertainty is less than about 1.0 mm, 0.5 mm and 0.4 mm for protons, helium ions and carbon ions, respectively. The estimations for protons are significantly lower than the 1–3 mm reported by Schaffner and Pedroni (1998), which is based on the method of Schneider et al (1996) using SECT.

6. Discussion and conclusion This paper provides an accurate formalism for characterizing tissue properties relevant to radiotherapy from DECT images. The method is based on a novel definition of EAN used to adapt the stoichiometric method of Schneider et al (1996) for SECT. A theoretical calculation, based on attenuation coefficient calculation using phantom and human tissue compositions, demonstrates that the present adaptation improves the accuracy of the method in the bone region as compared to the method of Schneider et al (1996) (i.e., residual errors within ±1.3 HU compared to errors up to 7–11 HU). For the thyroid, comparing the same methods also shows improvements (i.e., residual errors up to 9.3 HU compared to 37–65 HU). The adapted stoichiometric method is extended to DECT and allows determining ED and EAN simultaneously from DECT images. The novel definition of EAN is also used to obtain a new parametrization of the mean excitation energy (I-value) for human tissues, allowing to determine SP from EAN and ED. Results comparison demonstrate a significantly higher accuracy in determining ED and EAN from experimental measurements than the method of Bazalova et al (2008b) (i.e., mean error of 0.3% and 1.6% versus 1.8% and 2.8%, for ED and EAN respectively, using the consistent calibration phantom manufacturer and the 100–140 kVp energy pair). While the comparison of the results with the method of Bazalova et al (2008b) conclusively demonstrates the higher accuracy of the present formalism, this could potentially be explained by several differences in the methods. While both are mathematically similar (i.e., the functions F (Z) and G(Z) of Bazalova are also polynomials fitted on XCOM-based data), the present formalism finds the advantage of using a more rigorous definition of EAN than Mayneord’s, which seems to partially explain the better accuracy. However, one can argue that a fair comparison between both methods depends on the accuracy of composition data known at the time of the study, which relies on the manufacturer, and also on the reconstruction algorithm being used. In such case, an algorithm in which correction for beam hardening effects is 2080

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

superior could explain why the present method can achieve better accuracy without empirical beam hardening corrections, such as proposed by Bazalova et al (2008b). Nonetheless, the present method benefits from the same reliability as the stoichiometric calibration method of Schneider et al (1996), since the constants of the polynomial functions used in the formalism are determined experimentally. This seems to yield a better agreement with measurements than spectrum-based methods. Furthermore, it makes the procedure suitable for clinical implementation since the access to accurate spectral information can be difficult and experimental calibration is typically more reliable than theoretical calculations. While one can argue that the present method depends on the knowledge of the spectrum through initial calculations of EAN, it is shown that spectrum information does not have a critical impact on the value of Zmed (as shown in figure 2 and tables 1 and 2). Indeed, one could use software (or freeware) to estimate the spectrum shape from a basic description of the x-ray tube, as the anode angle and the target medium (Poludniowski et al 2009). A thorough uncertainty analysis of the present method is also proposed to obtain an overall uncertainty on the quantities determined experimentally, i.e., HU values, based on the adapted stoichiometric method. The uncertainty on the quantity of interest, such as ED, EAN, SP and ion range can then be estimated for clinical situations. The approach is fundamental to account for statistical correlations between ED and EAN which allows one to obtain an accurate uncertainty estimation on SP and ion range. These correlations are present due to the nature of the DECT formalism, even if one makes the assumption that HU values are statistically independent due to the high noise content. Based on the experimental measurements performed with the Gammex 467 phantom, the uncertainty on proton range calculations performed with the DECT formalism are found within less than 1.0–1.3 mm for all therapeutic energies studied and for both reconstruction algorithms used, with CSDA ranges up to 30 cm in water. This assumes the noise to yield uncertainties on single-voxel HU values up to 6–8 times higher than the uncertainty on HU values averaged over ten 496 voxels, which is being used during the stoichiometric calibration. For helium and carbon ions, the same analysis yields range uncertainties below 0.5–0.7 mm and 0.4–0.6 mm, respectively. While it was suggested by Yang et al (2010) that DECT should improve the accuracy of tissue characterization in comparison with SECT, the present method shows clear improvements in predicting ED and ion ranges. Proton range uncertainty based on CT data alone are almost three times more precise than the ones predicted by SECT, i.e, as reported by Schaffner and Pedroni (1998) (i.e, uncertainties up to 1.0–1.3 mm versus 1–3 mm). This study yields the same positive conclusion for determining ED with DECT versus SECT. On one hand, a simple comparison in determining water ED can be achieved by comparing the uncertainty in determining water ED, (0.3% on average versus 4–9 HU for SECT which is approximately equivalent to 0.4–0.9% in ED). On the other hand, systematic errors caused by the nature of HU–ED curves (i.e., not being bijective) could yield errors of the order of about 2–3% in ED in the soft tissue region (Yang et al 2010). A more thorough study could be performed to quantify this comparison more accurately. Despite the mathematical robustness of the method, clinical tissue characterization could be improved in several ways. A major source of error in HU data can be beam hardening effects. Although vendors typically provide reconstruction algorithm performing acceptable corrections for that matter, they are not perfect and they are typically proprietary. Therefore it is difficult to characterize the effect of errors potentially produced by such beam hardening correction algorithms on the phantom. In this paper, the strategy used is to assume that the correction is ideal on the first order and that the errors caused by beam hardening effects should be included in the overall estimate of the uncertainty. While this assumption is matter 2081

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

to discuss, the present results are shown to be more accurate than the results of Bazalova et al (2008b), which include an additional beam hardening correction. It is unclear if the reconstruction algorithm (SAFIRE) plays an important role in this study. Nonetheless, it is expected that the accuracy of tissue characterization depends on the accuracy to which reconstruction algorithms handle beam hardening effects. Another source of error in the method is in the knowledge of the mass density and the elemental compositions of the calibration phantom. As shown by Yang et al (2010), errors in mass density and elemental compositions can significantly affect the accuracy of calculation in radiotherapy. The present strategy is to assume ideal knowledge a priori while evaluating an overall estimate of the uncertainty, which should include these errors. However, improvements are expected by investigating the accuracy of the mass density and elemental compositions of the calibration phantom media. Finally, one important consideration in the clinical accuracy of the method relies on the accuracy to which human tissue elemental compositions are known. In this study, data from ICRP report 23 (Snyder et al 1975) was used to be consistent when comparing with the method of Schneider et al (1996) . However, other sources of data exists and could potentially impact the results in this paper (ICRU 1989, 1992, ICRP 2002). Nonetheless, this paper provides a valid method for any data and therefore a thorough investigation on the impact of choosing data for body compositions is left for future work. 7. Summary of the method To implement the dual energy computed tomography stoichiometric calibration method accurately, here is a list of steps to follow: 1. If the Gammex 467 phantom is used, use Zmed values reported in table 1. Otherwise, choose a phantom with tissue substitutes of known compositions and physical densities, and calculate the value of Zmed for every medium according to the definition of equation (15). Remember the importance of choosing the compositions to be such that the ranges of Zmed and ρe cover those of human tissue compositions, such as the ones reported in table 2. Note that an ideal calculation of Zmed requires the knowledge of all CT spectra used. However, because of the low sensitivity of Zmed with energy, generic spectra covering the energy range could be used. 2. With the scanner used for radiotherapy treatment planning, scan the tissue substitutes with a given energy pair to obtain the pairs of Hounsfield unit (HU) values and convert them into pairs of reduced HU values (i.e., uL and uH ) with equation (29). 3. Determine (either dual energy index (DEI) or dual energy ratio (DER)) of every tissue substitute using equation (30). Note that either DEI or DER yields the same accuracy. Find the least square solution of the coefficients ck with the matrix system of equation (32), given by equation (34). The estimator of the effective atomic number in clinical situations is then calculated as in equation (35). 4. Perform the single energy computed tomography stoichiometric calibration for both energies using the series of reduced HU for both energies (i.e., uL and uH ). Find the least square solution of the coefficient bm with the matrix system of equation (25), given by equation (27). The estimator of the relative electron density is then calculated with equations (36) and (37). 5. For stopping power estimations, find the estimator of the mean excitation energy (i.e., the I-value) from Zˆ med with the relation stated by equation (18). 6. To perform the variance analysis, follow the equations stated in section 3.4. 2082

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Acknowledgments The authors are grateful to Siemens and Philips for kindly providing spectrum information. The authors are also grateful to Gammex for kindly providing media information. We wish to acknowledge our colleague Danis Blais from CHUM for his involvement during experimental measurements. We also wish to acknowledge Julia Snaith and Rebecca Nutbrown from NPL for fruitful comments on the paper. Finally, one author (AEB) wishes to acknowledge the financial support by the Minist`ere de la sant´e et des services sociaux du Qu´ebec. Appendix A. Derivation of the uncertainty on the effective atomic number and the electron density During these derivations, it is implicit that uL = uL,med , uH = uH,med and = meas . Using the rule of error propagation on equation (35), the variance of Zˆ med is given by (here, we take Zˆ = Zˆ med to ease the notation) ˆ σZ2ˆ ≡ VAR(Z) #2 ∂ Zˆ VAR( ) ≈ ∂

#2  #    K  ∂ 2 2 ∂ 2 2 k−2 ≈ (k − 1)cˆk

σL + σH . (A.1) ∂uL ∂uH k=1 with σL2 = VAR(uL ) and σH2 = VAR(uH ). Since stands for two definitions (i.e., DEI or DER, see equation (30)), one writes ⎧ 2uH ⎪ ⎪ for = DEI ⎨ ∂

(uL + uH )2 = (A.2) 1 ⎪ ∂uL ⎪ ⎩ for = DER uH and ⎧ 2uL ⎪ ⎪ for = DEI ⎨− ∂

(uL + uH )2 = . (A.3) uL ⎪ ∂uH ⎪ for = DER − 2 ⎩ uH Therefore equations (A.1) become  ⎧   K k−2  ⎪ 2 (k − 1) c ˆ

  u2H σL2 + u2L σH2 ⎪ k k=1 ⎪ ⎪ ⎪ for = DEI ⎨ (uL + uH )2   σZˆ ≈ . (A.4) ⎪  K (k − 1)cˆ k−2  u2 σ 2 + u2 σ 2 ⎪ k ⎪ k=1 H L L H ⎪ ⎪ for = DER ⎩ u2H The variance of ρˆe can then be written as follows: σρ2ˆe ≡ VAR(ρˆe )   ∂ ρˆe 2 ∂ ρˆe 2 ≈ VAR(uL ) + VAR(uH ) ∂uL ∂uH   ∂ ρˆe,H 2 2 1 ∂ ρˆe,L ∂ ρˆe,H 2 2 1 ∂ ρˆe,L + σL + + σH . ≈ 4 ∂uL ∂uL 4 ∂uH ∂uH 2083

(A.5)

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

The derivative of ρˆe,L/H with respect to uL is given by % & M ∂ Zˆ ∂uL/H uL/H (m − 1)bˆ L/H,m Zˆ m−2 ∂u m=1 ∂ ρˆe,L/H L ∂uL = M − % & 2  m−1 ˆ ˆ ∂uL M ˆ ˆ m−1 m=1 bL/H,m Z m=1 bL/H,m Z % & % & M K k−2 ∂

∂uL/H ˆ L/H,m Zˆ m−2 uL/H (m−1) b (k−1) c ˆ

k m=1 k=1 ∂uL ∂uL = M − % &2 ˆbL/H,m Zˆ m−1 M ˆ ˆ m−1 m=1 m=1 bL/H,m Z 2 ρˆe,L/H ρˆe,L/H ∂uL/H ∂

− BL/H , uL/H ∂uL uL/H ∂uL given by

= with BL/H

BL/H ≡

M 

(A.6) #

(m − 1)bˆ L/H,m Z

ˆ m−2

m=1

K 

# (k − 1)cˆk

k−2

.

(A.7)

k=1

By symmetry with equation (A.6), the derivative of ρˆe,L/H with respect to uH is given by 2 ρˆe,L/H ρˆe,L/H ∂uL/H ∂ ρˆe,L/H ∂

= − BL/H . ∂uH uL/H ∂uH uL/H ∂uH Combining equation (A.5) with equations (A.6) and (A.8) yields # #2 2 2 ρˆe,L ρˆe,H 1 ρˆe,L ∂

2 σρˆe ≈ − BL + BH σL2 4 uL uL uH ∂uL # #2 2 2 ρˆe,L ρˆe,H 1 ρˆe,H ∂

+ − BL + BH σH2 . 4 uH uL uH ∂uH

Defining 2 ρˆe,L

A≡

uL

BL +

2 ρˆe,H

uH

(A.8)

(A.9)

# BH ,

(A.10)

one can write

  1 ρˆe,L ∂ 2 2 1 ρˆe,H ∂ 2 2 −A σL + −A σH 4 uL ∂uL 4 uH ∂uH and therefore equation (A.11) becomes ⎧   2 2 ⎪ 1 ρˆe,H ρˆe,L 2AuH 2AuL ⎪ 2 ⎪ ⎪ − σL + + σH2 ⎨ 2 uL (uL + uH )2 uH (uL + uH )2  σρˆe ≈  ⎪ ⎪ ρˆe,L 1 A 2 2 ρˆe,H AuL 2 2 ⎪ ⎪ ⎩ − σL + + 2 σH 2 uL uH uH uH σρ2ˆe ≈

(A.11)

for = DEI . for = DER (A.12)

Since ρˆe and Zˆ are statistically correlated, their covariance is given by: ˆ σρˆe ,Zˆ ≡ COVAR(ρˆe , Z) ! " ∂ ρˆe ∂ Zˆ ∂ ρˆe ∂ Zˆ ∂ ρˆe ∂ Zˆ 2 ∂ ρˆe ∂ Zˆ 2 σH + + σ COVAR(uL , uH ) + ≈ ∂uH ∂uH ∂uH ∂uL ∂uL ∂uH ∂uL ∂uL L ≈

∂ ρˆe ∂ Zˆ 2 ∂ ρˆe ∂ Zˆ 2 σ + σ . ∂uH ∂uH H ∂uL ∂uL L 2084

(A.13)

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Note here that uL and uH are assumed statistically independent and therefore COVAR(uL , uH ) = 0. For = DEI, equation (A.13) becomes # K  ρˆe,H −uL 2AuL k−2 σ2 (k − 1) c ˆ

+ σρˆe ,Zˆ ≈ k (uL + uH )2 k=1 uH (uL + uH )2 H # K  ρˆe,L uH 2AuH k−2 σ 2, + (k − 1)cˆk

− (A.14) (uL + uH )2 k=1 uL (uL + uH )2 L and for = DER, equation (A.13) becomes # K ρˆe,H uL  AuL k−2 (k − 1)cˆk

+ 2 σH2 σρˆe ,Zˆ ≈ − 2 uH 2uH k=1 uH # K ρˆe,L 1  A k−2 σ 2. + (k − 1)cˆk

− 2uH k=1 uL uH L

(A.15)

Appendix B. Derivation of ion stopping power uncertainties The variance on the SP is given by (here, we take Zˆ = Zˆ med to ease the notation) ˆ σS2ˆ ≡ VAR(S) #2 # #2 # ˆ ˆ ˆ ∂ S ∂ S ∂ Sˆ ∂ S σρ2ˆe + σZ2ˆ + 2 ≈ σ ˆ ∂ρe ∂Z ∂ρe ∂Z ρˆe ,Z #2   2  k0 z2 2me c2 β 2 ρe k0 z2 ∂ Iˆ 2 2 ≈ σρˆe + − σZ2ˆ ln −β β2 I(1 − β 2 ) I β 2 ∂Z #   2  2me c2 β 2 k0 z ρe k0 z2 ∂ Iˆ 2 ln −β +2 − σ ˆ β2 I(1 − β 2 ) I β 2 ∂Z ρˆe ,Z #2 # #2 # Sˆ Sˆ ρe k0 z2 ∂ Iˆ ρe k0 z2 ∂ Iˆ 2 2 ≈ σρˆe + − σZˆ − 2 σρˆe ,Zˆ . (B.1) ρe ρe Iˆ β 2 ∂Z Iˆ β 2 ∂Z 

Appendix C. Derivation of the uncertainty on ion range Assuming the ion to cross N voxels of dimensions xi (in the direction of the beam), the average energy loss by the particle is approximately equal to E ≈

N 

xi Si ,

(C.1)

i=1

with Si the SP in the ith voxel. For equally-size voxels, xi ≡ x and therefore # N 1 Si . E ≈ Nx N i=1 2085

(C.2)

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Let us choose the initial kinetic energy of the ion E0 such that E = E0 . Hence, the estimated range is given by Rˆ ≈ Nx and therefore #  1 E0 ≈ Rˆ Si , (C.3) N i hence E0 Rˆ ≈ ' 1  N

i

Si

(.

(C.4)

Therefore the variance on the estimated range is given by ˆ σR2ˆ ≡ VAR(R) N E2 1  ≈ ' 0 (4 2 VAR(Si ) 1 N i S i i N

≈

Rˆ 2 σS2 2 S N

ˆ ≈ Rx

σS2 S

2

.

(C.5)

References Ahnesj¨o A 1989 Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media Med. Phys. 16 577–92 Alvarez R E and Macovski A 1976 Energy-selective reconstructions in x-ray computerized tomography Phys. Med. Biol. 21 733–44 Bazalova M, Carrier J-F, Beaulieu L and Verhaegen F 2008a Tissue segmentation in Monte Carlo treatment planning: a simulation study using dual-energy CT images Radiother. Oncol. 86 93–8 Bazalova M, Carrier J-F, Beaulieu L and Verhaegen F 2008b Dual-energy CT-based material extraction for tissue segmentation in Monte Carlo dose calculations Phys. Med. Biol. 53 2439–56 Bedwani S, Carrier J-F and Bouchard H 2013 SUET501: a sensitivity study of tissue characterization for brachytherapy Monte Carlo dose calculation Med. Phys. 40 320 Berger M J, Hubbell J H, Seltzer S M, Chang J, Coursey J S, Sukumar R and Zucker D S 1998 XCOM: photon cross sections database NIST Stand. Ref. Database 8 87–3597 Bischel H 1972 Passage of Charged Particles Through Matter (New York: McGraw-Hill) pp 8–142 Bouchard H, Seuntjens J, Carrier J-F and Kawrakow I 2009 Ionization chamber gradient effects in nonstandard beam configurations Med. Phys. 36 4654 Chandra N and Langan D A 2011 Gemstone detector: dual energy imaging via fast kVp switching Dual Energy CT in Clinical Practice (Berlin: Springer) pp 35–41 Chetty I J et al 2007 Report of the AAPM task group no 105: issues associated with clinical implementation of Monte Carlo-based photon and electron external beam treatment planning Med. Phys. 34 4818–53 Du Plessis F C P, Willemse C A, L¨otter M G and Goedhals L 1998 The indirect use of CT numbers to establish material properties needed for Monte Carlo calculation of dose distributions in patients Med. Phys. 25 1195–201 Flohr T G et al 2006 First performance evaluation of a dual-source CT (DSCT) system Eur. Radiol. 16 256–68 Guan H, Yin F-F and Kim J H 2002 Accuracy of inhomogeneity correction in photon radiotherapy from CT scans with different settings Phys. Med. Biol. 47 N223–31 Hounsfield G N 1973 Computerized transverse axial scanning (tomography): Part 1. Description of system Br. J. Radiol. 46 1016–22 2086

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

ICRU 1984 Stopping powers for electrons and positrons ICRU report 37 (Bethesda, MD: International Commission on Radiation Units and Measurements) ICRU 1989 Tissue substitutes in radiation dosimetry and measurement ICRU Report 44 (Bethesda, MD: International Commission on Radiation Units and Measurements) ICRU 1992 Photon, electron, proton and neutron interaction data for body tissues ICRU Report 46 (Bethesda, MD: International Commission on Radiation Units and Measurements) Jackson D F and Hawkes D J 1981 X-ray attenuation coefficients of elements and mixtures Phys. Rep. 70 169–233 Jiang H and Paganetti H 2004 Adaptation of GEANT4 to Monte Carlo dose calculations based on CT data Med. Phys. 31 2811–8 Johnson T 2011 Dual Energy CT in Clinical Practice (Berlin: Springer) Kirby B J, Davis J R, Grant J A and Morgan M J 2003 Extracting material parameters from x-ray attenuation: a CT feasibility study using kilovoltage synchrotron x-rays incident upon low atomic number absorbers Phys. Med. Biol. 48 3389–409 Landry G, Reniers B, Granton P V, van Rooijen B, Beaulieu L, Wildberger J E and Verhaegen F 2011 Extracting atomic numbers and electron densities from a dual source dual energy CT scanner: experiments and a simulation model Radiother. Oncol. 100 375–9 Mahnken A H, Stanzel S and Heismann B 2009 Spectral ρZ-projection method for characterization of body fluids in computed tomography: ex vivo experiments 1 Acad. Radiol. 16 763–9 Mayneord W V 1937 The significance of the roentgen Acta Union Int. Against Cancer 2 271282 Midgley S M 2004 A parameterization scheme for the x-ray linear attenuation coefficient and energy absorption coefficient Phys. Med. Biol. 49 307–25 Midgley S M 2011 Feasibility study for a novel method of dual energy x-ray analysis Phys. Med. Biol. 56 5599–619 Paganetti H 2012 Range uncertainties in proton therapy and the role of Monte Carlo simulations Phys. Med. Biol. 57 R99–R117 Paganetti H, Jiang H, Parodi K, Slopsema R and Engelsman M 2008 Clinical implementation of full Monte Carlo dose calculation in proton beam therapy Phys. Med. Biol. 53 4825–53 Petersilka M, Bruder H, Krauss B, Stierstorfer K and Flohr T G 2008 Technical principles of dual source CT Eur. J. Radiol. 68 362–8 Poludniowski G, Landry G, DeBlois F, Evans P M and Verhaegen F 2009 SpekCalc: a program to calculate photon spectra from tungsten anode x-ray tubes Phys. Med. Biol. 54 N433 Press W H, Teukolsky S A, Vetterling W T and Flannery B P 2002 Numerical Recipes in C++: the Art of Scientific Computing (Cambridge: Cambridge University Press) Rutherford R A, Pullan B R and Isherwood I 1976 Measurement of effective atomic number and electron density using an EMI scanner Neuroradiology 11 15–21 Saito M 2012 Potential of dual-energy subtraction for converting CT numbers to electron density based on a single linear relationship Med. Phys. 39 2021–30 Schaffner B and Pedroni E 1998 The precision of proton range calculations in proton radiotherapy treatment planning: experimental verification of the relation between CT-HU and proton stopping power Phys. Med. Biol. 43 1579–92 Schlomka J P et al 2008 Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography Phys. Med. Biol. 53 4031–47 Schneider U, Pedroni E and Lomax A 1996 The calibration of CT Hounsfield units for radiotherapy treatment planning Phys. Med. Biol. 41 111–24 Schneider W, Bortfeld T and Schlegel W 2000 Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions Phys. Med. Biol. 45 459–78 Snyder W S, Cook M J, Nasset E S, Karhausen L R, Howells G P and Tipton I H 1975 Report of the Task Group on Reference Man (ICRP Publication 23) (Oxford: Pergamon) Spiers W 1946 Effective atomic number and energy absorption in tissues Br. J. Radiol. 19 52–63 Thomas S 1999 Relative electron density calibration of CT scanners for radiotherapy treatment planning Br. J. Radiol. 72 781–6 Torikoshi M, Tsunoo T, Sasaki M, Endo M, Noda Y, Ohno Y, Kohno T, Hyodo K, Uesugi K and Yagi N 2003 Electron density measurement with dual-energy x-ray CT using synchrotron radiation Phys. Med. Biol. 48 673–85 Tsunoo T, Torikoshi M, Ohno Y, Uesugi K and Yagi N 2008 Measurement of electron density in dualenergy x-ray CT with monochromatic x rays and evaluation of its accuracy Med. Phys. 35 4924–32 2087

A E Bourque et al

Phys. Med. Biol. 59 (2014) 2059

Valentin J 2002 Basic anatomical and physiological data for use in radiological protection: reference values ICRP Publ. 89 Ann. ICRP 32 1–277 Vanderstraeten B et al 2007 Conversion of CT numbers into tissue parameters for Monte Carlo dose calculations: a multi-centre study Phys. Med. Biol. 52 539–62 Verhaegen F and Devic S 2005 Sensitivity study for CT image use in Monte Carlo treatment planning Phys. Med. Biol. 50 937–46 Williamson J F, Li S, Devic S, Whiting B R and Lerma F A 2006 On two-parameter models of photon cross sections: Application to dual-energy CT imaging Med. Phys. 33 4115–29 Yang M, Virshup G, Clayton J, Zhu X R, Mohan R and Dong L 2010 Theoretical variance analysis of single-and dual-energy computed tomography methods for calculating proton stopping power ratios of biological tissues Phys. Med. Biol. 55 1343–62

2088