If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Department of Oncology, University of Oxford, Oxford OX3 7DQ, United KingdomDepartment of Physics, University of Adelaide, Adelaide, South Australia 5005, Australia
Department of Physics, University of Adelaide, Adelaide, South Australia 5005, AustraliaDepartment of Medical Physics, Royal Adelaide Hospital, Adelaide, South Australia 5000, Australia
Protons undergo multiple Coulomb scattering, so require Bayesian path estimation.
•
Most likely path currently calculated assuming a water scattering environment.
•
We account for inhomogeneity of tissues using scattering power ratios.
•
Water assumption is valid for majority of medical scanning scenarios.
Abstract
Purpose
Multiple Coulomb scattering (MCS) poses a challenge in proton CT (pCT) image reconstruction. The assumption of straight paths is replaced with Bayesian models of the most likely path (MLP). Current MLP-based pCT reconstruction approaches assume a water scattering environment. We propose an MLP formalism based on accurate determination of scattering moments in inhomogeneous media.
Methods
Scattering power relative to water (RScP) was calculated for a range of human tissues and investigated against relative stopping power (RStP). Monte Carlo simulation was used to compare the new inhomogeneous MLP formalism to the water approach in a slab geometry and a human head phantom. An MLP-Spline-Hybrid method was investigated for improved computational efficiency.
Results
A piecewise-linear correlation between RStP and RScP was shown, which may assist in iterative pCT reconstruction. The inhomogeneous formalism predicted Monte Carlo proton paths through a water cube with thick bone inserts to within 1.0 mm for beams ranging from 210 to 230 MeV incident energy. Improvement in accuracy over the conventional MLP ranged from 5% for a 230 MeV beam to 17% for 210 MeV. There was no noticeable gain in accuracy when predicting 200 MeV proton paths through a clinically relevant human head phantom. The MLP-Spline-Hybrid method reduced computation time by half while suffering negligible loss of accuracy.
Conclusions
We have presented an MLP formalism that accounts for material composition. In most clinical cases a water scattering environment can be assumed, however in certain cases of significant heterogeneity the proposed algorithm may improve proton path estimation.
The main advantage of proton therapy over conventional X-ray radiation therapy is the Bragg peak in energy deposition, a phenomenon whereby most of a heavy charged particle’s energy per unit length is deposited toward the end of its trajectory in a sharp peak. This quality is advantageous in radiotherapy treatment using such particles as beam energy and intensity can be manipulated to deposit a highly conformable dose to the tumor volume, with a low dose on entry and no exit dose.
The range of protons in matter is dependent on the stopping power of the material through which it travels. To accurately predict the location of the Bragg peak, therefore, accurate knowledge of the stopping power of the patient tissues is required. Tissues can be characterised by their relative stopping power (RStP). This is simply the stopping power of a material relative to that of water at a given energy. For human tissues, the RStP is approximately independent of proton energy for therapeutic energies.
In current clinical practice, RStP is estimated by conversion of X-ray CT Hounsfield units via an empirically derived calibration curve [
]. Proton CT (pCT) is an alternative approach in which RStP of the patient is measured directly with an energetic proton beam.
In conventional CT imaging, X-rays are assumed to move in straight lines, whereas protons are known to undergo multiple Coulomb scattering (MCS), which poses a challenge for image reconstruction that is absent from X-ray CT imaging. In order to perform a pCT reconstruction, knowledge of the position, direction and energy of individual protons before and after they traverse the patient is required. Hanson et al. first proposed the concept of estimating proton paths by measuring their exit position post-patient [
], measurements of individual proton positions and directions pre- and post-patient are obtained through two 2-dimension sensitive tracking modules upstream and downstream of the patient. Each tracking module consists of orthogonally oriented single sided silicon strip detectors (SiSD) that are position-sensitive in one dimension. The energy of an incoming proton is assumed to be equal to the energy at which it is ejected from the accelerator, while the residual energy post-patient is measured by a segmented scintillation crystal calorimeter.
Assuming a straight line path in proton radiography leads to poor spatial resolution [
], so it is therefore necessary to develop a probabilistic model for the movement of protons through matter, known as a most likely path (MLP) formalism. In this paper we use the matrix-based MLP formula proposed by Schulte et al. [
]. In the Fermi-Eyges framework, Schulte et al. have applied Bayesian techniques to a bivariate Gaussian distribution to calculate the most likely lateral position and angular deflection at an intermediate depth, given the entry and exit conditions of the proton. The distribution provides an error matrix given by Eq. (27) in their paper, from which the lateral error can be calculated at any depth in the medium to define a probability envelope surrounding the most likely path. Typically in pCT iterative reconstruction, the (deterministic) path of a proton is mapped out by assigning binary values to each voxel; 1 if it passes through that voxel, and 0 otherwise. Using a probability envelope surrounding the MLP, it is possible to assign to each voxel a continuous value between 0 and 1 which represents the probability of the proton traversing that voxel. Wang et al. [
] demonstrated that using a 90% probability map directly for pCT reconstruction can yield a smoother reconstructed image when compared to using deterministic (binary) path maps. It is worth noting that the computation time in the method proposed by Wang et al. [
] stems from the fact that the compact MLP formula is derived for homogeneous media, with the covariance matrices calculated assuming all scatter takes place in water. Such an approximation may lead to inaccuracies in the proton path estimate during pCT reconstruction and result in substandard spatial resolution in the reconstructed image. A study by Wong et al. investigated the effect of tissue inhomogeneity on the accuracy of the water-based MLP estimate [
]. They found the error in proton path reconstruction to be up to 20% larger in a 20 cm slab geometry which included small bone inserts and air cavities. There was no attempt, however, to adjust the MLP formalism to account for these materials.
In this paper, we propose an inhomogeneous MLP formalism, based on the ratio of scattering power in the material of interest to that in water. We refer to this quantity as the relative scattering power (RScP) hereon. Using the scattering power formulae of Gottschalk [
], the RScP at any depth is dependent only on the material at that depth, even when considering nonlocal effects of MCS buildup in the material preceeding the calculation point. It is therefore straightforward to calculate the Fermi-Eyges scattering moments [
] in inhomogeneous media when the composition is known a priori. Furthermore, we are able to investigate a relationship between RScP and RStP in different materials. In a 2017 study, Collins-Fekete et al. [
], however this requires integrating a material dependent quantity, the radiation length, over depth. Gottschalk’s formulae are both more accurate and better suited to mixed slab geometries than the differential Highland formula used by Kanematsu [
In an attempt to improve computational efficiency, several methods using cubic spline trajectories have been proposed as an alternative to Bayesian calculation of the MLP [
]. In this work, we introduce an inhomogeneous MLP-Spline-Hybrid method which involves the use of cubic splines fit through a limited number of Bayesian MLP points.
]. We show a piecewise-linear relationship between RStP and RScP for a range of human tissues and describe how this relationship may be implemented in iterative pCT reconstruction.
2. Methods
In the compact matrix-based MLP formula of Schulte et al. [
] the most likely lateral position and angular deflection at an intermediate depth are represented by the vector , given the respective entry and exit conditions, at depth , and at depth . The MLP is calculated by Eq. (24) in their work [
where is the scattering power, which is defined as the rate of increase, with depth u, of the mean square of the projected scattering angle . That is,
(7)
The variance in lateral displacement is given by the (1,1) matrix element of (8);
(8)
This variance reflects the square of the standard error in the MLP and may be used to define a probability envelope surrounding the most likely path. This probability envelope is defined such that it contains at least a specified proportion of all possible MCS influenced paths. For example, given a pair of entry and exit information, only 0.3% of all possible particle trajectories should exit a envelope surrounding the MLP. Note that (8), and therefore the size of the envelope, is dependent on the material composition of the medium traversed by the particle.
2.1 Scattering power calculations for inhomogeneous materials
Here, we present a method for directly calculating the depth dependent scattering power , and hence the scattering moments, (4) to (6), in various materials. This method does not require the empirical approximations from Lynch and Dahl [
In the Fermi-Eyges framework, in which a Gaussian scattering process is assumed, scattering power may be calculated from first principles using the scattering cross-section. Using the cross-section formula of Goudsmit and Saunderson [
] formulae, Gottschalk derived an expression for scattering power, which could be considerably simplified for radiotherapy protons in the energy range of 3 to 300 MeV [
], by first defining a quantity, he termed the scattering length, . For a scattering medium with mass density , atomic number Z, and mass number A, the scattering length is defined as
(9)
Here, is the fine structure constant, is Avogadro’s number, e is the charge of an electron, and is the rest mass of an electron. Using this scattering length, it can be shown that the scattering power obeys
(10)
where is the depth-dependent kinetic energy of the proton and its reduced kinetic energy, , is defined as
(11)
where is the rest mass of a proton. In addition to Gottschalk’s work [
] for clarification of the underlying theory leading to this result.
In order to calculate the scattering power for composite materials, such as those found in the human body, we notice that the only material dependent part of (10) is the scattering length. Thus, if the composite material consists of n elements, each with fractional weight per volume then
(12)
where is the mass density of the composite material.
We wish to catalogue materials based on a single energy-independent value, analagous to the RStP, which can be assumed valid at any depth. This can be achieved by defining the relative scattering power (RScP), , as the ratio of the scattering power for the material to the scattering power for water at the same energy. Using Gottschalk’s formula (10), the energy dependent terms cancel and the result is simply the ratio of the scattering lengths. It can be shown that
(13)
with weights and mass density (in units of g/cm3) as defined in (12). The factor of 0.01636 is the reciprocal of the summation term applied to liquid water. In practice, we can now calculate the scattering power at any depth, and hence the covariance matrices in (3) using
(14)
where the subscript ‘w’ refers to the value for water.
2.2 Stopping power calculations for inhomogeneous materials
If we wish to calculate scattering power, and hence the scattering moments, using (14) then we must have an estimate of the kinetic energy of the proton at depth u. This can be calculated through an appropriately weighted combination of the forward and backward Euler methods, using stopping powers.
One can rearrange the definition of stopping power,
(15)
to find the energy lost per unit length if the stopping power is known. In discretizing the depth we have for for some positive integer N. The spacing between successive depths is constant and equal to for . The forward Euler method gives
(16)
with a boundary condition of , the incoming energy of the proton as it enters the medium, assumed to be equal to the exiting energy from the accelerator. The backward Euler method gives
(17)
with a boundary condition of , the outgoing energy of the proton as it exits the medium, as measured by a detector.
Errors will accumulate with successive steps in each method. Hence it is reasonable to perform a weighted average of both methods. For example, if more steps must be taken in the backward method than in the forward, the forward method is expected to be more accurate and will be weighted more heavily. is multiplied by the normalised length between and , while is multiplied by the normalised length between and . That is,
(18)
with boundary conditions and .
Stopping powers can be determined using the Bethe-Bloch formula [
where is the electron number density of the material, I is the mean excitation potential, is the vacuum permittivity, and where v is the speed of the proton and c is the speed of light in a vacuum. Experimental values of the mean excitation potential for a single-element substance (e.g. O2 gas) can be obtained from a database, but for more complicated composite materials consisting of n elements, one may use the Bragg additive rule,
(20)
where
(21)
is the fractional weight per volume of the k-th element, as in (12).
Similar to relative scattering power, we may define the relative stopping power (RStP), , as the stopping power in the material, S, divided by the stopping power in water, , at the same energy. That is,
(22)
Note that the RStP is weakly energy dependent, however the effect for most materials and energies encountered clinically can be assumed negligible. This is explored in more detail in Section 2.4. Now, if the material through which the proton is passing is known, then one may calculate the stopping power by simply multiplying the stopping power in water by the RStP for that material,
(23)
In the proposed use of this method, the current estimate of RStP from the reconstructed pCT image will be used to approximate RStP in each voxel. The forward and backward Euler methods can now be alternatively written, respectively, as
(24)
(25)
where is the relative stopping power at depth .
The calculation of the MLP is implemented in a two-stage process, assuming there exists a catalogue of RScP and RStP values for all materials of interest (refer to Section 2.4). The first stage involves estimating the proton energy at discrete depths within the patient, using (18), which is then used to calculate the material-dependent scattering power, using (14). The second stage is to calculate the scattering moments, (4) to (6), using numerical integration, and therefore calculate the most likely position at each discrete depth using (1). More algorithmic details of our implementation can be found in Section 2.6 and in the Appendix.
2.3 MLP-Spline-Hybrid method
In an effort to increase computational efficiency, the inhomogeneous MLP may be calculated using the aforementioned methods only at subsampled depths within the phantom. The entire trajectory is then estimated by fitting a cubic spline through these points, given the initial and final directions of the proton as measured by detectors placed pre- and post-patient. The upper and lower bounds of the probability envelope are constructed by fitting cubic splines through the points generated by adding and subtracting, respectively, the lateral error, using (8), to the MLP estimate at the calculation points. A multiplier can be applied to the lateral error to set the significance level. For example, a multiplier of 3 creates a 3 envelope which should contain 99.7% of possible proton trajectories for given entry and exit information.
We propose that the sampled points be chosen as distinct material boundaries, so that the number of samples used for fitting the spline scales with the complexity of the heterogeneity. These boundaries may be identified by thresholded RStP differences between air, tissue and bone. We refer the reader to Section 4 for a discussion on how this may be implemented in the context of iterative pCT reconstruction.
The method described above shall be referred to as the inhomogeneous MLP-Spline-Hybrid method and denoted SH hereon.
2.4 Catalogue of anatomical materials
We have thus far explored methods for calculating the relative scattering and stopping powers in various materials. In iterative pCT reconstruction, one may first assume the imaged object is made purely of water, and on each successive iteration the internal composition may be updated. Prior to the reconstruction, the depth at which some material is located is unknown, and it is for this reason that we wish to characterise materials by energy independent RStP and RScP values (as proton energy is related to penetration depth). We have already observed that our expression for RScP in (13) is dependent only on material composition, and not on the proton energy, which is a direct result of using Gottschalk’s derived scattering power [
], which gave their elemental composition, mass density and electron number density. Values were also calculated for air, which has been approximated as consisting 79% of Nitrogen and 21% of Oxygen. ( g/cm3 has been used as the electron number density.) Mean excitation values used in the calculation of stopping power were obtained, for the same set of materials, from ICRU Report 49 [
]. The nominated RStP assigned to each material was arbitrarily defined as the average of the energy dependent RStP calculated at 1 MeV incremeents over the range MeV.
Arbor et al. quote a region of negligible energy dependence between 80 and 300 MeV, by comparing GEANT4 RStP values at various proton energies to that at 300 MeV [
]. They used a clinical Gammex 467 tissue characterisation phantom, for which the material contributing to the greatest RStP variation (approximately 1%) was the cortical bone equivalent, Gammex SB3. Our calculations show that our nominated RStP value does not differ from the true value by more than 1% for 30 to 300 MeV protons, with the largest error manifested in cortical bone. Further, the error is less than 2% over the entire range of 10 to 300 MeV, and the residual energy of a proton exiting a patient is generally larger than both of these lower limits.
2.5 Correlation of relative scattering and stopping powers in human tissues
As mentioned above, in an iterative pCT reconstruction algorithm, the RStP values in each voxel may be updated on successive iterations to build an image of the body. Implementation of a calibration curve that determines RScP from RStP could provide an improvement to the convergence of the algorithm and the final accuracy of the image.
In Fig. 1 we see strong evidence of a linear relationship between and for the solid tissues listed in ICRU report 46 [
]. We use the adipose tissue with the lowest stopping power to define the intersection point (shown by the red marker) between the low density materials and solid tissues. A linear relationship is also expected for the low density materials as the chemical composition will be roughly equivalent and thus it is only the density increase in both parameters defining the correlation.
Fig. 1Correlation between relative stopping and scattering powers ( and respectively) for the materials listed in ICRU Report 46
. Linear least squares fits are shown for solid tissues and low density materials separately, with the intersection point (in red) defined at adipose tissue (Adipose Adult #3 in
A linear least-squares model is applied to the data. In this method we assume
(26)
where are constant parameters, given in Table 1, and is constrained to zero for low density materials.
Table 1Parameters for the least squares fits, , on solid tissues and low density materials. is the correlation coefficient (with indicating a good quality of fit for a sufficiently large data set). is the number of samples in a group.
The performance of the inhomogeneous MLP formalism, , proposed in this paper was evaluated by comparison with proton tracks produced by Monte Carlo simulation using TOPAS, version 3.1.2 using the standard environment and physics lists [
] in three phantom setups, a homogeneous water cube, a water cube containing thick transverse slabs of bone, and a clinically relevant human head phantom. The effect on the performance of both algorithms as the nominal beam energy is decreased has also been investigated.
A monoenergetic planar fan beam of 200 MeV protons was delivered to a 20 cm cube of homogeneous water such that the distance between the particle source and the surface of the phantom was 160 cm. The fan beam was parameterised such that the beam exactly spanned the exit face of the phantom, as shown in Fig. 2. Proton histories, including energy, position and direction were collected at 5 mm depth increments. These query depths should not be confused with the step-size in the Monte Carlo simulation itself. The step-size and included physics processes were unaltered from the default settings in TOPAS, which have been optimised for proton therapy applications [
]. Increasing the sampling rate for the MLP calculation (at increments smaller than 5 mm), and therefore decreasing the step-size in the Euler method, (24) and (25), and in the numerical integration, (4) to (6), did not yield any significant improvements in accuracy in the case of the slab phantom.
Fig. 2Geometric setup in TOPAS of a fan beam of protons incident on a 20 cm cube of homogeneous water. The beam exactly spans the distal surface of the phantom.
The entire setup was placed in a vacuous world volume and two vacuum sensitive detectors were placed pre- and post-phantom to record the position, direction and energy of each proton before and after it traversed the cube. The recorded proton positions were projected onto the entrance and exit surfaces of the phantom prior to the MLP calculation. In an effort to retain only Gaussian-natured MCS events, data cuts were performed to remove tracks for which the total energy loss or projections of the relative exit angle onto both the t-u and u-v planes (refer to Fig. 2) differed from their mean values by more than a chosen number of standard deviations. Once the data cuts were performed, the MLP was only calculated for the projection of the path onto the t-u plane, as the scattering in orthogonal directions form two independent probabilistic processes [
The above procedure was repeated for a 20 cm inhomogeneous cube consisting of 2 cm of water, 7 cm of cranium bone, 2 cm of cortical bone, another 7 cm of cranium bone and another 2 cm of water (in this order). Nominal beam energies of 230, 225, 220, 215 and 210 MeV were tested to investigate the robustness of and as the overall likelihood of scattering is increased as a result of decreasing the beam energy. Implementation of for beam energies other than 200 MeV requires recalculation of the fifth order polynomial for the energy factor, , as detailed by Schulte et al. [
]. The average proton energy was recorded at 5 mm increments in the 20 cm water cube for each nominal beam energy, which was then used to calculate the energy factors.
Finally, the slab phantom was replaced with a human head phantom to better resemble a clinical scenario. A 10-year-old male phantom (Computational Human Phantom Series, U.S. National Institutes of Health, National Cancer Institute, 2019) [
]. The beam and detector setups were unchanged from the water phantom setup. A monoenergetic 200 MeV beam was delivered to the phantom and proton histories were inquired at 1 mm depth increments, due to the increased complexity of the geometry. The phantom spanned up to 143 mm in the beam direction. Each slice along the longitudinal axis is separated by 2.425 mm and each pixel in the transverse plane has dimensions of 0.99 mm by 0.99 mm. The full setup is shown in Fig. 3, which includes sample proton tracks (prior to statistical cuts) for illustrative purposes. Further data cuts were performed to only include particle tracks that traversed the phantom.
Unlike the transverse slab geometry, for which there is only one material present at depth u, the material composition of a realistic phantom is voxelised and must be treated as such in the implementation of and SH. For each proton registered at the prior and posterior detectors, the entry and exit positions and directions were used to fit a three-dimensional cubic spline. This approximate path was subsequently used to interpolate the RStP in the traversed voxel at depth u. The RStP was pre-calculated for each CT coordinate using the HU-RStP calibration of Schneider et al. [
]. The RScP was then calculated using (26), allowing the calculation of the inhomogeneous MLP. (Refer to Algorithm1 and Algorithm2 in the Appendix for more details on the implementation.) We note that while better accuracy may be attained by using the material composition directly in the calculation of RStP and RScP, our procedure has been designed to mimic the approximation of both values through the iterative reconstruction process, as we suggested in Section 2.5.
Our application of the Spline-Hybrid method to iterative reconstruction of voxelised, non-slab geometries involves thresholding the current image estimate with three (or more) predefined regions corresponding to bone-like materials, tissue-like materials and air. An edge detection filter can be applied to the current iterate of the image based on the differences between RStP values in these regions, in much the same way that the outer patient boundary is detected on the first iteration in current pCT reconstruction methods [
]. For instance, the categories of solid tissues and low density materials defined in Fig. 1 may be subdivided further to define a greater number of regions as the number of iterations increases. Here, our implementation of an edge-detection filter resulted in an average subsampling rate of 10% of the available SH calculation points. An example of this subsampling is shown in Fig. 3 (bottom right).
Fig. 3Sample proton tracks (prior to statistical cuts) through a human head phantom (Computational Human Phantom Series, U.S. National Institutes of Health, National Cancer Institute, 2019) [
]. Bottom right: Hounsfield units and RStP profiles for the particle track pictured in blue. An example of query points for implementing the MLP Spline-Hybrid method is shown. The depth coordinate X here is equivalent to u, which is used elsewhere in this paper. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
When comparing our inhomogeneous formalism to the homogeneous water approach, a hull-detection method was used prior to the calculation of , as this is routinely performed in practice. The inhomogeneous formalism required no hull-detection step, as the air outside the patient is automatically considered during the MLP calculation.
3. Results
3.1 Water phantom
The error in the MLP as a function of depth, shown in Fig. 4, demonstrates that the inhomogeneous formalism performs as required, giving the same outcome in water as the homogeneous formalism when data cuts are applied to minimise the effects of elastic nuclear collisions and large-angle MCS. After cuts are performed on the relative exit angle and total energy loss, the maximum RMS error in lateral displacement is approximately 0.56 mm. A suitably accurate MLP estimation should give not only a good estimation of the trajectory, but also of the surrounding probability envelope. Indeed, Table 2 shows that at least 99.7% of particle tracks fall within the 3 envelope, as expected, for both formalisms.
Fig. 4Root mean square (RMS) error in lateral displacement in a 20 cm water cube over 3000 proton histories. is compared to for no data cuts, 2 cuts and 3 cuts in both relative exit angle and total energy loss.
Table 2Percentage of proton tracks that escape the MLP envelope for the water and slab phantoms. This measure provides a point of comparison between the homogeneous and inhomogeneous MLP path predictions for different levels of data cuts in both relative exit angle and total energy loss.
The inhomogeneous formalism when applied to 210 MeV protons traversing the slab phantom achieved a 17% reduction in maximum RMS lateral position error when compared to the homogeneous formalism. This can be seen in Fig. 6(a) and Fig. 5. Both algorithms had a runtime of approximately 3 ms per particle track using MATLAB (MATLAB R2017b, The MathWorks, Inc., Natick, Massachusetts, United States) on a standard Intel Core i7 processor. The spline-hybrid method, SH, regularly achieved very similar results up to 3 times faster on average. It can be seen in Fig. 5 that at lower energies, at which scattering is more pronounced, SH retains greater accuracy than . After applying 2 statistical cuts, the percentage of tracks falling outside the 3 envelope, listed in Table 2, was consitently less than 0.3% using the inhomogeneous formalism, indicating that the probabilistic path estimation is suitably accurate. On the contrary, significantly more tracks escaped the envelope, with only 90.4% being contained for the 210 MeV beam. An example proton track is shown in Fig. 6(c) along with all three MLP estimates. It is clear that and SH give a larger and more right-skewed probability envelope surrounding the proton path when compared to for this particular track.
Fig. 5(a) Percentage reduction in the root mean square (RMS) error in lateral position when SH is chosen over for various nominal beam energies incident on the slab phantom. (b) Maximum percentage error reduction with each nominal beam energy.
Fig. 6(a,b) Root mean square (RMS) error in lateral displacement in (a) the slab Phantom over 3000 proton histories, and (b) the human head phantom over 10,000 proton histories. is compared to and SH after applying 2 cuts in both relative exit angle and total energy loss. In the case of the human head phantom, has also been applied following a hull-detection step. (c,d) Examples of proton tracks in (c) the slab phantom, and (d) the human head phantom, showing the actual Monte Carlo (MC) path, the path as inferred by the different MLP methods, including their surrounding 3 probability envelopes. Note: error bars have been omitted for (b) as they are neglible.
There was no improvement in accuracy by employing the inhomogeneous MLP formalism for the human head phantom when compared to the homogeneous method with a prior hull-detection, denoted (Hull). This is evident in Fig. 6(b). It can be seen in Fig. 6(d) that there is also neglibile difference between the surrounding envelopes of the inhomogeneous formalism and (Hull). Both and (Hull) recorded runtimes consistently below 24 ms per track, however SH was able to complete the task in half this time.
4. Discussion
A new inhomogeneous most likely path formalism using the matrix methods of Schulte et al. [
] has been used to calculate the relative scattering power (RScP) values of protons between 3 and 300 MeV for various human tissues. The relative stopping power (RStP) was calculated for this same energy range, showing negligible energy dependence (less than 1% for energies above 30 MeV, and less than 2% for those above 10 MeV). We must note, however, that the water equivalent thickness of the slab phantom is not much smaller than the range of 210 MeV protons. Consequently, an average residual energy of 15.7 MeV was recorded post-phantom and 19.7% of protons had a residual energy of less than 10 MeV (after 2 data cuts were applied). No other beam energies resulted in a residual energy of less than 10 MeV, however 5.25% of recorded protons from the 215 MeV beam had an energy of less than 30 MeV. It was necessary to determine to what degree the improvement in MLP accuracy was due to the inhomogeneous scattering moments, rather than simply the energy estimate via the Euler method given in ((24), (25)). We therefore recalculated (and SH) using a 5-th order polynomial estimate of the mean energy, akin to the energy factor used in [
]. There was no noticable difference (less than 0.3%) in the RMS MLP error between these two methods, which leads us to conclude that the improvement in MLP accuracy was almost entirely due to the inclusion of inhomogeneity in the scattering moments.
A bi-linear correlation between RScP and RStP has been shown, which may be implemented in iterative pCT reconstruction as the RStP in each voxel is updated on successive iterations. Using the inhomogeneous formalism in a 20 cm cube of water, gave the same path as but with a slightly larger probability envelope when comparing estimates to Monte Carlo data obtained using TOPAS [
]. Note that the implementation of requires calculation of the energy factor polynomial , and therefore simulation of protons traversing the required depth at the correct nominal beam energy. Many polynomial coefficients may be pre-calculated and stored, however this is not required when calculating .
In the inhomogeneous slab geometry consisting of large slabs of cortical and cranium bone, achieved greater accuracy than by up to 17% at the lowest beam energy of 210 MeV. It was shown that the improvement in accuracy gained by using over increased as the beam energy was decreased and scattering in the phantom was more pronounced. This robustness in the inhomogeneous formalism has potential to be valuable in situations where lower energy beams are used. For instance, lower energy pencil beams are advantageous when aimed toward small volumes of matter near the body boundary, to ensure that the Bragg peak occurs within the absorbing detector post-patient. Furthermore, density resolution in pCT imaging has been shown to increase when using lower energy protons [
has in general outperformed in cases of significant inhomogeneity. It is important to note that a suitably accurate MLP formalism should accurately infer both the MLP itself and the probability envelope surrounding it. After applying 2 data cuts, all inhomogeneous MLP estimates in this study contained at least 99.7% of tracks within the 3 probability envelope, as expected. This was not the case when applying the water assumption to the slab geometry, as the number of tracks falling outside the envelope was significantly greater than 0.3% for all nominal beam energies.
Throughout this paper we have referred to Gottschalk’s method for calculating material-dependent scattering power, assuming locality [
]. That is, the scattering power depends on only local properties of the material and the proton energy. In reality, one must consider the buildup of MCS events in material traversed prior to the point of calculation. Gottschalk defined a nonlocal extension to (10) which differs by a simple logarithmic single scattering correction factor, , depending only on the initial energy of the proton, , and its energy at depth u. Namely,
(27)
We have adopted the subscript “dM” (for “differential Molière”) from Gottschalk’s work. Gottschalk showed that this nonlocal formula outperformed the differential Highland formula presented by Kanematsu [
] is based. It was also shown to be more easily implemented in mixed slab geometries as the non-locality parameter is not material dependent. Contrastingly, Kanematsu’s method requires numerical integration over a material-dependent quantitiy, and its accuracy is thus subject to one’s confidence in the material composition preceeding the point of interest, and the chosen step-size for numerical integration. It was for these reasons that Gottschalk’s method was chosen for calculating the scattering power in our proposed MLP formalism.
Using (27) we still obtain the same expression (13) for the RScP, however the scattering power in water differs in (14). We found negligible (less than 1%) difference in terms of MLP RMS error and the percentage of paths falling outside the 3 envelope when implementing the non-local correction to the scattering power. Thus, the local version (10) was chosen for simplicity. This is due primarily to the fact that the proton entry and exit information are known within the Bayesian MLP formalism. It was found that, while the non-local correction is substantially large at energies that are either very low or very close to the input energy, the depths at which a proton exhibits these energies are near either the point of entry or exit. Elsewhere in the phantom, the effect of the correction is very close to null.
In terms of designing a practical pCT scanner, the intended energy range of the protons determines the required water equivalent thickness (WET) of the residual range telescope (RRT), or calorimeter. For example, a residual energy of 115 MeV post-patient will result in a WET of approximately 10 cm. In initial designs, a crystal calorimeter was prescribed a depth appropriate for stopping protons with an initial energy of 200 MeV during a head scan [
], an RRT with a WET of 55 mm was used, corresponding to a maximum residual energy of approximately 80 MeV. Increasing the incident proton energy, in addition to requiring a physically larger RRT, can lead to increased uncertainty in the measured range [
] included only simulations of 330 MeV protons with a head, thorax and abdominal phantom. While this proton energy is favourable in terms of reduced MCS, accurate residual energy detection, particularly for a head-sized object, would be challenging with current design concepts. We have assessed our inhomogeneous MLP algorithm with proton energies more suited to residual range detector designs for proton imaging of head-like objects.
The detectors in our slab phantom setups were placed 1 cm from the scoring geometry on either side. A linear projection from the detector to the phantom in the direction of motion was made to determine the entry and exit positions at the scoring boundary, as negligible scatter contribution is expected as the proton traverses air.
It is important to note that large volumes of high density material such as those in the slab phantom will rarely be encountered in clinical practice, and so a gain in lateral position accuracy alone from may be insignificant in many cases. However, as the variance in scattering probability is increased in higher density materials, the width and skewness of the probability envelope surrounding the MLP within the inhomogeneous formalism are variable and better representative of the confidence in the expected lateral position. This quality has the potential to be useful in probabilistic reconstruction (see, for example, [
]) of high density materials, if the volume is sufficiently large or the beam energy is sufficiently low. A right-skewed envelope is evident for the slab geometry tested in this work, which can be explained by the fact that protons, on entry, have larger momentum and therefore lower scattering power. Scattering is more pronounced toward the distal part of the phantom where the particle has lower kinetic energy.
In this study, we have only considered inhomogeneous slabs transverse to the initial beam direction which extend to the boundaries of the phantom. It would be beneficial to explore the effect of structures both finite in size and parallel to the beam direction on the accuracy of the inhomogeneous MLP. In such a setup, discontinuities in the scattering of protons at the boundary could be expected to introduce a systematic error in the MLP under the homogeneous water assumption [
Clinical applicability of our formalism was investigated by applying these methods to a human head phantom. There was no noticeable increase in path estimation accuracy using the inhomogeneous formalism for a beam of monoenergetic 200 MeV protons, when compared to the standard water assumption with a prior hull-detection step. Our results indicate that in the majority of medical scenarios, the assumption of a water scattering environment is valid.
The MLP calculation is the most computationally intensive aspect of iterative pCT reconstruction. Both and (Hull) performed similarly in terms of computational time, however SH was twice as fast on average, with no appreciable loss in accuracy. The speed of the algorithm may be improved further through redefining the subsampling technique, using an alternative hull-detection technique, or interpolating the RStP voxel values using a straight line path instead of a cubic spline. However it should be noted that techniques such as parallelisation are able to greatly improve the computational efficiency of the MLP calculation under the assumption of a water scattering environment, and a prior hull-detection step significantly reduces the number of calculation points. Our proposed inhomogeneous formalism has the potential advantage of improving convergence for air inside the body, and while a hull is suited well to a convex patient exterior, our method may be able to handle non-convexities on the exterior better than current methods such as filtered back-projection.
5. Conclusions
A new formalism for calculating the most likely path of protons through inhomogeneous matter has been proposed, based on the compact matrix methods of Schulte et al [
]. The inhomogeneous formalism was shown to predict Monte Carlo proton paths to within 1.0 mm on average for beams ranging from 230 MeV down to 210 MeV in a 20 cm phantom consisting of thick slabs of cortical and cranium bone. The improvement in accuracy was most noticeable at lower energies, ranging from 5% maximum RMS error in the MLP for a 230 MeV beam to 17% for 210 MeV, when compared to their corresponding Monte Carlo tracks. Implementation of a new MLP-Spline-Hybrid method offered reduced computation time while suffering negligible loss of accuracy. There was no accuracy improvement for simulated 200 MeV protons through a more clinically relevant human head phantom, suggesting that a water scattering envrionment can be assumed in most clinical cases. However, in certain cases of significant heterogeneity the proposed algorithm may improve proton path estimation.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
M. Brooke is supported by a John Monash Scholarship and a Clarendon Fund Scholarship. This work was supported by Cancer Research UK Grant No. C2195/A25197, through a CRUK Oxford Centre DPhil Prize Studentship. The authors acknowledge Dr. Choonsik Lee and the National Cancer Institute for providing data sets from the Computational Human Phantom Series (U.S. National Institutes of Health, National Cancer Institute, 2019).
Appendix A. Implementation of inhomogeneous most-likely-path algorithm
Algorithm1 followed by Algorithm2 provides one method of implementing the inhomogeneous MLP formalism proposed in this work. Relative scattering powers and relative stopping powers for materials of interest may be calculated using the methods outlined in Sections 2.1, 2.2, respectively. A vector of ordered discrete depths is passed to the MLP calculation function, along with corresponding vectors of RStP and RScP values at each depth, and the entry and exit information of the proton. Algorithm1 calculates the material-dependent scattering power, which is then used by Algorithm2 to calculate the most likely lateral position at each depth and the associated standard deviation in this position. The scattering moments in (4) to (6) are calculated using the trapezoidal rule for numerical integration, however alternate numerical integration methods are possible. The step-size for the numerical integration or for the Euler method can also be increased, for improving calculation time, or decreased, for improving accuracy. Experimentation showed insignificant accuracy improvement from using step-sizes smaller than 5 mm for the inhomogeneous slab phantom in this work. Therefore, both algorithms use a step-size that is the same as the depth sampling rate.
Tabled
1
Algorithm1 Scattering power calculation along path.
30028-4&id=gr1.jpg){kind=link}
30028-4&id=gr2.jpg){kind=link}
30028-4&id=gr3.jpg){kind=link}
30028-4&id=gr4.jpg){kind=link}
30028-4&id=gr5.jpg){kind=link}
30028-4&id=gr6.jpg){kind=link}