A maximum likelihood method for high resolution proton radiography/proton CT Charles-Antoine Collins-Fekete1,2,3 , S´ ebastien Brousmiche4 , Stephen K. N. Portillo5 , Luc Beaulieu1,2 , Joao Seco3,6,7 1

D´epartement de physique, de g´enie physique et d’optique et Centre de recherche sur le cancer, Universit´e Laval, Qu´ebec, Canada 2 D´epartement de radio-oncologie et CRCHU de Qu´ebec, CHU de Qu´ebec, QC, Canada 3 Department of Radiation Oncology, Francis H. Burr Proton Therapy Center Massachusetts General Hospital (MGH), Boston, MA, USA 4 Ion Beams Application - IBA, Louvain-la-Neuve, Belgique 5 Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, USA 6 Deutsches Krebsforschungszentrum Heidelberg, Baden-W¨ urttemberg, DE 7 University of Heidelberg, Department of Physics and Astronomy Heidelberg, Baden-W¨ urttemberg, DE E-mail: [email protected] Abstract. Multiple Coulomb scattering (MCS) is the largest contributor to blurring in proton imaging. In this work, we developed a maximum likelihood least squares estimator that improves proton radiography’s spatial resolution. The water equivalent thickness (WET) through projections defined from the source to the detector pixels were estimated such that they maximizes the likelihood of the energy loss of every proton crossing the volume. The length spent in each projection was calculated through the optimized cubic spline path estimate. The proton radiographies were produced using Geant4 simulations. Three phantoms were studied here: a slanted cube in a tank of water to measure 2-D spatial resolution, a voxelized head phantom for clinical performance evaluation as well as a parametric Catphan phantom (CTP528) for 3-D spatial resolution. Two proton beam configurations were used: a parallel and a conical beam. Proton beams of 200 and 330 MeV were simulated to acquire the radiography. Spatial resolution is increased from 2.44 lp/cm to 4.53 lp/cm in the 200 MeV beam and from 3.49 lp/cm to 5.76 lp/cm in the 330 MeV beam. Beam configurations do not affect the reconstructed spatial resolution as investigated between a radiography acquired with the parallel (3 .49 lp/cm to 5.76 lp/cm) or conical beam (from 3.49 lp/cm to 5.56 lp/cm). The improved images were then used as input in a photon

Maximum likelihood proton CT reconstruction

2

tomography algorithm. The proton CT reconstruction of the Catphan phantom shows high spatial resolution (from 2.79 to 5.55 lp/cm for the parallel beam and from 3.03 to 5.15 lp/cm for the conical beam) and the reconstruction of the head phantom, although qualitative, shows high contrast in the gradient region. The proposed formulation of the optimization demonstrates serious potential to increase the spatial resolution (up by 65%) in proton radiography and greatly accelerate proton computed tomography reconstruction.

Keywords:

proton radiography, proton computed tomography, Monte Carlo

Submitted to: Phys. Med. Biol.

1. Introduction The idea of proton computed tomography (pCT) was first theorized when Cormack (1963) proposed a novel tomographic method to estimate the density of matter through the energy loss of charged particles traversing a material of constant atomic composition. The experimental demonstration was made in a later publication (Cormack and Koehler (1976)). However, the use of non-deflected particles such as X-rays soon proved to be much more efficient with a reduced cost and the research in pCT halted. In the last decade, protons were brought forth as a therapeutic option with great potential for tumor treatment and control. The unique advantage of proton therapy over conventional photon-based therapy is the sharp dose peak, namely the Bragg peak, at a precise position along the beam path, with very low exit dose. With adequate planning, the delivered dose can theoretically be confined within the tumor zone. To fully utilize the advantages of proton therapy, a precise calculation of the Bragg peak distal position is required. This is done via the Bethe-Bloch equation (Berger (1993)) with ICRU corrections (Seltzer and Berger (1982)). This equation calculates the range of charged particles in a material using its atomic composition (Z), mean excitation energy (I-value) and density (ρ). The range calculation precision relies on accurate knowledge of those quantities. Yet, they are difficult to acquire without considerable uncertainties (Hanson et al. (1981)). To minimize the effect of these uncertainties, the concept of relative stopping power (RSP) was introduced by

Maximum likelihood proton CT reconstruction

3

Hanson et al. (1981). The relative stopping power is the stopping power of a material divided by the stopping power of water. It is equivalently the ratio of the water equivalent thickness (WET) to the material thickness. The RSP is a combination of the previous information (Zef f , ρ and I-value) in one quantity that can be measured precisely. It is approximately constant with energy, up to a variation of 0.7% in the 80-300 MeV range (Arbor et al. (2015)). To acquire RSP maps for the proton range calculation, the actual clinical method converts the X-ray CT Hounsfield Unit (HU) to RSP through stoichiometric calibration, a technique that separates the Z dependence of the material from its ρ dependence (Schneider et al. (1996)). However,this method has been shown to introduce error on the converted RSP of 0.8% on average (Matsufuji et al. (1998); Chvetsov and Paige (2010)) and 3% if the I-value uncertainty is accounted for (Schaffner and Pedroni (1998); Yang et al. (2012)). Zygmanski et al. (2000) demonstrated that pCT should be used in proton therapy treatment because it reconstructs the stopping power of the materials with better accuracy than the conventional conversion of HU to stopping power. In addition, proton imaging has a higher density resolution, 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. Finally, pCT suffers from different artifacts than X-ray CT and has diagnostic qualities (Depauw and Seco (2011)). However, proton imaging suffers from poor spatial resolution due to the multiple Coulomb scatterings (MCS) undergone by the proton throughout its path (Schneider and Pedroni (1994); Li et al. (2006)). Therefore, when using the unaltered proton radiography, the conventional X-ray CT tomographic algorithm struggles with the MCS and the output images have a very low spatial resolution (2 lp/cm - Li et al. (2006)). Proton path estimate methods have been proposed to account for the MCS effect. The most widely used is the most likely path. It is a probabilistic method to calculate the proton path given directional information and proton beam scattering based on the Fermi-Eyges moments (Schneider and Pedroni (1994); Williams (2004)). To facilitate the technique, a streamlined matrix-based Bayesian formalism has been introduced (Schulte et al. (2008); Erdelyi (2009)) that derives the proton path. These methods require sophisticated detectors to achieve proton-by-proton detection (Schulte et al. (2004); Bashkirov et al. (2007); Hurley et al. (2012)). Fi-

Maximum likelihood proton CT reconstruction

4

nally, Fekete et al. (2015) proposed a phenomenological approach to retrieve the most likely-path from a cubic spline. With a robust path estimator in hand, the pCT reconstruction can be done using a modified photon tomographic algorithm. Amongst the different techniques, the algebraic iterative reconstruction (Hurley et al. (2012)) and the adapted filtered back projection (Rit et al. (2013)) are the most used, both directly yielding the proton RSP map. These methods requires a large number of projections around the phantom to reconstruct the pCT. Furthermore, they do not deal with the problem of poor spatial resolution in the proton radiography itself. To acquire a precise proton radiograph from these algorithms, a full pCT reconstruction is required. The novel method proposed here uses a maximum likelihood estimator (MLE) to reconstruct each projection independently. By doing so, the proton radiograph can be utilized to fulfill two different goals. First, the improved proton radiography can be used for accurate patient alignment prior to the treatment due to the high spatial resolution achievable. It can also be used to trigger recalculation based on drastic change in the patient projected WET. Second, the improved proton radiograph can be used as an input in any classical tomography algorithm. These algorithms would then rapidly yield the patient-specific RSP map. The novel method investigated here demonstrates significant spatial resolution improvement both for proton radiography and pCT as well as great acceleration in the reconstruction.

2. Theory This section will go over the mathematics of the method developed here. Briefly, this work intends to deblur the proton radiography by making use of the estimated proton path. To do so, the object is discretized into voxels and the water equivalent thickness through a voxel projection, defined from the source to the detector pixels, is optimized such that it maximizes the likelihood of the proton energy loss. A voxel projection is labeled as a channel and Figure 1 illustrates the technique . First, let us define concepts that will be used in the following section. S(t) represents the estimated curved path of the proton as a function of the depth in the phantom. The parameter t ∈ [0, 1] represents the fraction of the traversed thickness

Maximum likelihood proton CT reconstruction

5

to the total length. `c,n is the length crossed in channel c by the proton with index n. It is computed using the spatial intersection of the path estimate with the channel’s edge. The index k represents the channel of interest and C the set of all channels including k. The sum of `c,n over all channels crossed defines the total length Ln . The total water equivalent path length of a proton is defined by the WEPLn . The concepts are summarized in Figure 1.

Front tracker! Beam origin!

Rear tracker !

lk,n! Channel K!

Phantom! N-th Proton path!

z

Cubic spline S(t)!

x

y

Figure 1: Schematic view of the scanned phantom. The front and rear tracker are shown with the phantom discretized in channels. The red line represents the proton path, the dotted-line represents the cubic spline estimate and the beige arrows represent the projections from each channel to the rear tracker. On the right, output from the classical technique and our technique are shown.

6

Maximum likelihood proton CT reconstruction

This research aims to estimate the W ETk that maximizes the likelihood of the protons energy. To do so, we define an equation that relates the W EP Ln measured by the n-th proton to the WET of the k-th channel (W ETk ) through an error variable υk,n (Equation 1). υk,n = W EP Ln − W ETk

(1)

Intuitively, a proton that has a longer length in the channel of interest (`k,n ) will have a smaller error υk,n . In the limit of `k,n = Ln , the error is null. In other words, the error υk,n has a high heteroskedasticity (the standard deviation of the error varies as a function of `k,n ). For these reasons, it is assumed that Var[υk,n ] ∝ σk2 /`2k,n . The standard deviation of υk,n is inversely proportional to the length spent in the channel k. To retrieve homoskedasticity (identical standard deviation of errors over all possible value of `k,n ), Equation 1 is multiplied by the `k,n . k,n =

`k,n (W EP Ln − W ETk ) Ln

(2)

where k,n = `k,n υk,n . The error distribution k is specific to the k-th channel and k,n is a single realization of this error distribution for the n-th proton. The initial step to obtain a high resolution pRad is to assume k is a zero-mean normal distribution. (  ) 2 ` 1 k,n N (k ; 0, σk2 ) = p exp − √ (W EP Ln − W ETk )2 (3) 2 2σ L 2πσk k n The problem can then be rewritten as a likelihood estimator. The total likelihood of k and σk2 for a given sample of size N is stated in Equation 4.

L=

N Y n

N (k ; 0, σk2 )

(4)

The value that maximizes the likelihood of that distribution, i.e. the MLE, is found by taking the first derivative of log(L) with respect to W ETk . The derivative of the log likelihood is set to zero, leading to Equation 5. PN W ETk =

n

`2k,n W EP Ln σ12 L2n k `2k,n 1 n L2n σk2

PN

PN =

n

`2k,n W EP Ln L2n PN `2k,n n L2n

(5)

Maximum likelihood proton CT reconstruction

7

3. Material and Methods Three phantoms were studied: a cubic slanted edge phantom for 2-D spatial resolution, an anthropomorphic head phantom for clinical geometry assessment both in 2-D and 3-D and finally the Catphan (The Phantom Laboratory; Salem, New York) phantom for calibration with module CTP528 for 3-D line-pair spatial resolution. Proton radiographies of each phantom were acquired through Geant4 Monte Carlo (Agostinelli et al. (2003)) (V4.10.00) simulations. The detailed technique to produce the radiography is presented in Section 3.1. The radiographies were uniformly simulated around the phantom with a step of 0.5◦ (n=720). 3x106 protons were simulated for the head phantom to demonstrate qualitative improvement, and 107 protons for the slanted edge and the Catphan phantom to quantify the spatial resolution gain. 200 and 330 MeV protons were used to constrain the exit energy to the constant-RSP domain (Arbor et al. (2015)) and to evaluate energy dependence. In order to best represent clinical situations, two types of proton beams were used to produce the radiography, a parallel beam and a conical beam. Tomographic techniques were then used to reconstruct the phantom using the improved proton radiography. The spatial resolution, computed using a slanted-edge phantom (2-D reconstruction) or using the Catphan 528 linepair phantom (3-D reconstruction) is the appropriate metric to assess the technique’s performance. 3.1. Monte Carlo proton radiography Proton radiographies were acquired through Monte Carlo simulations. No detector effects were simulated. Proton propagation through the defined object was simulated using Geant4 (Agostinelli et al. (2003)) Monte Carlo code (v4.10.00). The standard processes includes energy loss and straggling, multiple Coulomb scattering based on Lewis theory (Goudsmit and Saunderson (1940)) and parameterized proton interactions with nucleus and electrons. In precise terms, the following physics lists were enabled: the standard electromagnetic option 3 for higher accuracy of electrons, hadrons and ion tracking without a magnetic field, the hadrons elastic model and the binary ion models both for elastic and inelastic collisions. Cuts were set to 0.1 mm.

Maximum likelihood proton CT reconstruction

8

The head geometry was created from a voxelized X-ray CT. Geant4 tissue material requires both the density and atomic composition of each voxel to produce accurate stopping power tables. The mass density was extracted from the X-ray CT voxel HU, first converted to the relative electron density (ρerel ) using a calibrated look-up table. ρerel was then converted to mass density using Equation 6 (Schneider et al. (1996)). X X wi Zi ρerel · ρwater Ngwater ρ= → Ng = Ngi = NA (6) Ng A i i i The material composition of individual voxel was defined by the generic list of materials detailed in Table 1. Each material has an average HU value that was used as a pivot for assigning atomic composition of voxels. The composition for a particular HU value was defined by taking a linear combination between the two closest pivots, weighted by distance. The Catphan geometry was defined in a parametrical world. The atomic composition, density and the precise geometric definition were summarized in Hansen et al. (2014). The slanted edge phantom is a 20 cm cubic water phantom with a 5 cm aluminum cube in the middle, slanted by 2.5◦ around the beam axis. The output protons were filtered based on their position, direction and exit energy using the 3σ cut suggested in Schulte et al. (2008). 3.2. Proton trajectory estimate: Optimized cubic spline path The proton path estimate is paramount to proton radiography and pCT reconstruction. The path reconstruction is usually based on the knowledge of the initial/final proton position X and direction P. In a recent publication, Fekete et al. (2015) proposed a cubic spline to estimate the proton path. They showed that the cubic spline direction vector magnitude can be optimized to minimize the root-mean square position error with respect to the proton path. To retrieve this direction vector magnitude, they proposed the addition of a factor Λopt that is function of both the proton initial WEPL and the WET crossed. The cubic spline (S(t)) is calculated using Equation 7. S(t) = (2t3 − 3t2 + 1)X0 + (t3 − 2t2 + t)P0 + (3t2 − 2t3 )X1 + (t3 − t2 )P1

(7)

9

Maximum likelihood proton CT reconstruction

Table 1: Density and atomic composition of the materials used as reference for the Monte Carlo simulations (Woodard and White (1986); White et al. (1987)). Each material atomic composition is based on a set of element. Materials for HU within the described range are interpolated in between the closest defined materials. Material Name Lung Deflated Adipose Tissue 3 Adipose Tissue 1 Mean Male soft tissue Muscle Skeletal 1 Muscle Skeletal 2 Skin 1 Connective Tissue Sternum Humerus Spherical head Femur Total Bone Cranium Cortical Bone

HU -741 -98 -55

H 10.4 11.6 11.2

C 10.6 68.3 51.9

N 3.1 0.2 1.3

O 75.7 19.9 35.6

Ca 0.0 0.0 0.0

P 0.2 0.0 0.0

I [eV] 74.54 63.06 66.14

ρe 0.258 0.933 0.970

ρ 0.26 0.93 0.97

5 40 43 72 100 385

10.5 10.2 10.3 10.1 9.5 7.8

25.6 17.3 14.4 25.2 21.0 31.8

2.7 3.6 3.4 4.6 6.3 3.7

60.2 68.7 71.6 59.9 63.1 44.1

0.0 0.0 0.0 0.0 0.0 8.6

0.2 0.2 0.2 0.1 0.0 4.0

72.3 73.69 74.03 72.25 73.79 81.97

1.016 1.040 1.043 1.0 84 1.103 1.211

1.02 1.05 1.05 1.09 1.12 1.25

538 688 999 1524

7.1 6.3 5.0 3.4

38.1 33.4 21.3 15.6

2.6 2.9 4.0 4.2

34.3 36.3 43.8 43.8

12.2 14.4 17.7 22.6

5.6 6.6 8.1 10.4

85.43 90.24 99.69 111.63

1.279 1.355 1.517 1.781

1.33 1.43 1.61 1.92

Λopt can be retrieved from a parameterized fit of the WET/WEPL ratio with the factor C0 , C1 , C2 , C3 and α detailed in Fekete et al. (2015). In this work, the optimized cubic spline was used for proton path reconstruction. 3.3. Water equivalent proton radiography The main quantity of interest in this study is the water equivalent path length (WEPL) of the proton through the medium rather than its energy lost. The energy loss was converted to WEPL by integrating over the stopping power of water (Berger (1993)). The stopping power was calculated using the Bethe-Bloch equation (Equation 8) with the ICRU Units and Measurements (1993) I-value (Iw ) for water (78 eV) as well as their suggested corrections.     dE 4πe4 2 Z 2me c2 β 2 C δ 2 =ρ z ln − β − ln Iw − − (8) − dX me c2 uβ 2 A (1 − β 2 ) Z 2

Maximum likelihood proton CT reconstruction

10

The WEPL is represented by the integral of the inverse of the stopping power, over the energy loss ( Equation 9). Z Ein dE (9) W EP L = Eout −dE/dX(Iwater , β) 3.4. Proton radiography from conical or parallel beam Two beam configurations exist to acquire a proton radiography. First, a proton pencil beam can be used to scan the phantom in a procedure that will henceforth be defined as parallel beam proton radiography or parallel proton radiography. The second method, which is faster, consists of using a broad beam originating from a collimated point source, which will be referred as conical beam proton radiography or conical proton radiography. In our study, the parallel beam point source σF W HM is 5 mm. The conical point source distance is 2 m and the initial angles (φ, θ) are restricted to fully cover the phantom geometry. To optimize the results of the method proposed here, the channel in which the W ETk is extracted was defined so that it follows the average proton path defined from the entrance to the exit detector. The magnification introduced by a conical beam needs to be taken into account when defining the channels geometry for such type of beams. The method used here is shown in Figure 2, where the channels were defined from the point source to the rear detector.

3.5. Tomographic reconstruction algorithm The optimized proton radiography values represent the WEPL projection through the channel. The MLE-enhanced pRad were then used as input in a FDK algorithm (Feldkamp et al. (1984)) to reconstruct the pCT. The implementation was done through the RTK open-source library (Rit et al. (2014)).

11

Maximum likelihood proton CT reconstruction

Beam origin! z!

y!

x! Dout!

θ,φ

Front tracker! Di!

Pi!

Channel K!

Dj!

Pj!

Rear tracker!

Figure 2: Representation of the conical beam pCT reconstruction by taking into account the magnification forms a cone-beam. The distances Di and Dj represent different depths for reconstruction of the detector pixels with initial angles φ, θ and source distance Dout . 4. Results 4.1. W ETk prediction accuracy The assumption that the model proposed in Section 2 produces accurate WET is investigated. The results of the parallel proton beam configuration for WET

12

Maximum likelihood proton CT reconstruction

WET Difference [cm]

Water equivalent thickenss [cm]

30

Proton Radiography DRR CT Improved proton radiography

25 20 15 10 5 0 2

50

100

150 X axis [mm]

200

250

50

100

150 X axis [mm]

200

250

1 0 −1 −2

Figure 3: The upper panel shows the WET thickness projection of the head proton radiography at 0o . The lower panel shows the difference relative to the reference for both the proton radiography and the improved proton radiography. accuracy are shown on Figure 3. The WET originating from a single column of the improved proton radiography was compared to the reference, i.e. the X-ray CT digitally reconstructed radiography (DRR) converted to WET using MC material knowledge. The initial proton radiography was displayed as well for comparison. Figure 3 shows that the WET accuracy, relative to the reference, is ameliorated in the improved proton radiography compared to the normal proton radiography. The WET RMS error between the reference and the improved proton radiography is 0.19 cm, improved from 0.83 cm in the original proton radiography. The errors between the measured and predicted W ETk within the radiography are computed to investigate the bias introduced when the assumption that E[k ] = 0 is not valid. The left panel of Figure 4 displays this error throughout the head phantom between the MC reference and the optimized WET. The edge region shows that the absolute highest error value is ≤0.2 WET-cm. Elsewhere the error approaches zero as theorized. Therefore, the total error introduced by the E[k ] = 0 assumption does

13

Maximum likelihood proton CT reconstruction

0.10 0.05 0.00 −0.05 −0.10 −0.15

WETmeasured - WETpredicted [WET cm]

0.15

0.28 0.24 0.20 0.16 0.12 0.08 0.04

Standard deviation in the WET error (σ) [WET cm]

0.32

0.20

0.00

Figure 4: Left : Two-dimensional map of the mean difference between the measured and predicted water equivalent thickness. Right : Weighted standard deviation of the predicted water equivalent thickness. Both are in unit of WET cm, with 72 entries per channel on average. not cross the 0.2 WET cm threshold. The right panel of Figure 4 represents the weighted standard deviation of the predicted WET for each channel. As expected, the weighted standard deviation increases close to the high-gradient edges, up to 0.3 WET-cm. 4.2. Proton radiography reconstruction 4.2.1. Initial spatial resolution optimization The spatial resolution is quantified through the modulation transfer function (MTF) at 10% of the initial amplitude (MTF10% . For proton radiography, a 20 cm thick water phantom was simulated with a 5 cm slanted aluminum cube in the middle position. The slanted-edge technique (Hansen et al. (2014); Fujita et al. (1992)) was used to evaluate the proton radiography MTF for every configuration. The method proposed here improves the spatial resolution both in the parallel configuration (from 3.49 lp/cm to 5.76 lp/cm) and in the conical configuration (from 3.49 lp/cm to 5.56 lp/cm), as shown in Figure 5. A clinical voxelized anthropomorphic phantom was analyzed to demonstrate the improvement achievable in the clinic. The voxelized head phantom radiography was

14

Maximum likelihood proton CT reconstruction

1.0 Improved conical pRad, (M T F10% = 5.56 lp/cm) Initial conical pRad, (M T F10% = 3.49 lp/cm) Improved parallel pRad, (M T F10% = 5.76 lp/cm)

Intensity [%]

0.8

Initial parallel pRad, (M T F10% = 3.49 lp/cm)

0.6

0.4

0.2

0.0

0

1

2

3 4 5 Spatial resolution [lp/cm]

6

7

8

Figure 5: Modulation transfer function for the normal proton radiography and the improved proton radiography using both the parallel and conical beam configurations. A linear regression to a sigmoid fit has been done to reduce noise. The MTF10% is the conventional high frequency reference to compare various technique and is displayed for each configuration. investigated and reconstructed using the MLE method. The WET X-ray CT DRR is used as a reference. Figure 6 displays the normal proton radiography, the improved proton radiography and the WET-DRR reference. A zoom box is shown on the high gradient region to display the improvement achieved. Results from the conical beam reconstructed proton radiography have not been included since they are qualitatively very similar to Figure 6.

4.2.2. Relationship between pRad spatial resolution and the initial proton energy The MLE method was then investigated as a function of the beam energy for 200 MeV and 330 MeV protons. Again, the slanted edge phantom has been used to measure spatial resolution. The MLP uncertainty is inversely proportional to the energy and consequently the path estimator is expected to perform better with higher energies,

15

Maximum likelihood proton CT reconstruction

A) pRad - SLP

B) pRad - Opt

C) X-ray CT DRR

Figure 6: A) Proton radiography of the head phantom. B) Improved proton radiography with zoom on a sharp gradient structure. C) Digitally reconstructed radiography of the head CT. effectively leading to higher spatial resolutions. Figure 7 demonstrate the spatial resolution using the slanted-edge phantom for both initial energies. As hypothesized, the spatial resolution increases with the beam initial energy. 4.3. Proton CT reconstruction The improved proton radiographies were used to perform a full 3-D reconstruction using the FDK algorithm for the head and the Catphan phantoms. First, the 3-D Catphan phantom has been reconstructed using 720 projections (n=107 protons/projection) to quantify the spatial resolution achievable. For brevity, only the parallel beam configurations are shown on Figure 8. The pCT MTF is evaluated through the line-pairs on the CTP528 module, both for the conical beam and the parallel beam. The reconstructed parallel beam pCT displays a slightly better spatial resolution (from 2.79 to 5.55 lp/cm) than the conical beam pCT (from 3.03 to 5.15 lp/cm) as seen on Figure 9). Figure 10 displays the clinical improvement achievable from this technique in pCT reconstruction by comparing a superior-inferior slice of the head pCT and X-ray CT. The pCT slice compares favorably with the X-ray CT, although some blurring of the

16

Maximum likelihood proton CT reconstruction

1.0 200 MeV improved pRad, (M T F10% = 4.53 lp/cm) 200 MeV initial pRad, (M T F10% = 2.44 lp/cm) 330 MeV improved pRad, (M T F10% = 5.76 lp/cm)

Intensity [%]

0.8

330 MeV initial pRad, (M T F10% = 3.49 lp/cm)

0.6

0.4

0.2

0.0

0

1

2

3 4 5 Spatial resolution [lp/cm]

6

7

8

Figure 7: Modulation transfer function for proton radiography of the slanted-edge phantom using two different initial energy, 200 and 330 Mev. 107 protons have been used to produce the curve and a sigmoid curve regression has been done to reduce noise. image persists. Again, the proton conical beam configuration was not displayed but yielded very similar improvement.

5. Discussion This study proved that it is possible to separate the MCS in the proton path from the projection of the materials stopping power, using the MLE method described in Section 2. Classically, a proton loses its energy proportionally to the length spent in each channel as well as their WET. However, due to the actual detector configurations, the energy loss is only accounted in the last channel traversed rather than throughout its path. This effect causes a loss of spatial resolution and produces the blurring seen in proton imaging. This project aims to redistribute the proton path into the appropriate channels, based on the fraction of time spent in each individually

Maximum likelihood proton CT reconstruction

17

Figure 8: FDK reconstructed parallel beam pCT of the Catphan CTP528 module as calculated with a path estimator. By doing so, the new MLE-enhanced proton radiography represents the projected stopping power of the material channels with higher spatial resolution. Although the method developed in Section 2 considers a zero-mean normal distribution of the errors, which is known to become inaccurate around high-gradient edges, the introduced bias in those regions is minimal (≤ 0.2 WET cm). An iterative reconstruction could help refine the spatial resolution at the cost of computational time. The proton radiographies have been reconstructed with a gain in spatial resolution from 3.49 lp/cm to 5.76 lp/cm for the parallel beam and from 3.49 lp/cm to 5.56 lp/cm for the conical beam. For the parallel beam reconstruction, the WET rootmean square error with the reference WET-DRR was reduced from 0.83 cm in the normal proton radiography to 0.19 cm in the improved radiography. The method has also been applied to a voxelized anthropomorphic X-ray CT with good improvement in the high-gradient region (Figure 6). Seco et al. (2013) have shown that the proton radiography spatial resolution is highly dependent on the proton initial energy as well as the thickness crossed. For 200 MeV protons, they obtained a proton radiography spatial resolution of 2.3 lp/cm, similar to what is found here before the improvement

18

Maximum likelihood proton CT reconstruction

Improved parallel pRad, M TF10% = 5.55 lp/cm

1.0

Improved conical pRad, M TF10% = 5.15 lp/cm Initial parallel pRad,

M TF10% = 2.79 lp/cm

Initial conical pRad

M TF10% = 3.03 lp/cm

Intensity [%]

0.8

0.6

0.4

0.2

0.0

0

1

2

3 4 5 Spatial Resolution [lp/cm]

6

7

8

Figure 9: Three-dimensional reconstructed modulation transfer function from the Catphan CTP528 line-pair phantom for both the conical and the parallel beam, in both the initial and optimized configurations. A regression to a sigmoid has been done to reduce noise.

A) pCT - SLP

B) pCT - Opt

C) X-Ray CT

Figure 10: A) Reconstructed pCT of the head phantom from the original proton radiography, compared to the reconstruction with B) the optimized proton radiography and C) the same slice of the X-ray CT phantom.

Maximum likelihood proton CT reconstruction

19

in the 200 MeV case (2.44 lp/cm). Throughout this work, energy straggling has been neglected when considering the proton energy loss. For low energy loss events, the energy straggling should have a Gaussian behavior and average out. However, in the case of high energy loss events, the energy straggling will adopt a Landau distribution that will directly affect the reconstruction by introducing a shift in the reconstructed WET value. This uncertainty is not accounted for in normal pCT reconstruction. The MLE-enhanced pRad represents the linear projection of the channel WET. They were used as input into the FDK tomographic reconstruction algorithm (Feldkamp et al. (1984)), with the RTK library (Rit et al. (2014)). The results demonstrate high spatial resolution (5.55 lp/cm parallel and 5.15 lp/cm conical, from 2.79 and 3.03 lp/cm in the non-optimized proton CT configurations) and high quality reconstruction for anthropomorphic images. The reconstructed Catphan pCT (Figure 8) results are qualitatively comparable to Rit et al. (2013). Still, their characterization of the spatial resolution through the 10%-90% profile of a gradient differs from the techniques used here and is not comparable. Their method uses a convex-Hull to improve the path estimate and does not rely on the zero-mean assumption. It should therefore output a more accurate pCT. However, a convex-Hull is not applicable here. The definition of the convex-Hull in a complicated phantom, e.g. the anthropomorphic head, would require some form of adjunct imaging technique or pre-reconstruction that would make the reconstruction of individual radiographs pointless. Hansen et al. (2014) have investigated improved pCT by the use of combined proton projection and cone beam CT imaging (CBCT). In their work, they reported an MTF50% rather than the MTF10% investigated in this work. Their report shows an MTF50% of 5.22 lp/cm, compared to our MTF50% of 3.24 lp/cm for the parallel beam and 2.82 lp/cm for the conical beam. Their method does produce better resolution using photon imaging, but requires the full apparatus of the CBCT imaging and has not been tested in a clinical scenario where the setup would change in between the two techniques. Li et al. (2006) have analyzed the MTF using the most likely path estimate and the algebraic reconstruction technique for an in-house phantom they designed. They stated a spatial resolution of 5 lp/cm, slightly below the resolution found here. However, they do not specify whether it is the MTF10% they analyzed which makes this comparison weak. If it is the MTF10%

Maximum likelihood proton CT reconstruction

20

they provide, the results are in the same domain as ours for the conical beam. Finally, Cirrone et al. (2011) have applied severe filters to the protons, keeping only those which go in straight-line through the channels. With their method, they obtained a maximum spatial resolution of 3.3 lp/cm for 200 MeV protons, below the results obtained here (4.53 lp/cm - Figure 7). In addition, their technique is expected to suffer from dose problem as they need a large number of protons that pass their severe filter in order to reach statistical reliability. The algorithm is highly parallelizable as every single projection and every single proton is independent. The reconstruction time the algorithm requires to improve a single proton radiography is 4 min 22 seconds per bunch of 107 protons. This number of protons is sufficient to produce accurate results. GPU acceleration could easily provide very fast and reliable improved proton images. For clinical use, the angular divergence of the clinical proton beam would need to be characterized. Parallel beam considered here are useful for numerical study, but unrealistic as even a pencil beam has a small deviation. The conical beam should be applied in a realistic scenario. A major benefit of the presented technique, compared to other reconstruction techniques, is that it can produce an improved proton radiography independent of the pCT reconstruction. First, this improved proton radiography can help to validate the channel WET for quality assurance of the planned beam. Second, the improved proton radiography would increase the accuracy of patient alignment to minimize the motion uncertainty. The fact that the image is acquired in the same position setup as the treatment helps greatly. Finally, a complete iterative algebraic proton CT reconstruction requires the construction of a system matrix. The system matrix corresponds to the integral of each particle track through each voxel, sampling their relative stopping power. The image is then reconstructed through an iterative weighted least square algorithm (argmin||Au − B||22 ), where A corresponds to the system matrix, B is a vector of the measured data and u is the required RSP map. It would be possible to adapt this approach to proton radiography by downsizing the system matrix to channels rather than voxel and using only the measured data of a single proton radiography to construct the B vector. This alternative method would then produce high spatial resolution proton radiography. It is expected to be the optimal solution to the

Maximum likelihood proton CT reconstruction

21

problem. It would however suffer from the same time challenges that exist in the proton CT iterative algebraic reconstruction. 6. Conclusion In this project, a fast and precise technique has been developed that improves proton radiography spatial resolution significantly, from 3.49 lp/cm to 5.76 lp/cm, a gain of 65%. The method has been shown to produce accurate radiographies that represent the WET seen by a particle having a linear trajectory when crossing the material. These improved radiographies have been used to reconstruct three-dimensional pCT images with resolution of 5.55 lp/cm. The methods presented here work similarly for beams coming from either parallel pencil or conical broad beams, with resolution of 5.56 lp/cm and pCT resolution of 5.15 lp/cm.

7. Acknowledgement The authors would like to acknowledge Zachary Slepian for helpful discussions and help in text revision. Charles-Antoine Collins-Fekete is supported by a scholarship from Fonds de Recherche du Qu´ebec – Nature et technologies. CACF acknowledges partial support by the CREATE Medical Physics Research Training Network grant of the Natural Sciences and Engineering Research Council (Grant number: 432290). Computations were made on the supercomputer Colosse from Universit´e Laval, managed by Calcul Qu´ebec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), NanoQu´ebec, RMGA and the Fonds de Recherche du Qu´ebec - Nature et Technologies (FRQ-NT).

REFERENCES

22

References Agostinelli S, Allison J, Amako K, Apostolakis J, Araujo H, Arce P, Asai M, Axen D, Banerjee S, Barrand G, Behner F, Bellagamba L, Boudreau J, Broglia L, Brunengo A, Burkhardt H, Chauvie S, Chuma J, Chytracek R, Cooperman G, Cosmo G, Degtyarenko P, Dell’Acqua A, Depaola G, Dietrich D, Enami R, Feliciello A, Ferguson C, Fesefeldt H, Folger G, Foppiano F, Forti A, Garelli S, Giani S, Giannitrapani R, Gibin D, G´omez Cadenas J, Gonz´alez I, Gracia Abril G, Greeniaus G, Greiner W, Grichine V, Grossheim A, Guatelli S, Gumplinger P, Hamatsu R, Hashimoto K, Hasui H, Heikkinen A, Howard A, Ivanchenko V, Johnson A, Jones F, Kallenbach J, Kanaya N, Kawabata M, Kawabata Y, Kawaguti M, Kelner S, Kent P, Kimura A, Kodama T, Kokoulin R, Kossov M, Kurashige H, Lamanna E, Lamp´en T, Lara V, Lefebure V, Lei F, Liendl M, Lockman W, Longo F, Magni S, Maire M, Medernach E, Minamimoto K, Mora de Freitas P, Morita Y, Murakami K, Nagamatu M, Nartallo R, Nieminen P, Nishimura T, Ohtsubo K, Okamura M, O’Neale S, Oohata Y, Paech K, Perl J, Pfeiffer A, Pia M, Ranjard F, Rybin A, Sadilov S, Di Salvo E, Santin G, Sasaki T, Savvas N, Sawada Y, Scherer S, Sei S, Sirotenko V, Smith D, Starkov N, Stoecker H, Sulkimo J, Takahata M, Tanaka S, Tcherniaev E, Safai Tehrani E, Tropeano M, Truscott P, Uno H, Urban L, Urban P, Verderi M, Walkden A, Wander W, Weber H, Wellisch J, Wenaus T, Williams D, Wright D, Yamada T, Yoshida H and Zschiesche D 2003 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506(3), 250– 303. Arbor N, Dauvergne D, Dedes G, L´etang J M, Parodi K, Qui˜ nones C T, Testa E and Rit S 2015 Physics in Medicine and Biology 60(19), 7585–7599. Bashkirov V, Schulte R, Penfold S and Rosenfeld A 2007 in ‘IEEE Nuclear Science Symposium Conference Record, 2007. NSS ’07’ Vol. 6 pp. 4685–4688. Berger M J 1993 Penetration of Proton Beams Through Water 1. Depth-dose Distribution, Spectra and LET Distribution U.S. Department of Commerce, National Institute of Standards and Technology. Chvetsov A V and Paige S L 2010 Physics in Medicine and Biology 55(6), N141–149. Cirrone G A P, Bucciolini M, Bruzzi M, Candiano G, Civinini C, Cuttone G, Guarino P, Lo Presti D, Mazzaglia S E, Pallotta S, Randazzo N, Sipala V,

REFERENCES

23

Stancampiano C and Talamonti C 2011 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 658(1), 78–83. Cormack A M 1963 Journal of Applied Physics 34(9), 2722–2727. Cormack A M and Koehler A M 1976 Physics in Medicine and Biology 21, 560–569. Depauw N and Seco J 2011 Physics in Medicine and Biology 56(8), 2407–2421. Erdelyi B 2009 Physics in Medicine and Biology 54(20), 6095–6122. Fekete C A C, Doolan P, Dias M F, Beaulieu L and Seco J 2015 Physics in Medicine and Biology 60(13), 5071–5082. Feldkamp L A, Davis L C and Kress J W 1984 Journal of the Optical Society of America A 1, 612–619. Fujita H, Tsai D Y, Itoh T, Morishita J, Ueda K, Ohtsuka A and others 1992 Medical Imaging, IEEE Transactions on 11(1), 34–39. Goudsmit S and Saunderson J L 1940 Physical Review 57(1), 24–29. Hansen D C, Petersen J B B, Bassler N and Sørensen T S 2014 Medical Physics 41(3), 031904. Hanson K M, Bradbury J N, Cannon T M, Hutson R L, Laubacher D B, Macek R J, Paciotti M A and Taylor C A 1981 Physics in Medicine and Biology 26(6), 965– 983. Hurley R F, Schulte R W, Bashkirov V A, Wroe A J, Ghebremedhin A, Sadrozinski H F W, Rykalin V, Coutrakon G, Koss P and Patyal B 2012 Medical Physics 39(5), 2438–2446. Li T, Liang Z, Singanallur J V, Satogata T J, Williams D C and Schulte R W 2006 Medical Physics 33(3), 699. Matsufuji N, Tomura H, Futami Y, Yamashita H, Higashi A, Minohara S, Endo M and Kanai T 1998 Physics in Medicine and Biology 43(11), 3261–3275. Rit S, Dedes G, Freud N, Sarrut D and L´etang J M 2013 Medical Physics 40(3), 031103. Rit S, Oliva M V, Brousmiche S, Labarbe R, Sarrut D and Sharp G C 2014 in ‘Journal of Physics: Conference Series’ Vol. 489 IOP Publishing p. 012079. Schaffner B and Pedroni E 1998 Physics in Medicine and Biology 43(6), 1579–1592. Schneider U and Pedroni E 1994 Medical Physics 21(11), 1657–1663.

REFERENCES

24

Schneider U, Pedroni E and Lomax A 1996 Physics in medicine and biology 41(1), 111. Schulte R, Bashkirov V, Tianfang Li, Zhengrong Liang, Mueller K, Heimann J, Johnson L, Keeney B, Sadrozinski H W, Seiden A, Williams D, Lan Zhang, Zhang Li, Peggs S, Satogata T and Woody C 2004 IEEE Transactions on Nuclear Science 51(3), 866–872. Schulte R W, Bashkirov V, Klock M C L, Li T, Wroe A J, Evseev I, Williams D C and Satogata T 2005 Medical Physics 32(4), 1035–1046. Schulte R W, Penfold S N, Tafas J T and Schubert K E 2008 Medical Physics 35(11), 4849–4856. Seco J, Oumano M, Depauw N, Dias M F, Teixeira R P and Spadea M F 2013 Medical Physics 40(9), 091717. Seltzer S M and Berger M J 1982 The International Journal of Applied Radiation and Isotopes 33(11), 1189–1218. Units I C o R and Measurements 1993 Stopping powers and ranges for protons and alpha particles International commission on radiation units and measurements. White D R, Woodard H Q and Hammond S M 1987 The British Journal of Radiology 60(717), 907–913. Williams D C 2004 Physics in Medicine and Biology 49(13), 2899. Woodard H Q and White D R 1986 The British Journal of Radiology 59(708), 1209– 1218. Yang M, Zhu X R, Park P C, Titt U, Mohan R, Virshup G, Clayton J E and Dong L 2012 Physics in Medicine and Biology 57(13), 4095–4115. Zygmanski P, Gall K P, Rabin M S Z and Rosenthal S J 2000 Physics in Medicine and Biology 45(2), 511.