Statistical limitations in proton imaging

Proton imaging is a promising technology for proton radiotherapy as it can be used for: (1) direct sampling of the tissue stopping power, (2) input information for multi-modality RSP reconstruction, (3) gold-standard calibration against concurrent techniques, (4) tracking motion and (5) pre-treatment positioning. However, no end-to-end characterization of the image quality (signal-to-noise ratio and spatial resolution, blurring uncertainty) against the dose has been done. This work aims to establish a model relating these characteristics and to describe their relationship with proton energy and object size. The imaging noise originates from two processes: the Coulomb scattering with the nucleus, producing a path deviation, and the energy loss straggling with electrons. The noise is found to increases with thickness crossed and, independently, decreases with decreasing energy. The scattering noise is dominant around high-gradient edge whereas the straggling noise is maximal in homogeneous regions. Image quality metrics are found to behave oppositely against energy: lower energy minimizes both the noise and the spatial resolution, with the optimal energy choice depending on the application and location in the imaged object. In conclusion, the model presented will help define an optimal usage of proton imaging to reach the promised application of this technology and establish a fair comparison with other imaging techniques.


Introduction
Within the last two decades, proton and heavier hadron radiotherapy have experienced a surge of interest as a novel, potentially favourable, treatment option for complicated tumours near critical structures. This is mainly due to the advantageous dose deposition profile of protons known as the Bragg peak. The Bragg peak produces a dose peak highly focused at the end of the proton range, leaving a low entry dose and almost no dose in the distal region beyond it. On the other hand, proton therapy is limited by (1) a larger beam penumbra than photons for deep-seated tumours, (2) range mixing caused by complex heterogeneity and (3) scattered secondary radiation (Gottschalk 2012, Gottschalk et al 2014. Nevertheless, the potential for increased survival rates through dose escalation in the target volume, as well as the decreased toxicity due to the potential for reduced dose in the organ at risk have sparked the interest with many new centres starting treatments (Jermann 2015).
However, proton therapy's main advantage comes at the cost of limited robustness against uncertainties in treatment planning. Range uncertainties may push back, or project forward the Bragg-peak, leading to excessive dose delivered to the organ at risk and/or under-dosage of the target volume. Inter-fraction morphological changes, anatomical deformation due to intra-fraction motion, and treatment planning uncertainties (due to the empirical conversion of x-ray CT attenuation coefficients to relative stopping power) are the major causes of the range uncertainties (Yang et al 2012, Paganetti 2012. To improve robustness, non-optimal beam directions are often chosen to spare critical organs at risk, partially negating the advantages of proton therapy. Proton imaging has been proposed as a potential solution to these problems (Hanson et al 1981). Proton images are acquired by measuring the residual energy of a beam of protons exiting a patient. Higher entrance energy is required than for treatment to cross the patient. Previous studies have demonstrated that protons may have a significantly lower noise level and lower dose to the patient (Schulte et al 2005, Depauw and Seco 2011) than the conventional x-ray CT imaging. It could, therefore, be used more frequently for the same end-point dose, minimizing the risk of unaccounted for inter-fraction morphological changes. Furthermore, proton therapy imaging is acquired with the same beam that is used for treatment and samples directly the relative stopping power, minimizing uncertainties in the conversion of an imaged quantity to stopping power (Schneider et al 2004). However, protons suffer a significant amount of elastic scattering with nuclei through their trajectory in the form of multiple Coulomb scattering (MCS), reducing severely the spatial resolution of proton imaging (Schneider and Pedroni 1994). Advanced trajectory estimation methods have successfully helped address the problem of MCS in proton imaging, ameliorating the spatial resolution (Collins-Fekete et al 2017, Schulte et al 2008, Collins-Fekete et al 2015, Williams 2004. Thus, daily proton tomography (Bashkirov et al 2007, Bashkirov et al 2009, or the combination of daily proton radiography with conventional x-ray CT (Collins-Fekete et al 2016, Takatsu et al 2016, Schneider et al 2005, are currently investigated to minimize the effect of range uncertainties, patient positioning uncertainties (Schneider and Pedroni 1995, Poludniowski et al 2015, Depauw et al 2014 and breathing motion effects (Han et al 2011).
To understand and optimize the use of proton imaging in proton therapy, Schulte et al (2005) have investigated the density resolution achievable by this modality. They have shown that the noise-to-dose is expected to be lower or equivalent to x-ray CT at low to middle energy. They, however, neglected the components of scattering in their noise estimate and assumed that the full image noise is due to straggling. The MCS component of the noise was investigated by Rädler et al (2018) to help better fluence-field modulated proton CT. They have shown that the MCS noise may be very important at high-gradient region within the images, often overtaking the straggling noise contributions.
Although the noise characteristics are now better understood, no full characterization of their relation to spatial resolution, energy and energy loss has been done. This work perform an end-to-end characterisation of the proton imaging signal-to-noise (SNR) ratio and spatial resolution against delivered dose and entrance energy for the different existing proton interactions. By investigating these characteristics, a model of proton imaging's statistical limitations is defined and a relation between these factors and range uncertainties is established. The endpoint of this work is to propose a coherent framework that integrates the different noise sources and predict the optimal spatial resolution and noise level for a given beam and detector characteristics.

Theory and model
The purpose of this work is to model a computed tomography scan using a set of projections from a proton beam which is passed through an object and to understand the relationship between the delivered dose, the noise, the signal and the spatial resolution. The expected energy loss and noise is defined through the transport theory of protons. The noise is found for a projection and then extrapolated by back-projection to a tomograph. The signal-to-noise ratio is investigated in the middle of a uniform cylinder and related to dose and pixel size. Last, the spatial resolution degradation due to the uncertainty in the path estimate is investigated. Monte Carlo simulations are performed to validate the model and the simulation's setup is detailed at the end of this section.

The inverse problem
In proton medical imaging, the relevant reconstructed quantity is the ratio of the stopping power (RSP) of a material to the stopping power of water since it can be used to predict particle's range and it varies by less than 0.7% above 70 MeV (Arbor et al 2015). The local energy loss at a point r in space can be expressed as: where RSP(r): R 3 → R represents the local RSP, S w represents the proton stopping power in water given by the Bethe-Bloch equation (Bethe and Ashkin 1953) for a given mean excitation energy and local energy (E(r): R 3 → R). A proton will suffer a series of Coulomb interactions that deflect it from a linear path. Therefore, the inverse problem must be defined along the proton's path, namely Γ ∈ R 3 , rather than along a Cartesian axis. By integrating both sides of the expression, the inverse problem is given by: where the variable WET represents the water equivalent thickness of the projection which is the quantity required in proton imaging to reconstruct the tomograph. The exact proton path Γ is generally unknown but may be estimated by approaches such as a cubic spline fitted path (Collins-Fekete et al 2015) or a Bayesian likelihood maximization (Schulte et al 2008).

Geometry and definition
The problem will be approached on a 2D plane geometry neglecting the third dimension without loss of generality. An object is described by a function RSP(r) on this plane, where r is a two dimensional vector with component x and y. Furthermore, it will become useful to define a rotated orthogonal basis with axis t and s at an angle θ to the original system. This basis transformation will be used to define the projection at a fixed angle. The geometry is schematized in figure 1.The transformation between the (x, y) and (t, s) space, which is a rigid rotation of angle θ, is given by: The determinant of the Jacobian between the two bases is unity, thus dtds = dxdy.

Transport theory for proton imaging
In proton imaging, the measured signal is the energy loss (f E (t)) which is converted to WET (f WET (t)) subsequently (equation (2)). A definition of the noise in term of energy loss will help separate the contributions of the scattering and straggling processes. The signal will be expressed in energy in the first part of this work. For a beam of particle, the measured signal and noise vary with the beam scattering distribution, which is a function of 1) the detector type and the number of detection point (Krah et al 2018) and 2) the particle charge/mass ratio and velocity (or initial energy) (Collins-Fekete et al 2017). The scattering distribution is illustrated in figure 2-black curve for a scenario with two detection points with no detection uncertainty. The measured signal in an exit pixel is the energy loss of particles that have followed that scattering distribution, interacted with the object it passes through, and were measured in that exit pixel. Schematic representation of the scattering function of a beam of protons crossing a phantom with two trackers at the entrance and exit of the object. The lateral scale is exaggerated for representation purposes. The blue shape represents the imaged object. The area delimited by the black lines represent one standard deviation of the scattering probability distribution. The distribution decreases to zero at both measurements point on the trackers. The scattering function is calculated at each depth, and a single realization is shown in the zoom box. The red line represents an iso-probable path. The probability bin is represented by the transparent red area in the zoom-box and the red dot represents the mean position from which the iso-probable path is extracted.
In this work, we introduce a method to rapidly calculate the convolution of the particles scattering and energy loss with the underlying object. It is important to first note that the scattering distribution is approximated as a Gaussian normal distribution for which the standard deviation is calculated by the Fermi-Eyges equation (Eyges 1948) at a specific depth. When normalized, the distribution represents the probability of finding a particle as a function of lateral deviation.
To estimate the measured signal in an exit pixel, we integrate the scattering distribution over the probability. For each depth in the object, we first construct a histogram which is the probability density function of lateral displacement at that depth. The probabilities is discretized in 200 bins covering a range defined by the ± 5σ region. Then, for each bin in the numerical integral representing a discrete probability interval (red transparent box in figure 2), the procedure goes as follow: for a specific depth, the lateral scattering distribution is calculated and normalized to represent probabilities. The lateral deviation for this probability bin is found from the Gaussian distribution. Then, a trajectory is constructed by pairing the depth position and the lateral position (red dot in figure 2) and repeating the process for each depth. The constructed trajectories are named iso-probable paths. In this work, 1 mm depth interval are used to calculate iso-probable paths. An example of an iso-probable path is represented in figure 2 by the red line. The iso-probable path is defined as Γ iso (p, Y 0 , Y 2 ). The expected output signal is then: In (equation (4)), p represents a fixed probability for which the iso-probable path Γ iso is derived. The lateral scattering distribution at each depth is calculated with the Fermi-Eyges theory of charged particle scattering (Eyges 1948). ∆E(Γ iso (p, Y 0 , Y 2 )) represents the energy loss along the path and is calculated with the continuous slowing-down approximation (CSDA). A feature of interest is that the expected energy loss measured is dependent upon the path uncertainty, which itself depends on the tracker precision and the properties of the particle (i.e. charge to mass ratio and energy). For example, higher uncertainty on the trackers may lead to a larger scattering function (Krah et al 2018). In an heterogeneous object, this larger uncertainty would yield a different expected energy loss compared to a less uncertain tracker.
The model is validated for this work on proton radiographs binned: at the entrance trackers, exit tracker and reconstructed using both tracker. The expected energy loss for the entrance and exit tracker reconstruction method are calculated using (equation (4)) for every single tracker position. The radiograph with both trackers is produced by calculating the probability integral over every possible entrance/exit position combination.
Collins-Fekete et al (2016) approach will be used here for tomography. They propose to use a normal back-projection reconstruction algorithm on the deblurred radiographs. As mentioned earlier, the back-projected quantity is the WET which is found from the energy loss through (equation (2)). A back-projection operator is defined as: where ∆θ = 2π/M is the increment in projection angle and M is the number of projections. (Equation (5)) represents a discrete angular integration. The range of the integral is 2π, because of the asymmetry of the scattering function of a proton with energy loss. From the theory of tomographic reconstruction, the back-projection of a line integral convolved with a suitable filtering function (h(t)) should return an estimate of the original quantity:R The difference between the original quantity (RSP(x, y)) and the measured one (R SP(x, y)) comes from the blurring of the signal distorted by the multiple Coulomb scattering of the particles. Importantly, this means the blurring alters the reconstructed signal in the form of a systematic shift and thus proton imaging is not a perfect quantity preserving tomographic reconstruction method.

Noise in proton radiography
The variance in energy loss due to the path's uncertainty is defined as The expectation operator is defined as in (equation (4)). A perfectly-known entrance energy is assumed (σ E0 = 0). Due to the linearity of the expectation operator, the variance in exit energy is equivalent to the variance in energy loss, i.e. σ 2 ∆E = σ 2 Eout . This result represents the noise inter-path. For two protons following an almost identical trajectory through a medium, the amount of energy loss is subject to two other sources of fluctuations. The number of proton-electron collisions and the energy loss in each collision can fluctuate statistically. The energy variation introduced by these fluctuations is called energy straggling. Energy straggling has been described thoroughly by Vavilov (1957) and more recently approximated by Tschälar (Tschalar 1968, Tschalar andMaccabee 1970). The energy straggling standard deviation is defined by the Tschalar equation (Tschalar and Maccabee 1970): where k 1 and k 2 are defined as: In this equation, c is the speed of light, β is the proton velocity relative to the speed of light, η e is the relative electron density of the medium to the electron density of water, m e is the relativistic electron rest mass, I(s, t) is mean excitation energy of the medium, and the constant K = 170 MeV cm −1 combines various fixed physical parameters. Straggling is calculated as a function of the exit energy. For the distribution of exit energies observed by an exit detector position, the straggling must be calculated and weighted for each path probability: The Coulomb scattering and straggling cross sections both depend on energy and energy loss and hence are correlated. In a first approximation, we assume that the covariance is null. The impacts of this assumption are discussed in section 4. For our assumption, the combined error is the root of the sum of the squared individual standard deviation (σ 2 Eout = σ 2 E,strag + σ 2 E,MCS ). The reconstructed quantity used in proton imaging reconstruction is the WET and the noise should also be expressed in these units. The WET is related to the energy of the exiting particle by (equation (2)), and the error is calculated through the error propagation formula. If higher-order terms are neglected, the WET variance can be expressed as: The noise is expressed as standard error on the mean, i.e. divided by the number of detected protons (N D ).

Noise in proton tomography
The noise in proton tomography has been recently described by Rädler et al (2018). Briefly, if one neglects interpolation effects required when a limited number of projections is acquired, the variance can be expressed as the back-projection of the variance convolved with the square filter:

Signal to noise ratio
In this section, the SNR is determined in the centre of a uniform radially symmetric object. We assume that the number of protons detected in the centre bin is constant over all projections as the object is uniform and cylindric, i.e. N D,θ (t) = N D . First, let us define the SNR as: The first estimate of the object RSP (R SP(t, s)) is obtained by back-projecting the radiographs using the operator presented in (equation (5)). The noise in the middle of the phantom can be calculated by acknowledging that the noise in the middle should be similar for all projections since the object is radially symmetric: where the sum of the filter's intensity, when using the Ram-Lak convolution filter (Ramachandran and Lakshminarayanan 1971), is equal to 1/12a 2 where a 2 is the pixel size (Gore and Tofts 1978). The SNR equation in the middle of the phantom is:

Dose in the middle of a cylindrical object
The dose in the middle of a uniform water cylinder of diameter d has been derived by Schulte et al (2005). However, for this study, the results are calculated through Monte Carlo simulation. This allows us to extend the dose results to heavier particles for which a defined model does not exist that accounts for all secondaries. Precisely, the average dose deposited per particle in a voxel at a position r is calculated and denoted by S MC (E, r). The equation for the dose in the middle of a phantom is: where E c is the energy of the proton in the middle of the phantom and Φ MC (E c , d/2) represent the flux of particles at the centre of the phantom for a single projection and all others variables are defined as before.
The flux halfway through the cylindrical object due to nuclear reactions can be described from the measured number of protons detected as: where N D is the number of detected protons at the exit detector and a is the pixel size. The flux loss through nuclear reactions is denoted by the variable g MC (σ nuc ) and is calculated through Monte Carlo simulations, by extracting the fraction of particles that crosses the middle voxel but suffer nuclear annihilation. This leads to a final expression for the dose in the centre of the phantom as a function of the detected number of protons (N D ): If the number of protons detected (N D ) is substituted into (equation (15)), the final relationship between the dose, the SNR, the energy and the pixel size (a) can be defined as: This equation is valid for any ion, given that the average dose deposited in the voxel per particle (S MC (E c , d/2)) and the attenuation (g MC (σ nuc )) are replaced by quantities specific for that ion species. If one considers only the electromagnetic interaction component of deposited dose, then (equation (19)) can be rewritten as: where the impact of the blurring of the measured quantity on the SNR is outlined by the ratio between each.

Spatial resolution: pixel size and scattering effect
Spatial resolution is an important imaging metric for treatment planning as well as diagnostics. As in other imaging techniques, the spatial resolution of proton imaging is affected by the detector, the reconstruction algorithm, as well as the sampling resolution of the signal. However, in addition to those effects, the spatial resolution in proton imaging is affected by the non-linear trajectory followed by the protons, which is caused by MCS. For the scope of this work, only the sampling resolution and MCS effect on the modulation transfer function (MTF) are investigated, both separately and in combination. Firstly, the MCS effect leads to trajectory estimation uncertainty. The process can be understood by considering the imaging of a Dirac function object. Due to the uncertainty region displayed in figure 2, the reconstructed image must be the convolution of the Dirac function with the scattering function width, evaluated at all depths in the object. The scattering function, following the Fermi-Eyges scattering theory, is a normal distribution with a standard deviation defined by σ scatt . The MCS component of the spatial resolution can be derived by taking the Fourier transform of the convolution of the Dirac function and the scattering. It will be represented by a centred Gaussian function with standard deviation σ MTF = 1 2πσscatt , as expressed by Plautz et al (2016). The spatial resolution will furthermore be limited by the pixel size denoted a. The pixel size effect may be represented by a rect function in the spatial domain and a sinc function in the frequency domain: The sampling will be represented by a multiplication of the rect function with a Dirac comb (III a ). Let us define the MCS MTF component as the function g(x) in the spatial domain and G(ε) in the frequency domain as calculated at the position r. The output MTF will be the result of the convolution between the Gaussian spread due to the MCS, the pixel size, and the sampling frequency: Except at the edge of the phantom, the blurring of the signal due to MCS will degrade the spatial frequency more rapidly than the sampling frequency (pixel size). Thus, the effect of the comb function is neglected.
The size of the scattering function spread will vary with depth and the optimal pixel size may depend on the reconstructed region within the object. The pixel size will be investigated where the scattering is largest in the phantom, therefore imposing a minimal spatial resolution correlated with minimal noise.

Geant4 MC simulations
Monte Carlo simulations were carried out to produce projection data to validate the model described in this paper. MC simulations in this work were implemented using Geant4 MC code version 10.2.1 (Agostinelli et al 2003).

Physics package
In this work, parameterized interactions with nuclei and elastic/inelastic ion interactions are considered exclusively for the extra dose involved, but are tagged and removed when evaluating the electromagnetic noise. The model aims to represent electromagnetic interactions only and the introduction of nuclear interactions would introduce unnecessary uncertainties against the goal of the model. Furthermore, it is expected that nuclear interactions can be filtered out of the signal using the recent dE-E filter developments proposed by Volz et al (2018) for proton imaging, and would affect only the noise. The processes considered include electromagnetic energy loss and straggling (following Bethe-Bloch theory) and MCS based on Lewis theory (Goudsmit and Saunderson 1940) using the Urban model (Urban 2006) as well as elastic/inelastic nuclear interactions. In precise terms, for all particles the following physics lists were used: the standard electromagnetic option 3 for high accuracy of electron and ion tracking and the ions elastic model (G4HadronElasticPhysics). For inelastic interactions, the binary cascade models (G4IonBinaryCascadePhysics) is used for protons and the quantum-molecular-dynamics (QMD) model for heavier ions (G4IonQMDPhysics). Radioactive decays module is used for all ions (G4RadioactiveDecayPhysics). These physics lists were chosen based on recommendations from Lechner et al (2010). The energy straggling and MCS processes were inactivated in some simulations to evaluate the noise coming from each component as well as the cumulative noise.
Step limiter cuts were set to 1 mm.

Beam setup
Protons (n = 10 7 protons/simulation) were simulated through water cylinders of 10 and 15 cm radius with initial energy varying between 110 and 300 MeV in steps of 10 MeV. The initial beam flux was distributed evenly along the lateral side of the simulation world, centred on the water cylinder. No initial angular deviation was given to the protons.

Acquisition setup
The simulation world is defined as either a 30 cm 3 (40 cm 3 ) air box in which a 10 cm (15 cm) radius water cylinder sits in the middle, with height matching the world size. Projections are binned at the rear tracker. Bin size of 1 mm was chosen for every reconstruction. Protons are recorded at the exit planar detector located at the distal edge of the simulation world. In both scenarios, this yields a minimal air gap of 5 cm between the frontal/distal detectors on either side of the phantom. A schema of the acquisition setup is shown on figure 3.

Results
In this section, the proposed framework is first validated against Monte Carlo simulations. From there onward, the focus is put on analytically produced results. The implications of the previously developed model is then investigated. The parameters that are studied in the noise profiles are the proton initial energy as well as the diameter of the water cylinder. Noise profiles are separated into noise due to the scattering and straggling physical processes, as well as the total combined noise. Spatial resolution is investigated by evaluating the relative effect on the pixel size and the proton MCS spread as a function of the entrance energy and cylinder size. (Equation (19)) shows the relationship between the SNR, the dose, the exit energy as well as the pixel size. The effect of the pixel size in this model is investigated by producing a nomogram of the model with a fixed entrance energy. The effect of the entrance energy is also investigated by fixing the pixel size.

Noise profile in radiography
The   The noise profile on the rear tracker displays a significantly increased noise at the edge, due to the MCS contributions. Apart from the edge of the phantom where the Fermi-Eyges model of the scattering is expected to fail, the maximal difference is found to be below 0.1 cm WET. This allows us to conclude that the model predicts accurately the expected exit energy and the noise. The noise profile in the front tracker is significantly different from the noise profile seen in the rear tracker. The noise is considerably lower, with minor contributions from the MCS noise seen on the two wings around the central position on figure 4. The middle bin noise level is equivalent to that of the rear tracker, but lower everywhere else. This is seen on both energies for both phantom radius. The noise profile for both trackers is similar to the noise seen on the front tracker only, although surprisingly slightly higher. (figure 4), with a much lower impact from the scattering and the noise being mostly dominated by the straggling. Figure 5. Noise profile as a function of the position on the lateral axis for proton radiography reconstructed at the rear tracker(left) and the front tracker (right). The noise is due to the different electronic interaction processes experienced by protons while crossing the cylinder. Shown here are the noise due to proton straggling, the noise due to the multiple Coulomb scattering and the combined total noise. The noise is shown for a water cylinder of 15 cm radius with proton initial energy of 300 MeV.

Noise due to individual physical processes
The radiographic noise profiles are shown in figure 5 as a function of the lateral axis. They are separated into three components: (1) the noise due to the straggling effect, (2) the noise due to the scattering effect, and (3) the combined total noise. The radiographic noise profiles are shown for a 300 MeV beam crossing respectively a 15 cm radius water cylinder.
A notable feature is a drastic increase of scattering noise closer to the edge of the cylindrical object for the rear tracker noise (figure 5). Similar features appear for the front tracker noise but are negligible against the straggling noise. Furthermore, the scattering noise is minimal in the middle whereas the straggling noise is maximal at the same position. In both cylinders, the straggling noise in the middle is very similar in size to the total noise, being the major component of it. This explains the accuracy of the results seen by Schulte et al (2005).

Noise profile in tomography
The next step is to study the variation in the noise as a function of the entrance energy at different radial positions in tomographic reconstruction. To do so, the radiographic analytical model of the noise is back-projected using a square filter (B{σ 2 θ,WET (t) ⊛ h 2 (t)}) for different entrance energies. Results are shown in figure 6 for the 10 and 15 cm radius cylinder respectively.
The noise profile for both the front tracker (not shown here) and the dual trackers methods (dotted line on figure 6) are highly similar, decreasing with decreasing energy.
Recent studies in particle imaging (Schulte et al 2005) indicated that a lower entrance energy correlates with lower noise, due to an increase in the signal caused by the higher energy loss in the cylinder which manifests itself in the energy loss-WET conversion (equation (2)). Our results corroborate those findings in the middle of the cylinder for all three methods evaluated here. In the case of front tracker binning and dual trackers, this conclusion holds true away from the center. In the case of rear tracker binning, the noise increases drastically away from the center of the cylinder, due to the scattering component becoming dominant.

Spatial resolution
The combined effects on the transfer function of the pixel size and the scattering function were evaluated. To fix the effect of the particle energy loss and the width of the scattering function, the pixel size is expressed in term of the standard deviation of the Gaussian scattering function at the center of the phantom. The transfer function is going to be dominated by the Gaussian scattering. The left panel of figure 7 shows the apodising effect of the pixel size on the Gaussian scattering component. The right panel shows the scattering components of the MTF as a function of the entrance energy for both cylinders. The results shown in figure 7 are produced by solving the σ scatt equation for two measurement points as detailed in section 2.3 in the middle position of the phantom. Figure 6. Tomographic noise as a function of beam entrance energy for a beam of proton crossing a 15 cm radius water cylinder binned at the rear (right) and using both trackers (left). The binning using the entrance detector is similar to the latter and omitted here for clarity. The noise is shown at different radial distances from the middle of the cylinder to demonstrate the effect of MCS noise as a function of the position within the cylinder. In light of the results shown in figure 7, it seems that a pixel size of one standard deviation of the scattering width will not severely affect the transfer function (∆ MTF 10% ≈ 4%). These results can be anchored to absolute spatial resolution by investigating the Gaussian scattering transfer function component at different energies. This is shown in figure 7 where the MCS component of the MTF 10% (independent of the pixel size) is shown for protons crossing a 10 or 15 cm radius water cylinder, as a function of the entrance energy. Since the MCS component is the major MTF limiting factor, particularly at the middle of the phantom, it is possible to approximate the spatial resolution of a reconstructed image. A tomographic reconstruction of a 200 MeV proton beam crossing a 10 cm radius water cylinder, choosing a pixel size of one SD (≈ 0.5 mm), would yield a MTF 10% 4% below the 5.6 lp/cm (right panel of figure 7) level in an otherwise perfect setup. This result is lower than the result found by Plautz et al (2016) for a 197 mm thickness phantom (6.3 lp/cm). It is slightly higher than the results produced in the recent work by Collins-Fekete et al (2016) (4.53 lp/cm) although their setup has a larger water equivalent thickness than the one studied here and they use a reconstruction algorithm that may reduce the spatial resolution. Similarly, Li et al (2006)  produces reconstructed images using the MLP around 5 lp/cm, again using a tomographic reconstruction technique that may lower the theoretical achievable spatial resolution. Those results represent only the achievable spatial resolution without any other effects such as detector resolution, reconstruction algorithm precision, contamination of the detected signal with large-angle scattering events and energy resolution of the detector. In a real scenario, a lower spatial resolution is expected.

Dose nomograms
Finally, the impacts of the pixel size and entrance energy on the Dose/SNR relationship are studied for a central pixel. To do so, nomograms are produced by fixing one of these parameters. The nomogram for a fixed energy and a fixed pixel size are shown in figure 8. Firstly, for a fixed dose, decreasing the pixel size has a high cost on the SNR with a minor impact on the spatial resolution (figure 7). This is due to the inverse fourth-power relationship with the pixel size seen in (equation (19)) and well known in x-ray imaging.
The energy nomogram (figure 8) shows results in-line with what has been found by Schulte et al (2005). This comes to no surprise as the difference with this framework and their model is the inclusion of MCS, which is minimal at the centre of the phantom. Alongside the proton nomograms is shown a mono-energetic 75 keV photon beam nomogram, as derived from the Barrett (Barrett et al 1976) model. The total noise in the middle of a phantom monotonically decreases with energy, consequently increasing the SNR. However, it is important to point out that this property will not be preserved at a different radial position and the nomograms will differ, following the relationship shown in figure 6.

Noise and spatial resolution limitations as a function of energy
Finally, given the definition of the noise, the spatial resolution, and their relation with dose through the model developed in this work, we combine these results in figure 9 to present the full picture of spatial resolution and noise as a function of the energy. The anti-correlated nature of spatial resolution and noise against energy is clearly outlined. In all cases, as the energy increases, the SNR decreases and the MTF increases. A lower proton beam energy (lower spatial resolution) will introduce range mixing (Espana and Paganetti 2011) and increase the blurring of the measured signal (Sawakuchi et al 2008) due to a larger scattering distribution, which will impact and distort the Bragg peak, thus affecting the range uncertainty. A higher beam energy will increase the noise level, which will impact the range margins required to ensure adequate target coverage in proton therapy.

Discussion
In this paper, we have presented a framework for modelling the effects of the statistical characteristics of the imaging process on the image quality in proton beam CT. The approach is informed by the model for Dose = 0.5 mGy Dose = 1.5 mGy Dose = 2.5 mGy Dose = 3.5 mGy Dose = 4.5 mGy Figure 9. The figure demonstrate the impact of beam energy against the SNR(noise) and spatial resolution for a 200 MeV beam crossing a 20 cm cylinder phantom. On the left-axis, the SNR is plotted against dose. The spatial resolution (red) is shown on the right-axis. The figures demonstrates that the inverse relationship of SNR and spatial resolution against energy. An higher energy leads to a higher spatial resolution and an higher noise. Pixel size is defined as the standard deviation of the scattering distribution as defined earlier (section 2.8). describing x-ray CT by Barrett et al (1976). In Barrett's work, they describe these effects for x-rays in term of the relationship between dose, spatial resolution, signal and noise (with the latter described as the SNR). Gore and Tofts (1978) explored the statistical limitations in x-ray CT and arrived at an expression for the standard deviation of the attenuation coefficient as a function of imaging parameters including pixel size and the number of incident particles. X-ray CT involves a binary process where x-rays either pass through the object unhindered or are absorbed or scattered and do not contribute to the imaged signal. In the case of protons, the particles used for image reconstruction are assumed to have passed through the object but are subject to multiple scattering. Although some particles suffer nuclear elastic and inelastic reactions, they are rejected using recent advanced filtering techniques . The scattering leads to spatial uncertainty in the path of the proton and hence a source of spatial uncertainty in imaging not seen in x-ray CT. Schulte et al (2005) built on the work of Gore and Tofts and arrived at a model describing density resolution in proton CT. Their work did not include consideration of multiple Coulomb scattering (MCS), which forms a significant effect. Recently, Rädler et al (2018) have characterized the MCS component of the noise and described how the total noise affects the reconstruction of the fluence-modulated pCT algorithm. Their work, however, did not include any spatial resolution effects in their model nor analyzed the degradation of the expected signal as a function of the scattering. One of the purposes of our work is to include MCS in a novel manner that incorporates the systematic shift and stochastic noise as well as the effects of the spatial resolution on the noise.
The framework developed has been used to generate noise estimates for measurement of radiographic projections of uniform cylindrical objects on the rear and front tracker and reconstructed using both trackers. These noise estimates are in the context of energy straggling noise and noise due to MCS. As of now, the different sources of noise in the imaging process are considered independent. The model could be developed by investigating the covariance between the straggling noise and MCS noise, however, the accuracy of the results (figure 4) seem to indicate that the correlation is low when propagating the straggling noise through the scattering distribution (equation (10)). The effects over M projections with a back-projection filter have been modelled to generate data relevant to both proton radiography and CT. A noise model was developed by Rädler et al (2018) to incorporate in a fluence modulated proton projection and CT. However, their model did not seem to account fully for the variation close to the cylinder edge as an underestimation of the noise is seen compared to MC noise away from the centre pixel.
The projection data are shown here for rear-binned trackers, front-binned trackers, and dual-trackers reconstruction. A tomograph is produced using dual-trackers radiographs as input in a classical back-projection x-ray algorithm (Collins-Fekete et al 2016). The inclusion of a MLP path algorithm in the tomographic reconstruction (Penfold et al 2010, Rit et al 2014 will reduce the scattering function (figure 2). However, since the scattering has a minimal impact on the total noise (figure 6), including advanced reconstruction technique should not alter significantly our results. An interesting feature shown here is the slightly lower noise in front tracker radiography compared to dual-trackers reconstructed radiograph. This potentially indicates that for a low noise proton tomograph reconstruction, a front tracker may be sufficient at the cost of a lower spatial resolution.
Equation 20 details the relationship between dose, SNR, signal degradation and pixel size, for a central pixel. Then, (equation (22)) describes the relationship between pixel size and spatial resolution. An important difference between proton and x-ray CT is that with x-ray CT, the spatial resolution can, in principle, be made arbitrarily high, at the expense of high dose, by having a very small pixel size. In the case of proton CT, the spatial resolution is limited by the finite uncertainty in particle trajectory caused by MCS, effectively producing a limiting spatial resolution. Furthermore, the scattering suffered by the particles induces a blurring in the reconstructed quantity which reduces the SNR for a fixed dose (equation (20)). In this work, we have suggested a limiting pixel size of one SD of the Gaussian spread of the MCS distribution (figure 7), which may be useful as a lower limit for future detectors development, when combined with detector position/sampling effects on the MLP/MTF described by Krah et al (2018).
In figure 8 we present nomograms relating dose and SNR as a function of both pixel size (expressed as a function of the SD of the MCS spread) and incident proton energy. Again, the general picture for proton CT is more complex than for x-ray with energy being a variable that affects the noise 1) directly through WET conversion effects and 2) indirectly through intrinsic spatial resolution/pixel size effects and thus number of particles contributing to noise in this pixel. This framework shows that there is a complex trade-off between the noise effects of energy straggling and those of MCS ( figure 9). This trade-off changes with incident energy, object size and location in the object being imaged (see figure 6). As a consequence, the optimum imaging parameters will depend on object size and the location of the area of interest in the object. Hence for a given patient diameter, the optimum imaging parameters will be different for a tumour region located in the centre and for one closer to the surface. Nevertheless, it is, as of now, difficult to define recommendations for parameters in proton imaging given that the impact of the imaging quality metrics on range uncertainties is not yet fully understood.
A major difference between this model and previous work (Schulte et al 2005) is the introduction of the SNR metric which includes the expected reconstructed signal. As shown in our framework, the expected exit energy is a superposition of the energy loss of protons following the scattering distribution across the objects (equation (4)) which results in a blurring of the measured quantity. This blurring alters the reconstructed signal, especially in the case of small inserts, a feature observed by Piersimoni et al (2018). The signal-to-noise ratio is reduced by 1) the noise in the image due to straggling and scattering, and 2) by the blurring of the measured quantity due to scattering.
The noise and spatial resolution characteristics presented in this work will be important when considering the use of proton imaging for its main purpose, i.e. radiotherapy and range uncertainties. The presence of noise is expected to increase the stochastic errors in the measured quantities requiring a wider margin to ensure coverage of the target. The noise should therefore be minimised, through the use of a lower energy for imaging. However, a lower energy is correlated with a higher blurring and a lower spatial resolution. The lower spatial resolution may induce partial volume effects that produce range mixing (Paganetti 2012) (figure 9). The blurred signal may induce systematic error, translating into a systematic shift of the expected proton range. It is therefore not clear which energy would benefit proton radiotherapy treatment the most. Nevertheless, the proposed framework may be useful to predict the limits of proton CT as a treatment planning imaging tool as a function of the imaging parameters but those limits will be highly application-dependent.
Furthermore, we have corroborated results found earlier by Rädler et al (2018) that the noise is minimal in the middle part of the cylindrical phantom ( figure 5). This is caused by the minimal variation in WET seen by protons crossing this section of the cylinder. However, a human body is hardly a perfect cylinder of water and the conclusions drawn here must be nuanced in the application to patient imaging. Human body regions that are highly heterogeneous such as head and neck and lung will have higher scattering variations and therefore higher total noise than what is seen here.
There is a range of options for CT reconstruction in proton tomography. Firstly, a range of back-projection filter types is available for tomography. In this work we have assumed the filter has the characteristics as described by others ( Barrett et al (1976), Gore and Tofts (1978), Chesler et al (1977)) leading to the key step that we may separate the noise component for radially symmetrical objects in the middle (equation (14)). All filters in CT are by nature closely related to a ramp filter and hence will have those characteristics as discussed by Chesler. Secondly there are a range of methods of determining the path of a particle through the object being imaged, including most likely path approaches (Schneider and Pedroni (1994), Williams (2004), Schulte et al (2008), Collins-Fekete et al (2015, 2017, Erdelyi (2009), Wang et al (2010)). In essence these approaches involve a set of fixed, known points where sensors are placed to measure each particle's position (and possibly energy) and a model to describe the path between these known points, with associated spatial uncertainty. This framework uses the uncertainty in the path between the known points to describe the spatial uncertainty caused by MCS. For different reconstruction models, this translates to different spreads in the scattering distribution.
This work has made no assumptions about the characteristics of the detection system. There are a range of proton CT systems in development (Scaringella et al 2013, Poludniowski et al 2014, Johnson 2018, Bashkirov et al 2009, Civinini et al 2013 with a range of strategies. However, most combine position sensitive sensors to measure the path into and out of the object, coupled with a range telescope or calorimeter to measure the particle energy. Our presentation assumes this general approach. Detector response of the imaging system could easily be incorporated into this framework. The choice of different arrangements of tracking and energy detectors could also be incorporated. For instance, to have a single spatial tracker at the exit of the object would produce zero spatial uncertainty at this point and a Gaussian spread of increasing width towards the entrance of the object. This has been investigated recently by Krah et al (2018).
A final point related to imaging systems for proton CT reported in the literature is that some systems modulate the energy across the object to get sufficient energy to cross the object, reducing the need for thick and costly residual energy/range detectors. Our results suggest that, in this scenario, the detector precision becomes paramount and highly related to the image quality. Whereas a lower energy drastically increased the noise for a rear tracker binning, the same scenario resulted in a lower noise for a tomographic reconstruction including the MLP. In the case of a dual tracking detector system, the minimal energy that crosses the phantom while keeping and exit energy in the RSP range of validity will minimize the noise, at the cost of spatial resolution. Whether this would lead to an overall lower range uncertainty is not yet validated.
The findings of this work should be considered when using proton CT as a gold-standard for imaging for proton radiotherapy. In most situations proton CT will have poorer spatial resolution than other imaging methods such as x-ray CT. Also, the noise characteristics and expected signal will change across the field of view and so the choice of a gold standard may be problematic with a combination of imaging modalities and/or the use of priors likely to be a good strategy. Future applications of this framework include optimisation of beam energy and profile for proton radiography and CT of an object of known size and a known location of the area of interest, quantitative comparison and optimisation of imaging with different charged hadron species and modelling the combination of hadron and x-ray projections in CT for priors-based tomography in hadron radiotherapy.

Conclusion
A computed tomography scan using a set of projections from a proton beam which is passed through an object has been modelled and the relation between the delivered dose, the noise, the signal and the spatial resolution has been summarized. The optimal energy to minimize noise is the minimal energy that crosses the thickest section of the phantom while yielding an exit energy above the 70 eV threshold. Using this minimal energy will come at the cost of spatial resolution. Energy modulation as a function of the phantom thickness decreases the noise but decrease the spatial resolution, with the consequences of that choices not fully understood.