Comparison of projection- and image-based methods for proton stopping power estimation using dual energy CT

Alvarez R.E.
Macovski A.

Energy-selective reconstructions in X-ray computerized tomography.

] and Rutherford et al. [

[4]

Rutherford R.A.
Pullan B.R.
Isherwood I.

Measurement of effective atomic number and electron density using an EMI scanner.

] observed that the linear attenuation coefficient,

μ

, could be separated into terms describing the two main interaction processes for X-ray photons in the energy range used in diagnostic CT: Photoelectric absorption and Compton scattering. Material decomposition can therefore be performed when having attenuation measurements at two different energies [

[6]

Heismann B.J.
Schmidt B.
Flohr T.

Spectral computed tomography.

]. The decomposition schemes proposed by Alvarez and Macovski [

Alvarez R.E.
Macovski A.

Energy-selective reconstructions in X-ray computerized tomography.

] and Rutherford et al. [

[4]

Rutherford R.A.
Pullan B.R.
Isherwood I.

Measurement of effective atomic number and electron density using an EMI scanner.

] differed both in the suggested basis functions and in the sense that Alvarez and Macovski based their estimation scheme on the projection data directly obtained from the CT scan and then reconstructed these basis images, whereas Rutherford et al. relied on reconstructed CT data of the linear attenuation for their estimations. These two approaches can be categorized as projection-based and image-based methods, respectively. The main difference in the approaches is therefore the order of the decomposition and the reconstruction step. In projection-based methods, material decomposition is performed prior to image reconstruction, and for image-based methods the decomposition takes place after image reconstruction.

Both projection-based and image-based approaches can be further divided into two main types: Type (1) Model-based: prior knowledge of source energy spectra and detector response is required to solve analytically the system of equations; Type (2) Calibration-based: X-ray measurements of materials with known properties are performed to obtain parameters that characterize the energy spectra and the detector response.

Several algorithms based on DECT have been proposed to predict the SPR or the intermediate variables required to compute the SPR, e.g., the relative electron density (RED), the effective atomic number

(Z_{eff})

or the logarithm of the mean excitation energy [

1

Yang M.
Virshup G.
Clayton J.
Zhu X.R.
Mohan R.
Dong L.

Theoretical variance analysis of single- and dual-energy computed tomography methods for calculating proton stopping power ratios of biological tissues.

,

7

Saito M.

Potential of dual-energy subtraction for converting CT numbers to electron density based on a single linear relationship.

,

8

Hünemohr N.
Krauss B.
Tremmel C.
Ackermann B.
Jäkel O.
Greilich S.

Experimental verification of ion stopping power prediction from dual energy CT data in tissue surrogates.

,

9

Hünemohr N.
Paganetti H.
Greilich S.
Jäkel O.
Seco J.

Tissue decomposition from dual energy CT data for MC based dose calculation in particle therapy.

,

10

Landry G.
Seco J.
Gaudreault M.
Verhaegen F.

Deriving effective atomic numbers from DECT based on a parameterization of the ratio of high and low linear attenuation coefficients.

,

11

Bourque A.E.
Carrier J.F.
Bouchard H.

A stoichiometric calibration method for dual energy computed tomography.

,

12

Farace P.

Comment: experimental verification of ion stopping power prediction from dual energy CT data in tissue surrogates.

,

13

van Abbema J.K.
van Goethem M.J.
Greuter M.J.W.
van der Schaaf A.
Brandenburg S.
van der Graaf E.R.

Relative electron density determination using a physics based parameterization of photon interactions in medical DECT.

,

14

Garcia L.I.R.
Azorin J.F.P.
Almansa J.F.

A new method to measure electron density and effective atomic number using dual-energy CT images.

,

15

Han D.
Siebers J.V.
Williamson J.F.

A linear, separable two-parameter model for dual energy CT imaging of proton stopping power computation.

,

16

Möhler C.
Wohlfahrt P.
Richter C.
Greilich S.

Range prediction for tissue mixtures based on dual-energy CT.

,

17

Taasti V.T.
Petersen J.B.B.
Muren L.P.
Thygesen J.
Hansen D.C.

A robust empirical parametrization of proton stopping power using dual energy CT.

]. At present, no consensus is reached within the community to prefer one method over the others, nor a decomposition domain (the projection- or the image-domain) for the SPR estimation. Nearly all DECT methods are image-based [

[18]

van Elmpt W.
Landry G.
Das M.
Verhaegen F.

Dual energy CT in radiotherapy: current applications and future outlook.

], possibly because image-based methods have the clear advantage that they can be performed on all DECT scanners, including conventional SECT scanners by acquiring two consecutive images at different kVp settings. In contrast, for projection-based methods the low and the high energy projections need to be acquired at the same angle or be accurately interpolated. Nevertheless, with the advent of dual-layer systems and photon-counting detector technology, projection-based methods are becoming more attractive since the projections from the different energy spectra are obtained without angular separation and there cannot be any motion between the two (or more) projections per source position. Another reason is that projection-based methods require access to dual-energy raw-data, which are not accessible in most commercial CT systems without the vendor agreement or reverse engineering.

It has been shown that DECT has a potential to improve the accuracy of the SPR estimation compared to conventional single energy CT (SECT) [

1

Yang M.
Virshup G.
Clayton J.
Zhu X.R.
Mohan R.
Dong L.

Theoretical variance analysis of single- and dual-energy computed tomography methods for calculating proton stopping power ratios of biological tissues.

,

19

Hansen D.C.
Seco J.
Srensen T.S.
Petersen J.B.B.
Wildberger J.E.
Verhaegen F.
et al.

A simulation study on proton computed tomography (CT) stopping power accuracy using dual energy CT scans as benchmark.

,

20

Hünemohr N.
Krauss B.
Dinkel J.
Gillmann C.
Ackermann B.
Jäkel O.
et al.

Ion range estimation by using dual energy computed tomography.

,

21

Zhu J.
Penfold S.N.

Dosimetric comparison of stopping power calibration with dual-energy CT and single-energy CT in proton therapy treatment planning.

], however, these studies were based on image-based methods. Tremblay et al. [

Tremblay J.-É.
Bedwani S.
Bouchard H.

A theoretical comparison of tissue parameter extraction methods for dual energy computed tomography.

] concluded that combining the output of projection-based methods with the image-based stoichiometric calibration proposed by Bourque et al. [

[11]

Bourque A.E.
Carrier J.F.
Bouchard H.

A stoichiometric calibration method for dual energy computed tomography.

] gave better accuracy for RED and

Z_{eff}

than a projection-based method alone. The aim of the present study was therefore to investigate the impact on the SPR accuracy when using a projection-based method. We compared the performances of the SPR estimation and the corresponding proton range errors of one projection- and three image-based DECT methods, through a simulation-based CT acquisition framework with the measured detector response and dual-energy spectra of a commercial DECT scanner and realistic noise levels.

2. Methods and materials

2.1 SPR estimation methods

Four DECT-based SPR estimation methods were investigated, including one projection-based and three image-based approaches. The projection-based method relied on the two-material decomposition proposed by Alvarez and Macovski [

Alvarez R.E.
Macovski A.

Energy-selective reconstructions in X-ray computerized tomography.

] to estimate the RED and the conversion from RED to SPR proposed by Kanematsu et al. [

Kanematsu N.
Inaniwa T.
Koba Y.

Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams.

]. We will refer to this method as “AMK” for Alvarez and Macovski, and Kanematsu. The first image-based method used a parametrization for the SPR proposed by Taasti et al. [

Taasti V.T.
Petersen J.B.B.
Muren L.P.
Thygesen J.
Hansen D.C.

A robust empirical parametrization of proton stopping power using dual energy CT.

]. This image-based method will be referred to as “SPP” for Stopping Power Parametrization. The second image-based approach was a combination of the methods proposed by Saito [

[7]

Saito M.

Potential of dual-energy subtraction for converting CT numbers to electron density based on a single linear relationship.

] to compute the RED and by Kanematsu et al. [

Kanematsu N.
Inaniwa T.
Koba Y.

Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams.

] to derive the SPR. We will refer to this method as “SK” for Saito and Kanematsu. The third image-based method is a two-material decomposition method with different basis materials for soft and bone tissues, suggested by Han et al. [

[15]

Han D.
Siebers J.V.
Williamson J.F.

A linear, separable two-parameter model for dual energy CT imaging of proton stopping power computation.

]. The four SPR estimation methods are described in Appendix A.

2.2 Virtual patient

The Adult Female (AF) reference computational phantom of the International Commission on Radiological Protection (ICRP) [

[24]

ICRP. Adult Reference Computational Phantoms. ICRP Publication 110. Ann. ICRP 39 (2).

Google Scholar

] was used as a virtual patient to evaluate the performance of the SPR methods. This phantom, with voxel dimensions of 1.775 × 1.775 × 4.84 mm³, represented an average female subject structured with 140 organs made of 53 standard human tissues. The mass density and elemental weights of each tissue were provided in the ICRP Publication 110 [

[24]

ICRP. Adult Reference Computational Phantoms. ICRP Publication 110. Ann. ICRP 39 (2).

Google Scholar

]. To evaluate a large number of tissue densities and compositions, four anatomical regions were selected: head, sternum, breast and pelvis (Fig. 1).

Figure thumbnail gr1 — Fig. 1Reference SPR for the four slices used in the comparison. Placement of ROIs are shown by blue circles, three to four ROIs are placed in each slice. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

2.3 Calculation of reference SPR

Reference SPR values of the ICRP phantom were computed from the chemical composition and density pixel-by-pixel using the Bethe equation without correction terms, as described by Schneider et al. [

[25]

Schneider U.
Pedroni E.
Lomax A.

The calibration of CT Hounsfield units for radiotherapy treatment planning.

]. A proton beam of 200 MeV kinetic energy was considered and the mean ionization energy of water was set to 78 eV [

[26]

Sigmund P, Schinner A, Paul H. Errata and Addenda for ICRU Report 73. Stopping of ions heavier than helium. J ICRU 2009;5(1).

Google Scholar

]. The mean ionization energies of the tissues were calculated using the Bragg additivity rule [

[25]

Schneider U.
Pedroni E.
Lomax A.

The calibration of CT Hounsfield units for radiotherapy treatment planning.

] and the I-values given in Tables 2.8 and 2.11 in the ICRU49 report [

[27]

ICRU. Stopping Power and Ranges for Protons and Alpha Particles. ICRU Report No. 49.

Google Scholar

].

2.4 CT imaging and CT reconstruction

Virtual CT projections of the Gammex and the ICRP phantom were generated in GATE [

[28]

Jan S.
Santin G.
Strul D.
Staelens S.
Assié K.
Autret D.
et al.

GATE: a simulation toolkit for PET and SPECT.

] v7.2 (based on Geant4 v10.1, physics-list: emlivermore) using the Fixed Forced Detection Actor. This module deterministically computes digitally reconstructed radiographs using the reconstruction toolkit (RTK) v1.3.0 [

[29]

Rit S.
Vila Oliva M.
Brousmiche S.
Labarbe R.
Sarrut D.
Sharp G.C.

The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK).

] and the Geant4 database of X-rays cross sections. Fan-beam projections with 600 views of 2052 pixels subsequently re-binned to 1026 pixels of 1 mm² were acquired. The thickness of the detector row was 1 mm for the Gammex acquisitions and 2 mm for the ICRP phantom. The source-to-isocenter distance was 626 mm and the source-to-detector distance was 1026 mm. For the simulations of the CT projections, the measured dual-energy spectra (low energy (LE): 100 kVp, high energy (HE): 150 kVp + 0.6 mm Sn) and the measured detector response for the SOMATOM Force dual-source CT scanner were used, kindly provided by Siemens Healthcare (Forchheim, Germany). To have a fair comparison between the projection-based and the image-based methods, energy spectra with high mean energies were chosen to reduce the beam-hardening for the image-based methods. To represent a realistic scenario, Poisson distributed noise was applied to the projections. For each slice, a total central dose of 20 mGy was delivered, with an equal dose split between the two energy spectra, thus a central dose of 10 mGy per energy spectrum; the calculation of the number of photons needed to deliver a given dose is described in Appendix C. No bowtie-filter was simulated.

Prior to CT reconstruction, the LE and the HE sinograms of the Gammex and the AF phantom were corrected for beam-hardening using the method described in Ref. [

[30]

Brooks R.A.
Di Chiro G.

Beam hardening in X-ray reconstructive tomography.

], since this was found to be mandatory for the image-based methods. A look-up table linearizing the input projection values for water [

[30]

Brooks R.A.
Di Chiro G.

Beam hardening in X-ray reconstructive tomography.

] was created using the energy spectra and the detector response and increasing water thicknesses.

CT images of the Gammex phantom and the AF slices were reconstructed using the LE and the HE sinograms, for the image-based methods; and soft tissue and cortical bone density images of the AF slices, for the projection-based method. All reconstructions were performed using the filtered-backprojection (FBP) reconstruction of RTK; for the reconstructions of the AF slices preserving the initial pixel dimension of 1.775 × 1.775 mm² and for the Gammex the pixel dimension of 1 × 1 mm².

2.5 SPR comparison

The performance of the SPR estimation was evaluated based on different criteria: relative SPR differences over defined ROIs, each covering a single tissue, and relative range errors computed over a whole slice. The placement of the ROIs can be seen in Fig. 1, and their reference SPRs are listed in Table B.5 in Appendix B. For each ROI, the mean error was calculated and its precision was assessed using the standard error of the mean (SEM). To have a direct comparison of the three methods, the root-mean-square error (RMSE) over the mean relative SPR difference for the thirteen ROIs was taken; this gave a combined error measure for the thirteen defined ROIs.

To estimate the influence of the SPR deviations on the range calculation, we calculated the range errors along the proton beam path through the entire slice using the Radon transform implemented in MATLAB (The MathWorks Inc., Natick, MA). The Radon transform was computed on the absolute SPR difference images, for angles in the interval from 0–179° in steps of 1°. The SPR difference for pixels outside the body outline was set to zero. The Radon transform was also performed on the reference SPR maps to calculate the water-equivalent path length (WEPL). To exclude beam paths entirely outside the body, an exclusion penalty was applied to pixels in the reference SPR sinogram with a WEPL of zero. The results for the range errors were given relative to reference WEPL maps, to take into account that the beam paths in different directions and slices were not of equal length.

3. Results

The RMSE over the mean relative error for each of the thirteen ROIs was 0.54% for the AMK method, 0.68% for the SPP method, 0.61% for the SK method, and 0.70% for Han’s method (Fig. 2). The bias (given by the mean error) for the head, sternum and breast slices was smallest for the SPP method, however, for the pelvis slice (ROIs 4a-d) this method gave the largest errors (Table 1). Considering all slices, the SK method had the largest bias. For all methods, the relative errors for ROIs in the pelvis slice (ROIs 4a-d) were larger than for the other slices.

Figure thumbnail gr2 — Fig. 2Relative SPR errors for each of the thirteen ROIs (for the placement see Fig. 1). The center squares show the mean of the relative SPR error over the ROI, and the error-bars show the standard error of the mean (SEM).
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

Table 1Comparison of relative SPR deviations over the defined ROIs. The upper part of the table gives the results averaged over the individual slices, while the lower part gives the results averaged over all four slices together. (PB: Projection-based, IB: Image-based).

Per slice		Head	Sternum	Breast	Pelvis
PB (AMK)	Mean	0.28%	0.32%	0.29%	−0.45%
PB (AMK)	Uns. mean	0.28%	0.32%	0.64%	0.66%
IB (SPP)	Mean	−0.06%	0.03%	0.17%	−0.99%
IB (SPP)	Uns. mean	0.46%	0.30%	0.29%	0.99%
IB (SK)	Mean	−0.17%	−0.06%	−0.09%	−0.85%
IB (SK)	Uns. mean	0.20%	0.20%	0.40%	0.85%
IB (Han)	Mean	0.01%	0.26%	0.79%	−0.60%
IB (Han)	Uns. mean	0.34%	0.26%	0.79%	0.86%

All slices		RMSE	Mean	Uns. mean	$σ$
PB (AMK)		0.54%	0.07%	0.49%	0.56%
IB (SPP)		0.68%	−0.27%	0.55%	0.65%
IB (SK)		0.61%	−0.33%	0.44%	0.53%
IB (Han)		0.70%	0.06%	0.59%	0.73%

Open table in a new tab

The AMK method gave the lowest RMSEs for the relative range errors. For the AMK and SPP methods, the positive and negative range errors were nearly averaged out over the slice, giving mean relative range errors close to 0%; only for the head slice estimated by the SPP method the 0% range error was not within the interquartile range (IQ) (Fig. 3). In contrast, for the image-based SK method, the range error distribution had a negative bias for all slices, and Han’s method gave a slight positive bias for some slices. Even though the mean relative range errors for the SPP method were fairly low, this method produced the widest error distributions (Fig. 3 and Table 2).

Figure thumbnail gr3 — Fig. 3Box-plot of relative range error distribution, as calculated by the Radon transform for each anatomical slice using each of the four SPR estimation methods (PB: Projection-based, IB: Image-based). The boxes show the interquartile (IQ) range, i.e. from the 25%-percentile (Q1) to the 75%-percentile (Q3), and the black dots represent the median. The whiskers go from $Q 1 - 2.5 \times IQ$ up to $Q 3 + 2.5 \times IQ$ . The outliers are not shown in the figure, but the percentages (excluding the sinogram entries corresponding to a WEPL of zero) are given in the right side of the figure.
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

Table 2Comparison of the relative range errors over each image slice. The results are given as the signed mean errors

(μ \pm σ)

and the root-mean-square errors (RMSE) in percentage. The range errors are taken relative to the reference water equivalent path length (WEPL); the mean reference WEPL for each slice is given in parentheses beneath the slice name. (PB: Projection-based, IB: Image-based).

	Head (122.8 mm)		Sternum (162.5 mm)		Breast (162.1 mm)		Pelvis (181.7 mm)
Method	$μ \pm σ$ (%)	RMSE (%)	$μ \pm σ$ (%)	RMSE (%)	$μ \pm σ$ (%)	RMSE (%)	$μ \pm σ$ (%)	RMSE (%)
PB (AMK)	−0.24 ± 0.81	0.84	−0.01 ± 0.64	0.64	−0.04 ± 0.58	0.58	−0.14 ± 0.80	0.82
IB (SPP)	−0.41 ± 0.86	0.95	0.01 ± 0.81	0.81	0.04 ± 0.72	0.72	0.03 ± 1.10	1.10
IB (SK)	−0.51 ± 0.80	0.95	−0.28 ± 0.63	0.69	−0.33 ± 0.57	0.66	−0.44 ± 0.80	0.91
IB (Han)	−0.10 ± 0.82	0.82	0.36 ± 0.75	0.83	0.35 ± 0.65	0.74	0.47 ± 1.19	1.28

Open table in a new tab

4. Discussion

In this study, we have compared one projection-based and three image-based SPR methods. We found that the projection-based method produced slightly better results than the three image-based methods used in this study. The differences were not found to be statistically significant. The image-based SPP method gave mean errors comparable to the results for the projection-based method, but larger standard deviations. In contrast, the image-based SK method resulted in slightly biased results. It has been shown earlier that the SPR accuracy using Saito’s method in combination with an estimation of

Z_{eff}

is dependent on the calibration and evaluation materials being very similar [

[19]

Hansen D.C.
Seco J.
Srensen T.S.
Petersen J.B.B.
Wildberger J.E.
Verhaegen F.
et al.

A simulation study on proton computed tomography (CT) stopping power accuracy using dual energy CT scans as benchmark.

].

In our simulations, clinical polychromatic energy spectra, a realistic detector response and noise levels corresponding to a 20 mGy central dose were used. The material basis for the projection-based approach was optimized to give the best results. To have a direct comparison between projection-based and image-based calculations, we used an image-based method (the SK method) where the conversion from RED to SPR were performed in the same way as for the projection-based method, using the polyline fit proposed by Kanematsu et al. [

Kanematsu N.
Inaniwa T.
Koba Y.

Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams.

], and we also compared to a two-material decomposition image-based method, though not using the exact same materials as for the projection-based method, as we chose to use the optimal basis materials in each case. Further, we also used one of the latest image-based methods [

Taasti V.T.
Petersen J.B.B.
Muren L.P.
Thygesen J.
Hansen D.C.

A robust empirical parametrization of proton stopping power using dual energy CT.

], as this method was proven to give smaller SPR errors and to be more robust to noise than two other existing DECT methods [

15

Han D.
Siebers J.V.
Williamson J.F.

A linear, separable two-parameter model for dual energy CT imaging of proton stopping power computation.

,

19

Hansen D.C.
Seco J.
Srensen T.S.
Petersen J.B.B.
Wildberger J.E.
Verhaegen F.
et al.

A simulation study on proton computed tomography (CT) stopping power accuracy using dual energy CT scans as benchmark.

]. For these reasons, we consider that the comparison was carried out under fair conditions.

The image-based SK method showed a net negative bias on the range accuracy that was not observed for the image-based SPP method. Range errors seem to correlate well with SPR errors obtained for some tissues, such as breast tissue (ROI 3c), adipose tissue (ROI 4d) or pelvis spongiosa (ROI 4c), which were found in large proportions in the selected slices. This suggests that accuracy errors add up along the proton beam path. The SPP method resulted in noisier SPR images and, thus, wider error distributions. In terms of RMSE, the projection-based approach produced better results for the SPR and the range with respect to the three image-based methods. Therefore, based on these results, we recommend increasing the consideration of projection-based approaches for proton SPR estimation. However, the inherent limitations of each method should be taken into account, and therefore the most well-suited SPR estimation method for a specific facility might depend on the DECT technology available at the facility. It can therefore not be generally stated which method will produce the most accurate results.

Beam-hardening artifacts occur because the energy dependency of the linear attenuation coefficients is neglected. To overcome this approximation, image-based methods require a beam-hardening correction. In contrast, material decomposition in projection-based approaches is performed prior to reconstruction, and as the energy dependency is factored out in this decomposition process (see Eq. (A.1)) the basis-images should not be affected by beam-hardening. But this is not true unless a good choice of the basis-material decomposition is performed, such that the assumption that

μ

can be separated into an energy-dependent and an energy-independent part holds true. In the energy range of clinical CT scanners (i.e. from 40 to 150 keV) and for materials with

Z < 15

this assumption is acceptable [

[31]

Jackson D.F.
Hawkes D.J.

X-ray attenuation coefficients of elements and mixtures.

]. In this study, the implemented strategy to correct for beam-hardening in the image-based method works well for water-equivalent tissues but might not be accurate enough for low- and high-density tissues. A more sophisticated beam-hardening correction might therefore improve the results for the image-based methods. If a better beam-hardening correction had been applied, a LE spectrum with a lower mean energy might be more relevant than the 100 kVp spectrum used in this study.

The three image-based methods used in this study were calibration-based (type 2), and a calibration was needed to find the effective energies for the LE and HE spectrum. This calibration procedure is sensitive to the calibration phantom being comparable to the object for which the SPR is to be estimated; i.e. the calibration phantom should be of a similar size and consist of materials with a composition comparable to human tissues, such that the phantom hardens the X-ray energy spectrum in the same way as the investigated anatomical site. For this study, only a single calibration phantom with a diameter of 16.5 cm was used for all four anatomical slices. This may be one cause for the larger relative SPR deviations in the pelvis slice for the image-based methods (Table 1).

For model-based (type 1) image-based methods, the energy spectrum must also be properly hardened to reproduce the beam-hardening of the scanned object. This issue is taken into account by additionally filtering the X-ray source for instance by water [

1

Yang M.
Virshup G.
Clayton J.
Zhu X.R.
Mohan R.
Dong L.

Theoretical variance analysis of single- and dual-energy computed tomography methods for calculating proton stopping power ratios of biological tissues.

,

32

Bazalova M.
Carrier J.-F.
Beaulieu L.
Verhaegen F.

Dual-energy CT-based material extraction for tissue segmentation in Monte Carlo dose calculations.

]. Here, the amount of water may also be dependent on the anatomical region being investigated, since the performance of these methods will depend on the accuracy of this prior knowledge on the energy spectrum and the detector response [

[33]

Yang M.

Dual energy computed tomography for proton therapy treatment planning.

Google Scholar

]. Therefore, the input parameters must be tuned with care for both types of image-based methods. However, for calibration-based methods (type 2), this is more straight-forward since it only requires using a well-suited calibration phantom of the proper size.

For the projection-based method applied in this study, the system of Eqs. ((A.4), (A.5)) was solved numerically assuming that the energy spectra and the detector response were well-known (type 1). The same was assumed for the beam-hardening correction of image-based methods. It was shown that projection-based methods are sensitive to the amount of noise [

[34]

Brendel B, Bergner F, Brown K, Koehler T. Penalized likelihood decomposition for dual layer spectral CT. In: Proc 4th intl mtg on image formation in X-ray CT; 2016. p. 41–4.

Google Scholar

] and to alterations of the energy spectrum [

Tremblay J.-É.
Bedwani S.
Bouchard H.

A theoretical comparison of tissue parameter extraction methods for dual energy computed tomography.

]. Tremblay et al. found that for projection-based methods the accuracy of the RED was more degraded when the spectrum was altered than when using an image-based method [

Tremblay J.-É.
Bedwani S.
Bouchard H.

A theoretical comparison of tissue parameter extraction methods for dual energy computed tomography.

]. The projection-based results in the present work are therefore dependent on accurate knowledge of the energy spectrum and detector response. An alternative would be to approximate the system of Eqs. ((A.4), (A.5)) by a polynomial expansion of two variables [

3

Alvarez R.E.
Macovski A.

Energy-selective reconstructions in X-ray computerized tomography.

,

35

Chuang K.-S.
Huang H.K.

Comparison of four dual energy image decomposition methods.

]. The unknown coefficients could be determined experimentally through a calibration procedure with increasing thicknesses of two well-known materials (type 2). However, there would also be disadvantages with this approach, including sensitivity to the order of the polynomial, to the use of cross-terms in the polynomial, to the exact thicknesses of the interposed slabs of two materials and to the number of points of the calibration curve. It should also be noted that type 1 approaches always rely on a pre-calibration of the model, e.g., the one presented in [

[36]

Vilches-Freixas G.
Létang J.M.
Brousmiche S.
Romero E.
Oliva M.V.
Kellner D.
et al.

Technical note: procedure for the calibration and validation of kilo-voltage cone-beam CT models.

].

The basis functions used for the projection-based method were found to give the best SPR results, after testing several different basis functions. Other basis functions have been suggested in the literature, e.g. Han et al. [

[15]

Han D.
Siebers J.V.
Williamson J.F.

A linear, separable two-parameter model for dual energy CT imaging of proton stopping power computation.

] used two different sets of basis functions for soft and bone tissues, respectively, where the separation was based on the ratio of the CT numbers which is not available for projection-based methods. This method was tested in image-space in this study. Tremblay et al. [

Tremblay J.-É.
Bedwani S.
Bouchard H.

A theoretical comparison of tissue parameter extraction methods for dual energy computed tomography.

] found that in most cases a two-material decomposition rather than decomposition into physical processes performed better.

Another ongoing discussion within the community is how to determine the reference SPR values. In this study, a 200 MeV proton energy was used to determine the reference SPR, as for proton energies in the range between 80 and 300 MeV, the variation of SPR with proton energy is negligible

(< 1 %)

[

[37]

Arbor N.
Dauvergne D.
Dedes G.
Létang J.M.
Parodi K.
Quiñones C.T.
et al.

Monte Carlo comparison of X-ray and proton CT for range calculations of proton therapy beams.

]. In addition, Yang [

[33]

Yang M.

Dual energy computed tomography for proton therapy treatment planning.

Google Scholar

] conducted a study with a 175 MeV initial kinetic energy beam and concluded that neglecting the SPR dependence with energy introduced a 0.5% range error.

5. Conclusion

In this simulation-based study, the dual-energy decomposition domains – the projection-domain (decomposition prior to image reconstruction) and the image-domain (decomposition after image reconstruction) – were compared in terms of the performance of extracting the SPR from DECT and the resulting WEPL. In terms of SPR accuracy and range errors, the observed differences between projection- and image-based methods were not statistically significant. If properly calibrated, the two SPR estimation approaches can therefore be assumed to give similar results. With the advent of technological developments, projection-based methods are becoming more attractive. Most of the existing methods to compute the SPR are image-based but, based on these results, projection-based approaches should also be considered for proton SPR determination.

Disclosure of conflicts of interest

The authors have no relevant conflicts of interest to disclose.

Acknowledgments

This work was partially supported by grant ANR-13-IS03-0002-01 (DEXTER project) from the French National Research Agency (ANR), and by the Danish Cancer Society (Grant No. A5992).

Appendix A. SPR estimation methods

A.1 Projection-based method (AMK)

The method used to compute the SPR from dual-energy projection data is an adaptation of the two-step procedure proposed by Farace [

[12]

Farace P.

Comment: experimental verification of ion stopping power prediction from dual energy CT data in tissue surrogates.

] in the image-domain. This method was chosen for its simplicity and because it produced comparable results with respect to Hnemohr’s approach [

[8]

Hünemohr N.
Krauss B.
Tremmel C.
Ackermann B.
Jäkel O.
Greilich S.

Experimental verification of ion stopping power prediction from dual energy CT data in tissue surrogates.

]. First, RED was derived from the reconstructed mass density images of soft tissue (ST) and cortical bone (CB) [

[38]

P. Linstrom, W. Mallard, NIST Standard Reference Database Number 69, [Online] Available:http://webbook.nist.gov, National Institute of Standards and Technology, Gaithersburg, MD. (retrieved December 6, 2016).

Google Scholar

], issued from the implementation of the two-material decomposition method in the projection domain proposed by Alvarez and Macovski [

Alvarez R.E.
Macovski A.

Energy-selective reconstructions in X-ray computerized tomography.

]. Second, RED was converted into SPR using the poly-line relations proposed by Kanematsu et al. [

Kanematsu N.
Inaniwa T.
Koba Y.

Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams.

].

The basis material decomposition (BMD) method assumed that the linear attenuation coefficient of the scanned object at a location

x

can be expressed as a linear combination of two energy-dependent basis functions of two materials with energy-independent coefficients:

μ (E, x) = {(\frac{μ}{ρ})}_{ST} (E) \cdot a_{ST} (x) + {(\frac{μ}{ρ})}_{CB} (E) \cdot a_{CB} (x)

(A.1)

where

(μ / ρ)

represented the mass attenuation coefficient of ST and CB, and a is the mass fraction of ST and CB in the volume at location

x

– in this BMD these mass fractions represented the energy-independent coefficients.

The projection of the mass fractions

a_{ST}

and

a_{CB}

was expressed as follows:

A_{ST} (u, θ) = \int_{L (u, θ)} a_{ST} (ℓ) d ℓ

(A.2)

A_{CB} (u, θ) = \int_{L (u, θ)} a_{CB} (ℓ) d ℓ

(A.3)

where

ℓ \in L (u, θ)

was the line-segment between the source and a detector pixel located at position

u

for a given projection angle

θ

.

When performing a dual-energy acquisition, two sinograms of the same object at two different energy spectra were available:

I_{LE} (A_{ST}, A_{CB}) = \int_{E_{min}}^{E_{max}} S_{LE} (E) \exp (- A_{ST} {(\frac{μ}{ρ})}_{ST} (E) - A_{CB} {(\frac{μ}{ρ})}_{CB} (E)) d E

(A.4)

I_{HE} (A_{ST}, A_{CB}) = \int_{E_{min}}^{E_{max}} S_{HE} (E) \exp (- A_{ST} {(\frac{μ}{ρ})}_{ST} (E) - A_{CB} {(\frac{μ}{ρ})}_{CB} (E)) d E

(A.5)

where

I_{LE}

and

I_{HE}

were the measured intensities for the LE and the HE spectrum, respectively, for a given projection angle;

S_{LE}

and

S_{HE}

were the normalized energy spectra weighted by the detector response.

The unknowns

A_{ST}

and

A_{CB}

were determined by solving the system of Eqs. ((A.4), (A.5)) using the implementation of Nelder and Mead [

[39]

Nelder J.A.
Mead R.

A simplex method for function minimization.

] in RTK [

[29]

Rit S.
Vila Oliva M.
Brousmiche S.
Labarbe R.
Sarrut D.
Sharp G.C.

The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK).

], assuming that the energy spectra and the detector response were known (type 1 decomposition).

Then, the mass fraction per volume of soft tissue

(a_{ST})

and cortical bone

(a_{CB})

was determined by reconstruction of

A_{ST}

and

A_{CB}

(see Section 2.4). Afterwards, the RED was derived for each pixel by:

RED (x) = \frac{a_{ST} (x) {[\sum_{i} w_{i} \frac{Z_{i}}{A_{i}}]}_{ST} + a_{CB} (x) {[\sum_{i} w_{i} \frac{Z_{i}}{A_{i}}]}_{CB}}{ρ_{W} {[\sum_{i} w_{i} \frac{Z_{i}}{A_{i}}]}_{W}}

(A.6)

where

Z_{i}

was the atomic number,

A_{i}

was the atomic mass and

w_{i}

was the elemental weight fraction for element i of the tabulated compounds ST, CB and water (represented with the index W) [

[38]

P. Linstrom, W. Mallard, NIST Standard Reference Database Number 69, [Online] Available:http://webbook.nist.gov, National Institute of Standards and Technology, Gaithersburg, MD. (retrieved December 6, 2016).

Google Scholar

];

ρ_{W}

was the mass density for water.

To reproduce the piece-wise linear relations suggested by Kanematsu et al. [

Kanematsu N.
Inaniwa T.
Koba Y.

Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams.

], the SPR and RED were calculated for 92 ICRU46 tissues of known chemical composition and mass density [

[40]

ICRU, Photon, electron, proton and neutron interaction data for body tissues, ICRU Report No. 46.

Google Scholar

] applying the Bethe equation (see Section 2.3). The fitting lines are presented in Fig. A.1.

Figure thumbnail gr4 — Fig. A.1Kanematsu fitting lines to convert from RED to SPR/RED.
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

A.2 Image-based method – SPR parametrization (SPP)

The details of the image-based SPP method were described by Taasti et al. [

Taasti V.T.
Petersen J.B.B.
Muren L.P.
Thygesen J.
Hansen D.C.

A robust empirical parametrization of proton stopping power using dual energy CT.

], but the main principles are described here.

CT scans of a calibration phantom were required for the characterization of the X-ray energy spectra (type 2 decomposition). This calibration process was performed on simulated CT images of the Gammex 467 electron density calibration phantom (Gammex, Middleton, WI). The dimensions of the phantom were down-scaled to half the size of the real Gammex phantom which has a diameter of 33 cm; this was done to obtain the same beam hardening effect on this calibration phantom as for the ICRP image slices. Average CT numbers over ROIs inside the inserts in the reconstructed CT images were used.

Based on the known density and elemental composition of the inserts of the Gammex phantom, the effective energies of the used energy spectra were found using the CT numbers

(H_{j})

for the phantom inserts. The effective energy was defined as the energy which maximized the coefficient of determination,

R^{2}

, for the linear fit:

\frac{μ (E_{eff, j})}{μ_{W} (E_{eff, j})} = \frac{H_{j}}{A_{j}^{t}} + B_{j}^{t}

(A.7)

Here, subscript j refers to the energy spectrum

(j = LE,HE)

, and A and B are fitting parameters. The linear attenuation coefficients,

μ (E)

, for the Gammex inserts were calculated based on XCOM data [

[41]

Berger MJ, Hubbell JH, Seltzer SM, Chang J, Coursey JS, Sukumar R, et al., XCOM: photon cross section database (version 1.5), [Online] Available: http://physics.nist.gov/xcom, National Institute of Standards and Technology, Gaithersburg, MD; 2010.

Google Scholar

]. When the effective energies were found, the Gammex inserts were divided into two categories, soft and bone tissues based on their CT numbers in the LE CT image. Then Eq. (A.7) was refitted for each tissue group, to find two sets of fitting parameters,

A_{j}^{t}

and

B_{j}^{t}

, where superscript

t \in \{soft,bone\}

indicates the tissue group, using the effective energies

E_{eff, j}

for each energy spectrum j. In this study we used a separation point between the tissue groups set to

H_{LE} = 150

HU.

The attenuation ratios, u, for the 92 reference human tissues [

[40]

ICRU, Photon, electron, proton and neutron interaction data for body tissues, ICRU Report No. 46.

Google Scholar

] were calculated at the effective energies based on their density and elemental composition:

u \equiv \frac{μ (E_{eff})}{μ_{W} (E_{eff})} = ρ \frac{\sum_{i} w_{i} {(\frac{μ}{ρ})}_{i} (E_{eff})}{{[\sum_{i} w_{i} {(\frac{μ}{ρ})}_{i} (E_{eff})]}_{W}}

(A.8)

These attenuation ratios were fitted to the SPR for the tissues based on the following equations:

{SPR}_{soft}^{est} = (1 + x_{1}) u_{HE} - x_{1} u_{LE} + x_{2} u_{LE}^{2} + x_{3} u_{HE}^{2} + x_{4} (u_{LE}^{3} + u_{HE}^{3})

(A.9)

{SPR}_{bone}^{est} = (1 + x_{5}) u_{HE} - x_{5} u_{LE} + x_{6} \frac{u_{LE}}{u_{HE}} + x_{7} (u_{LE}^{2} - u_{HE}^{2}) + x_{8} (u_{LE}^{3} + u_{HE}^{3})

(A.10)

where the

x_{i}

’s are fitting parameters. The fitting parameters used in this study can be found in Table A.1. When these expressions were used to estimate the SPR for the AF phantom, the attenuation ratios were calculated using the fitting parameters found together with the effective energies,

u_{j}^{t} = H_{j} / A_{j}^{t} + B_{j}^{t}

. The same separation between soft and bone tissue was used for the SPR estimation as for the calculation of the attenuation ratios.

Table A.1Calibration parameters used in the SPR parametrization; for the energy spectra characterization (Eq. A.7)), and for the SPR estimation (Eq. (A.9), (A.10)). The effective energies of the low energy (LE) and high energy (HE) spectra are given in parentheses.

Energy spectra characterization			SPR fitting parameters
	LE (64 keV)	HE (96 keV)	Soft tissues		Bone tissues
$A^{soft}$	988.8	991.3	$x_{1}$	3.161	$x_{5}$	0.8251
$A^{bone}$	971.8	984.8	$x_{2}$	1.176	$x_{6}$	0.03853
$B^{soft}$	1.006	1.007	$x_{3}$	−1.136	$x_{7}$	0.1150
$B^{bone}$	0.9803	1.004	$x_{4}$	−0.01883	$x_{8}$	−0.008910

Open table in a new tab

A.3 Image-based method – Saito and Kanematsu’s (SK) method

This image-based DECT method is a combination of two methods originally presented by Saito [

[7]

Saito M.

Potential of dual-energy subtraction for converting CT numbers to electron density based on a single linear relationship.

] (calculation of the relative electron density, RED), and by Kanematsu et al. [

Kanematsu N.
Inaniwa T.
Koba Y.

Relationship between electron density and effective densities of body tissues for stopping, scattering, and nuclear interactions of proton and ion beams.

] (conversion from RED to SPR).

In Saito’s method the RED was estimated as

RED = a \frac{(1 + α) H_{HE} - α H_{LE}}{1000} + b

(A.11)

The constants a, b and

α

were found by making calibration fits to the theoretical

ρ_{e}

values for the 92 reference human tissues and their CT numbers calculated from Eqs. ((A.7), (A.8)). To take the low RED values for lung tissue into account the constants were found by minimizing the relative deviations

({RED}^{theo} - {RED}^{est}) / {RED}^{theo}

. The parameters used in this study can be seen in Table A.2. The CT numbers for the 92 reference human tissues used in the calibration were calculated from Eq. (A.8) and the effective energies given in Table A.1.

Table A.2Calibration parameters for Saito’s method, found by minimization of Eq. A.11. These parameters were used for calculation of RED in the SK method applied in this study.

a	b	$α$
1.0085	1.0091	0.5202

Open table in a new tab

The RED estimates were converted into SPR using the same calibration curve as described in Section A.1 for the Kanematsu method (see Fig. A.1).

A.4 Image-based method – Han’s (Han) method

The image-based two-material decomposition method used in this study was proposed by Han et al. [

[15]

Han D.
Siebers J.V.
Williamson J.F.

A linear, separable two-parameter model for dual energy CT imaging of proton stopping power computation.

], and was also explained in a compressed way in Appendix A.2 in Ref. [