Optimum size of a calibration phantom for x-ray CT to convert the Hounsfield units to stopping power ratios in charged particle therapy treatment planning.

In charged-particle therapy treatment planning, the volumetric distribution of stopping power ratios (SPRs) of body tissues relative to water is used for patient dose calculation. The distribution is conventionally obtained from computed tomography (CT) images of a patient using predetermined conversion functions from the CT numbers to the SPRs. One of the biggest uncertainty sources of patient SPR estimation is insufficient correction of beam hardening arising from the mismatch between the size of the patient cross section and the calibration phantom for producing the conversion functions. The uncertainty would be minimized by selecting a suitable size for the cylindrical water calibration phantom, referred to as an 'effective size' of the patient cross section, Leffective. We investigated the Leffective for pelvis, abdomen, thorax, and head and neck regions by simulating an ideal CT system using volumetric models of the reference male and female phantoms. The Leffective values were 23.3, 20.3, 22.7 and 18.8 cm for the pelvis, abdomen, thorax, and head and neck regions, respectively, and the Leffective for whole body was 21.0 cm. Using the conversion function for a 21.0-cm-diameter cylindrical water phantom, we could reduce the root mean square deviation of the SPRs and their mean deviation to ≤0.011 and ≤0.001, respectively, in the whole body. Accordingly, for simplicity, the effective size of 21.0 cm can be used for the whole body, irrespective of body-part regions for treatment planning in clinical practice.

INTRODUCTION

In charged-particle therapy treatment planning, the accurate prediction of particle range in patients is essential for conformal dose delivery to the target. Particle range is determined by integrating the stopping power ratios (SPRs) of body tissues relative to water along the beam path in a patient. The volumetric distribution of the SPRs in a patient is conventionally obtained from the x-ray computed tomography (CT) data, using a predetermined polyline relationship between CT number and SPR of the body tissues, referred to as the CT number-to-SPR conversion function [1-3]. Uncertainties of SPR estimation can induce range uncertainties of up to 3.5% [4, 5], which in current clinical practice are considered by adjusting the corresponding distal and proximal margins to the target. Yang et al. [6] grouped the uncertainties in the SPR estimation into several categories according to their sources and estimated their degrees for (lung tissues, soft tissues, bone tissues): uncertainties in patient CT imaging (3.3%, 0.6%, 1.5%), uncertainties related to the CT-number-to-SPR conversion functions (3.8%, 1.4%, 1.7%), uncertainties in mean excitation energies (0.2%, 0.2%, 0.6%) and uncertainties due to the energy dependence of SPR not commonly accounted for by a dose algorithm (0.2%, 0.2%, 0.4%). To fully exploit the advantages of charged-particle therapy, these uncertainties should ultimately all be minimized, no matter how large or small they actually are.

In this study, we focused on the uncertainties in patient CT imaging mainly caused by the acquisition of the CT numbers themselves. The CT number is directly related to the linear attenuation coefficient of the object for x-rays, and is usually calibrated to 0 for water and −1000 for air. The CT number is strongly affected by the x-ray energy spectrum at the point of measurement. The initial energy spectrum depends on the scanner properties, such as tube voltage, target material, filter, and detector sensitivity. The spectrum varies in the materials traversed by the x-rays up to the point of interest due to energy-dependent attenuations, namely beam hardening. The variation in CT number between scanners can be handled by creating a CT-number-to-SPR conversion function specifically for each scanning condition of each scanner. The variation in CT number due to the beam hardening is considered to be minimized by creating the CT-number-to-SPR conversion functions using x-ray calibration phantoms with typical body sizes, e.g. 10–40 cm [6-8]. However, the mismatch between the size of the object and the calibration phantom induces non-negligible effects in SPR of body tissues. Schaffner and Pedroni [9] reported that the CT-number variation of bone tissues between calibration phantoms of 15 cm and 30 cm diameters leads to an uncertainty in SPR of 1.5%. Yang et al. [6] reported that CT-number variations in bone and lung tissues between phantoms of 16 cm and 32 cm diameters lead to uncertainties of 1.9% and 2.6%, respectively. These studies, however, were solely based on the CT-number measurements of the tissue substitutes in cylindrical calibration phantoms. The impact of the beam hardening on the uncertainty of SPR estimation in a patient has never been investigated using a realistic patient geometry. In addition, the optimum size of the x-ray calibration phantom for minimizing the uncertainty has not been determined or proposed. The optimum phantom size, in sum, represents an ‘effective size’ of the patient cross section for patient SPR estimation.

In this study, we evaluated the uncertainty of SPR estimation due to beam hardening, and investigated the effective patient-cross-section size for patient SPR estimation in charged-particle therapy using an ideal CT system and realistic human tissue computational phantoms.

MATERIALS AND METHODS

X-ray CT

We modeled an ideal CT scanner consisting of an x-ray source generating parallel broad beams and a detector-array with a 100% detection efficiency with perfect antiscatter grids. In addition, CT scanning was simulated based on a theoretical x-ray spectrum with an infinite imaging dose. The reconstructed CT image was not affected by factors such as statistical noise, electronic noise or scatter, so we focused our investigation on the error in SPR estimation due to the beam-hardening effects.

The x-ray energy spectrum of a CT scanner Φ(E) was generated using the SpecCalc x-ray generator program [10]. For the generation, we used the following parameters: a normal tube voltage of 120 kVp with variations of ±10 kVp, tungsten as the target material, a 7° anode angle and a 7.4-mm-thickness aluminum filter. The generated energy spectra of the CT scanner are shown in Fig. 1.

Fig. 1.

The x-ray energy spectra Φ(E) of 120 kVp (solid curve), 110 kVp (dashed-curve) and 130 kVp (dotted curve) beams from the simulated CT scanner.

The simulated detector signal of the incident x-rays without an object was calculated as: where Emax is the maximum energy of the x-rays, i.e. 120 ± 10 keV. The measured signal of the x-rays transmitted through an object that enter the detector was calculated as: where μ(E, t) is the linear attenuation coefficient of the object at a position t along the projection line with length l. Projection of the object by the x-rays in a given direction was defined as:

To eliminate the effect of beam hardening by the polychromatic x-rays, λP of each projection line was corrected to the projection by a monochromatic x-ray with equivalent attenuation at the incident, λM, using a look-up table describing water-equivalent thickness [11]. Tomographic image reconstruction from a sinogram of scanned and corrected projections gave the mean attenuation coefficient of the object over the incident energy spectrum, . In this study, projections λM were generated for every 1° step, and a 2D filtered backprojection algorithm with a ramp filter was used for the reconstruction.

The values of were converted into CT number H in Hounsfield units (HUs): where and are the mean attenuation coefficients of water and air, respectively.

H-to-(S/Sw) conversion functions

To produce the CT number to SPR conversion functions, H-to-(S/Sw), we followed the procedure reported by Kanematsu et al. [1], based on the standard tissue data in the International Commission on Radiological Protection Publication 110 [12]. They defined 11 representative tissues of the human body with mass density ρ and elemental weights w and compiled their SPRs. When the CT numbers H of the 11 tissues are known, based on the hypothesis that an arbitrary tissue is a binary mixture of the representative tissues adjacent with higher and lower H, the SPR of the arbitrary tissue can be derived from its CT number H by polyline interpolation.

To determine the CT number H of the 11 tissues, we simulated cylindrical water phantoms of diameter L with a 2.5-cm diameter insert of one of these tissue materials. The energy-dependent attenuation coefficient μ of water, air, and the tissue materials were determined by interpolation from XCOM, a database of photon cross sections provided by the National Institute of Standards and Technology (NIST) at http://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html. The inserts at the center of the water phantoms were scanned by the CT scanner individually, and their CT numbers H were derived from the reconstructed CT images. For the CT-number calculation with equation (4), the values of and determined with a 25-cm diameter cylindrical water phantom were used throughout this study. The H-to-(S/Sw) conversion function for the phantom size L was produced by interpolating the CT number H and SPR of the 11 tissue materials.

To confirm the accuracy of the produced H-to-(S/Sw) conversion functions, 53 standard body tissues with ρ and w listed in the ICRP report were applied to the simulation of CT scanning as inserts of water phantoms with L = 10, 30 and 50 cm. The H-to-(S/Sw) relations of the standard body tissues determined for respective phantom sizes L were compared with the corresponding H-to-(S/Sw) conversion functions.

ICRP computational phantom

We used volumetric models of the reference male and female phantoms provided in the ICRP Publication 110 [12]. The height and mass of the models were (176 cm, 70 kg) and (167 cm, 59 kg), respectively. Their slice thickness (voxel height) and voxel in-plane resolution were (8.0 mm, 2.137 mm) and (4.84 mm, 1.775 mm). For better spatial resolution, in this study, the voxel in-plane resolutions were reduced to one-third of the original resolutions, i.e. 0.7123 mm and 0.5917 mm, respectively. Since charged-particle therapy is applied to a variety of tumors [13], the phantom data for pelvis, abdomen, thorax, and head and neck regions were included in the analysis, as shown in Fig. 4a. The phantom arms were removed to simulate the arm abduction used in treatment.

Fig. 4.

(a) Water-equivalent thickness distribution of the ICRP female phantom in the AP (PA) direction. White solid lines indicate the boundaries of anatomical regions, i.e. pelvis, abdomen, thorax, and head and neck. (b) The root mean square deviation and (c) the mean deviation of the SPR , at each slice position by the H-to-(S/Sw) conversion functions with L = 10 cm (black curves), 30 cm (red curves), 50 cm (green curves), 21 cm (blue curves), the region-specific Leffective (light blue curves) and the slice-specific Leffective (pink curves). (d) The effective patient-cross-section size Leffective determined at each slice position.

Reference SPR

The volumetric models of the male and female phantoms were converted to the SPR distributions by the following steps. The electron density ne of body tissues with a given ρ and w was derived by: where u = 931.5 MeV/c2 is the atomic mass unit, and Z and Ar are the atomic number and relative atomic mass for element i, respectively. The SPR, S/Sw, of the tissues was calculated using the Bethe–Bloch equation, which can be approximated by: where new = 3.343 × 1023 cm−3 is the electron density of water, mec2 = 0.511 MeV is the electron rest energy, β is the particle velocity relative to the speed of light in a vacuum, and I is the mean excitation energy for element i of solid or liquid compounds [14]. This led to Iw = 75.3 eV for water [15]. In equation (6), we used a fixed value of β2 = 0.135, which minimized the range errors in a patient for proton radiotherapy [16]. We refer to the SPR calculated by equation (6) as the ‘reference SPR’, and we symbolize it as (S/Sw)ref hereafter.

Predicted SPR

The volumetric models of male and female phantoms were scanned by the CT scanner slice-by-slice. The energy-dependent attenuation coefficients of 53 tissue materials μ were determined with XCOM. The reconstructed CT images represented in HUs were then converted to the SPR maps using the H-to-(S/Sw) conversion function for a phantom size L. We refer to such a derived SPR as the ‘predicted SPR’ and we symbolize it as (S/Sw)prd hereafter.

SPR error analysis

From the reference and predicted SPR distributions of the male and female phantoms, voxel-by-voxel absolute deviation of the SPR, δS = (S/Sw)prd − (S/Sw)ref, was calculated. The root mean square deviation of the SPR, , was derived to investigate the absolute error in the SPR estimation with the H-to-(S/Sw) conversion functions. The mean deviation of the SPR, , was also derived to estimate the range error of proton beams in a patient.

Effective size of patient cross section for H-to-(S/Sw) conversion functions

We varied the diameter of the x-ray calibration phantom L from 6 cm to 50 cm in 1 cm steps, and 45 H-to-(S/Sw) conversion functions were produced. For each conversion function, the δS map was obtained on each slice. We derived the root mean square deviation of the SPR in each body-part region (pelvis, abdomen, thorax, and head and neck) of the male and female phantoms. The average value of between the two phantoms was used to determine the effective size of the patient cross section for SPR estimation for each body-part region, Leffective. The effective size Leffective was also determined for respective slice positions of the male and female phantoms. To investigate the slice-by-slice variation of Leffective across the body-part regions, the standard deviation of Leffective, σeffective, was determined for the respective regions.

The uncertainty in the determined Leffective in each body-part region due to variations in body size was evaluated by one-half of the difference between the effective size of the male phantom Lmale and that of the female phantom Lfemale, i.e. ∆L ≡ (Lmale − Lfemale)/2. In addition, the uncertainty in the determined Leffective due to variations in initial energy spectrum of the CT scanner was evaluated by the deviations in the effective phantom sizes for the energy spectra with tube voltages of 110 kVp and 130 kVp from that with 120 kVp, ∆L110 and ∆L130.

RESULTS

Figure 2 shows the H-to-(S/Sw) conversion functions produced by the x-ray calibration phantoms of L = 10, 30 and 50 cm. Although there was a small discrepancy around −150 HU for adipose tissues, three conversion functions practically coincided with each other up to 100 HU for soft tissues. Subsequently, these functions deviated gradually with HU. Since beam hardening up to the point of the insert was milder in the smaller phantom, the slope of the conversion function for L = 10 cm was gentler than the slopes for L = 30 or 50 cm. The SPRs at 1400 HU, i.e. for the bone tissues, derived by the functions for L = 10, 30 and 50 cm were 1.61, 1.70 and 1.76, respectively. This implied that the error in the predicted SPR of bone tissues could be ±4% due to beam hardening, even if a moderate size of L = 30 cm was used for the conversion function. Further, the error in the predicted SPR could be ±6% for the teeth at ≈2700 HU.

Fig. 2.

H-to-(S/Sw) conversion functions produced with the x-ray calibration phantoms of L = 10 cm (black line), 30 cm (red line) and 50 cm (green line). Correlations between H and S/Sw of 53 ICRP body tissues determined with the phantoms were plotted by plus symbols (+) of different colors.

Correlations between H and S/Sw of 53 ICRP body tissues determined with the x-ray calibration phantoms of L = 10, 30 and 50 cm are plotted by colored plus symbols in Fig. 2. The ICRP tissues distribute around the polylines with RMS below 0.009 in S/Sw for all phantom sizes L.

SPR error

Figure 3 shows axial (S/Sw)ref distributions at typical slice positions in the body-part regions of the female phantom, and corresponding δS distributions by the H-to-(S/Sw) conversion functions of L = 10, 30 and 50 cm. The δS deviated from zero in bone and adipose tissues, and the degree of deviation was strongly dependent on the selected phantom size L. On the head plane, SPR of bone tissues was slightly underestimated by the conversion function of L = 10 cm, while it was overestimated by the functions of L = 30 and 50 cm. For L = 50 cm, δS amounted to +0.16 at the occipital region. The root mean square deviation and the mean deviation of the distributions were [0.021: −0.0029], [0.018: 0.0066] and [0.030: 0.0139] for L = 10, 30 and 50 cm, respectively. On the thorax plane, δS was the smallest for the function with L = 10 cm, which was identified by the smallest value of . The largest value of was observed for L = 50 cm. However, the absolute value of for L = 50 cm was smaller than that for L = 10 cm, since the overestimation of the SPR in bone tissues was well compensated for by the underestimation in adipose tissues of breasts. On the abdominal plane, δS was equally small for the functions with L = 10 and 30 cm, while it was the smallest for the function with L = 30 cm on the pelvic plane. and for each plane are described as legends in the figures.

Fig. 3.

Axial (S/Sw)ref distributions (upper row) on the head (left column), thorax (2nd left column), abdominal (2nd right column) and pelvic (right column) planes. The corresponding δS distributions by the H-to-(S/Sw) conversion functions with L = 10 cm (2nd upper row), 30 cm (middle row) and 50 cm (2nd lower row), and the effective size Leffective of the slices (lower row). The root mean square deviation and the mean deviation of the δS distribution are shown in square brackets, .

Figure 4a shows the water-equivalent thickness distribution of the ICRP female phantom in the AP (PA) direction. Figure 4b and c show variations of and with slice position in the female phantom, respectively. At all slice positions, was the largest for the conversion function with L = 50 cm. In the head and neck region, amounted to 0.047 due to the high absorption of x-rays in teeth and bone tissues. For the same reason, was also large in the head and neck region, and it amounted to 0.019 at the slice positions including the teeth.

Figure 5 shows the within the whole body and the body-part regions of the male and female phantoms for different L. The average between the two phantoms is also shown there. The variation of the average was different for different body-part regions, e.g. it reached 0.030 in the head and neck region for L = 50 cm, while it was ≤0.012 in the abdominal region for 6 cm ≤ L ≤ 50 cm. However, in all regions, the average showed a convex shape with respect to L. The phantom size corresponding to the minimum root mean square deviation was defined as the effective patient-cross-section size Leffective for constructing the H-to-(S/Sw) conversion function. Table 1 shows Leffective corresponding to the minimum for the whole body and the body-part regions. The effective sizes Leffective were slightly different for different body-part regions, while their variations were <±2.5 cm around Leffective = 21.0 cm of the whole body. The effective size Leffective determined for each slice position of the female phantom is shown in Fig. 4d. The variations of and with slice position of the female phantom by the conversion functions with L = 21.0 cm, with the Leffective of the respective body-part regions, and with the Leffective of the respective slice positions are shown in Fig. 4b and c, respectively. The conversion functions by the effective sizes considerably reduced the values of and , especially in the head and neck region. Figure 3(1e)−(4e) (the bottom row images) show axial δS distributions by the H-to-(S/Sw) conversion functions with Leffective determined for the slice positions. The σeffective quantifying the slice-by-slice variation of Leffective is also shown in Table 1. The x-ray attenuation in the lung tissues was insignificant due to their low mass density, while the attenuation in the blade bones was significant. The thorax region contained both of these tissues, inducing a large σeffective, i.e. ≥6 cm, in the region as shown in Fig. 4d. Figure 6 shows the within the whole body and the body-part regions of the male and female phantoms, and their average. The variations of the average for 6 cm < L < 50 cm were ≤0.004 for all of the body-part regions except for the head and neck region. The fraction of the bone tissues was high in the head and neck region, which resulted in the overestimation of SPR for large L, e.g. = 0.011 for L = 50 cm. In other regions, the overestimation of SPR in bone tissues was compensated for by the underestimation in adipose tissues, resulting in a moderate value of . The effective size deviations due to variation in body size, ∆L, and due to the variation in initial energy spectrum of the CT scanner, ∆L110 and ∆L130, are shown in Table 1. The deviation ∆L in the thorax region was ≥4 cm due to the significant differences in shoulder width and volume of the breasts between male and female models. In contrast, the deviations ∆L110 and ∆L130 were insignificant compared with ∆L.

Fig. 5.

Root mean square deviation of SPR, , for different phantom sizes L in (a) the whole body, (b) head and neck, (c) thorax, (d) abdominal and (e) pelvic regions of the male phantom (blue curves), female phantom (red curves), and their average (black curves). The solid, dashed, and dotted curves are for an x-ray tube voltage of 120 kVp, 110 kVp and 130 kVp, respectively.

Table 1.

Effective sizes of the patient cross section, Leffective, in the whole body, head and neck, thorax, abdominal and pelvic regions

	Leffective [cm]	σeffective [cm]	∆L [cm]	∆L110 [cm]	∆L130 [cm]
Whole body	21.0	4.8	2.5	−0.1	+0.1
Head/neck	18.8	4.7	1.5	+0.2	−0.2
Thorax	22.7	6.2	4.1	−0.2	+0.1
Abdomen	20.3	2.3	1.7	−0.3	+0.3
Pelvis	23.3	4.5	2.5	−0.3	+0.2

The standard deviation of Leffective across each body-part region, σeffective. The deviation in Leffective due to variation in body size, ∆L, and the deviations due to the variation in initial energy spectrum of the CT scanner, ∆L110 and ∆L130.

Fig. 6.

Mean deviation of SPR, , for different phantom sizes L in (a) the whole body, (b) head and neck, (c) thorax, (d) abdominal and (e) pelvic regions of the male phantom (blue curves), female phantom (red curves), and their average (black curves). The solid, dashed and dotted curves are for an x-ray tube voltage of 120 kVp, 110 kVp and 130 kVp, respectively.

DISCUSSION

The effective size of the patient cross section for minimizing the uncertainty in patient SPR estimation was determined for charged-particle therapy treatment planning using human-tissue computational phantoms provided by ICRP. The effective size was 21.0 cm for the whole body. Using the determined effective size, we were able to reduce the root mean square deviation of SPR to ≤0.011, and the mean deviation to ≤0.001 in the whole body.

In the human body, we found that the overestimation (or underestimation) of the SPR in bone tissues is often compensated for by the underestimation (or overestimation) in adipose tissues. The absolute value of realized by Leffective = 21.0 cm was in fact one order of magnitude smaller than the corresponding value of . The absolute value of was <0.001 in the pelvic region, while it was <0.004 in the head and neck region, with Leffective = 21.0 cm. These results indicate, for instance, that for a proton beam with 25-cm water equivalent length (WEL) range in the pelvic region and a proton beam with 12-cm WEL range in the head and neck region, the expected mean range errors are <0.03 and <0.05 cm WEL in their regions, respectively. Schaffner and Pedroni [9] reported larger range errors of 0.33 cm WEL in a prostate patient case and 0.14 cm WEL in a brain patient case for proton beams with similar ranges, i.e. 25 and 12 cm WEL, where they derived the range errors by linearly adding the expected range uncertainties of soft tissue and bone, corresponding to the calculation of , for their patient cases.

The slice-by-slice variation in Leffective in each body-part region quantified by σeffective was larger than the deviations ∆L, ∆L110 and ∆L130 as shown in Table 1. Further, σeffective was larger than the variation in Leffective among the body-part regions. These results suggest that the preparation of the H-to-(S/Sw) conversion functions according to the body-part regions is unnecessary for practical purposes. For simplicity, an x-ray calibration phantom with a fixed size of L = 21.0 cm should rather be used to produce the H-to-(S/Sw) conversion function, irrespective of the body-part regions, by which and are reasonably reduced in all body-part regions, as shown in Fig. 4. At the beginning of this study, we expected that the impact of the beam hardening on the uncertainty of SPR estimation in a patient can be minimized by selecting the optimum phantom sizes for calibrating H-to-(S/Sw) conversion functions. However, even with the conversion function derived from the optimum calibration for each slice, the predicted SPRs deviated from the reference SPRs. This may be the intrinsic limitation of a single-energy CT with polychromatic x-rays. The H-to-(S/Sw) conversion functions constructed with the selected representative tissues specifically for the body-part regions may potentially reduce the deviations. The beam-hardening correction based on the phantoms with more realistic compositions and configurations may also reduce the deviations. A dual-energy CT for patient SPR estimation may be another method for reducing the deviations [17-19].

There are several sources of uncertainty in patient SPR estimation, and one of them is insufficient correction of beam hardening arising from the mismatch between the size of the object and the calibration phantoms investigated in this study. To fully exploit the advantages of charged particle therapy, the remaining uncertainty sources should also be minimized.

CONCLUSION

One of the uncertainty sources in patient SPR estimation in charged particle therapy treatment planning is insufficient correction of the beam-hardening effect arising from the mismatch between the size of the object and the calibration phantom used. The uncertainty would be minimized by selecting a suitably sized x-ray calibration phantom in construction of the H-to-(S/Sw) conversion function, namely the effective size of the patient cross section. We determined the effective size using an ideal CT system and realistic human tissue computational phantoms provided by ICRP. The effective size was 21.0 cm, with which the root mean square deviation of SPRs and the mean deviation of SPRs could be reduced to ≤0.011 and ≤0.001, respectively, in the whole body. The effective patient-cross-section size of 21.0 cm can be used for the whole patient body, irrespective of body-part regions, for treatment planning in clinical practice.