An inhomogeneous most likely path formalism for proton computed tomography

Highlights • 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.


Introduction
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 [1]. This approach can lead to errors in stopping power of up to 3% [2,3]. 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 [4][5][6]. In the pCT system concept and design by Schulte et al. in 2004 [7], 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 [8], 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. [9], which followed the methods of Williams [10]. 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. [11] 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. [11] is increased significantly due to calculations involving non-sparse matrices.
The incompleteness of the formalism by Schulte et al. [9] 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 [12]. 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 [13], 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 [14] 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. [15] presented an inhomogeneous MLP formalism in which the Fermi-Eyges moments are calculated using a method proposed by Kanematsu [16], 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 [13].
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 [17][18][19]. 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 have used the TOPAS tool [20], which wraps and extends the GEANT4 simulation toolkit [21], to compare our new inhomgeneous MLP formalism to that of Schulte et al. [9]. 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.

Methods
In the compact matrix-based MLP formula of Schulte et al. [9] the most likely lateral position t 1 and angular deflection 1 at an intermediate depth u 1 are represented by the vector = t y ( ) T at depth u 2 . The MLP is calculated by Eq. (24) in their work [9], where R 0 and R 1 are change-of-basis matrices, and 1 and 2 are the covariance matrices, The elements of the covariance matrices in (3), known as the scattering moments, are given by (for where T u ( ) 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, The variance in lateral displacement is given by the (1,1) matrix element of (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 3 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.

Scattering power calculations for inhomogeneous materials
Here, we present a method for directly calculating the depth dependent scattering power T u ( ), and hence the scattering moments, (4) to (6), in various materials. This method does not require the empirical approximations from Lynch and Dahl [22] that are used in the homogeneous formalism of Schulte et al. [9].
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 [23,24], which is better behaved at small angles than the well-known Mott [25] and Rutherford [26] 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 [13], by first defining a quantity, he termed the scattering length, X 1/ s . For a scattering medium with mass density , atomic number Z, and mass number A, the scattering length is defined as Here, is the fine structure constant, N A is Avogadro's number, e is the charge of an electron, and m c e 2 is the rest mass of an electron. Using this scattering length, it can be shown that the scattering power obeys where E u ( ) is the depth-dependent kinetic energy of the proton and its reduced kinetic energy, , is defined as where m c p 2 is the rest mass of a proton. In addition to Gottschalk's work [13], the reader is referred to Rossi's book [27] 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 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), T , 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 with weights w k and mass density (in units of g/cm 3 ) 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 where the subscript 'w' refers to the value for water.

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, to find the energy lost per unit length if the stopping power is known. In discretizing the depth we have u j for = … j N 0, , for some positive integer N. The spacing between successive depths is constant and equal to with a boundary condition of = E E F 0 in , 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 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. E j F is multiplied by the normalised length between u j and u N , while E B is multiplied by the normalised length between u 0 and u j . That is, where e is the electron number density of the material, I is the mean excitation potential, 0 is the vacuum permittivity, and v c / 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. O 2 gas) can be obtained from a database, but for more complicated composite materials consisting of n elements, one may use the Bragg additive rule, w k 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), S , as the stopping power in the material, S, divided by the stopping power in water, S w , 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, 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 where S j is the relative stopping power at depth u j . 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.

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 MLP x SH hereon.

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 [13]. The validity of energy independent RStP is, however, not so straightforward. We explore this below.
RScP (T ) and RStP (S ) were calculated for the forty-six human tissues listed in ICRU Report 46 [29], 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. ( = 0.001225 e g/cm 3 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 [30]. The nominated RStP assigned to each material was arbitrarily defined as the average of the energy dependent RStP E ( ) calculated at 1 MeV incremeents over the range E 10 300 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 [31]. 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.

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 T and S for the solid tissues listed in ICRU report 46 [29]. 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.
A linear least-squares model is applied to the data. In this method we assume where p p , 1 2 are constant parameters, given in Table 1

Relative Scattering Power
Linear fit (solid tissues) Linear fit (low density materials) Fig. 1. Correlation between relative stopping and scattering powers (S and T respectively) for the materials listed in ICRU Report 46 [29]. 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 [29]). Parameters for the fits are detailed in Table 1.

Table 1
Parameters for the least squares fits, = + S p T p  constrained to zero for low density materials.

Monte Carlo simulations
The performance of the inhomogeneous MLP formalism, MLP x , 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 [20]. MLP x was tested against the homogeneous formalism, MLP H2O by Schulte et al. [9] 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 [32]. 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.
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 [9].
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 MLP x and MLP H2O as the overall likelihood of scattering is increased as a result of decreasing the beam energy. Implementation of MLP H2O for beam energies other than 200 MeV requires recalculation of the fifth order polynomial for the energy factor, p 2 2 , as detailed by Schulte et al. [9]. 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) [33,34] was imported into TOPAS as standard Digital Imaging and Communication in Medicine (DICOM) CT files using the internal GDCM package [20]. 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 MLP x and MLP x 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. [1]. The RScP was then calculated using (26), allowing the calculation of the inhomogeneous MLP. (Refer to Algorithm 1 and Algorithm 2 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 [9,35]. 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 edgedetection filter resulted in an average subsampling rate of 10% of the available MLP x SH calculation points. An example of this subsampling is shown in Fig. 3 (bottom right).
When comparing our inhomogeneous formalism to the homogeneous water approach, a hull-detection method was used prior to the calculation of MLP H2O , 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.

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 2 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.

Slab phantom
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 (MA-TLAB R2017b, The MathWorks, Inc., Natick, Massachusetts, United States) on a standard Intel Core i7 processor. The spline-hybrid method, MLP x 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, MLP x SH retains greater accuracy than MLP H2O . 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,  [33,34]. 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.) significantly more tracks escaped the MLP H2O 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 MLP x and MLP x SH give a larger and more right-skewed probability envelope surrounding the proton path when compared to MLP H2O for this particular track.

Human head phantom
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 MLP H2O (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 MLP H2O (Hull). Both MLP x and MLP H2O (Hull) recorded runtimes consistently below 24 ms per track, however MLP x SH was able to complete the task in half this time.

Discussion
A new inhomogeneous most likely path formalism using the matrix methods of Schulte et al. [9] has been proposed for use in iterative image reconstruction in proton computed tomography. Gottschalk's method [13] 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) and (25). We therefore recalculated MLP x (and MLP x SH) using a 5th order polynomial estimate of the mean energy, akin to the energy factor used in MLP H2O [9]. 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 Table 2 Percentage of proton tracks that escape the MLP 3 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.  inhomogeneous formalism in a 20 cm cube of water, MLP x gave the same path as MLP H2O but with a slightly larger probability envelope when comparing estimates to Monte Carlo data obtained using TOPAS [20]. Note that the implementation of MLP H2O requires calculation of the energy factor polynomial p 2 2 , 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 MLP x . In the inhomogeneous slab geometry consisting of large slabs of cortical and cranium bone, MLP x achieved greater accuracy than MLP H2O by up to 17% at the lowest beam energy of 210 MeV. It was shown that the improvement in accuracy gained by using MLP x over MLP H2O 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 [36].

MLP
x has in general outperformed MLP H2O 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 [13]. 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, f E E u ( , ( )) dM 0 , depending only on the initial energy of the proton, E 0 , and its energy at depth u. Namely,  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 [16], off which the inhomogeneous MLP proposed by Collins-Fekete et al. [15] 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 [37]. In a solid-state pCT system design by Esposito et al. [38], 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 [39]. Collins-Fekete et al. [15] 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 MLP x 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, [11]) 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 [31]. The reader is referred to the recent investigation by Khellaf et al. [40] on the effects of transverse tissue heterogeneities on the accuracy of the conventional water MLP formalism by Schulte et al [9].
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 MLP x and MLP H2O (Hull) performed similarly in terms of computational time, however MLP x 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.

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 [9] and the scattering power formulae of Gottschalk [13]. 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.