Feriel Khellaf
email: feriel.khellaf@creatis.insa-lyon.fr

Nils Krah

Jean Michel Létang

Simon Rit

2D directional ramp filter

2D directional ramp filter

INTRODUCTION

The filtered backprojection (FBP) is the reference algorithm to reconstruct tomographic data. It consists in (1) filtering the projections with a 1D ramp filter and (2) backprojecting the result in the image domain. An alternative method, the backproject-filter or backproject-first (BPF), starts with the backprojection before filtering in the 2D image domain. BPF approaches can be useful when dealing with certain modalities to avoid resampling. For example, in proton CT (pCT) (Figure 1 center), starting by backprojecting the protons along their non-linear paths could improve spatial resolution, compared with binning the list-mode data into projections [START_REF] Zeng | Image reconstruction via the finite Hilbert transform of the derivative of the backprojection[END_REF][START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojection-then-filtering approach[END_REF][START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF]. Similarly, time-of-flight information of positron emission tomography (PET) can be accounted for in the image domain [START_REF] Watson | An evaluation of image noise variance for time-of-flight PET[END_REF] (Figure 1 right). Our objective is to develop such an approach by extending the usual 1D ramp filter in 2D. BPF methods usually rely on 2D radial filters [START_REF] Zeng | Can the backprojection filtering algorithm be as accurate as the filtered backprojection algorithm?[END_REF][START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojection-then-filtering approach[END_REF]. However, the backprojected image has an infinite support which is problematic during convolution since the backprojection matrix must have a finite size. The truncation will cause an offset in the reconstructed values. This effect can be minimized by using a very large backprojection matrix. Another possible backproject-first approach is based on the Hilbert transform inversion [START_REF] Noo | A two-step Hilbert transform method for 2D image reconstruction[END_REF][START_REF] Zeng | Image reconstruction via the finite Hilbert transform of the derivative of the backprojection[END_REF][START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF], and involves a differentiation step.

In this paper, the idea is to (1) smear each projection in the image space, (2) filter the smeared projection with a 2D version of the 1D ramp filter and (3) integrate over all projection angles. A similar 2D version of the ramp filter has been proposed by [START_REF] Watson | An evaluation of image noise variance for time-of-flight PET[END_REF] in the context of time-of-flight positron emission tomography (TOF PET), although the filter was sampled directly in the frequency domain causing sampling artifacts. A variation of this filter, also proposed by [START_REF] Watson | An evaluation of image noise variance for time-of-flight PET[END_REF] and called CRG-CW (Convolved Ramp and Gaussian with Confidence Weighting), has been used for TOF direct reconstruction [START_REF] Kao | Windowed image reconstruction for time-of-flight positron emission tomography[END_REF][START_REF] Guo | Evaluation of windowed image reconstruction for time-of-flight PET[END_REF][START_REF] Conti | Estimating image quality for future generations of TOF PET scanners[END_REF][START_REF] Kao | Analytic reconstruction methods for list-mode time-of-flight PET[END_REF]. Similarly to what is done for the 1D ramp filter, we sample the filter in the spatial domain to avoid sampling artifacts. We start by defining the 2D ramp filter in the case of parallel X-ray CT (xCT), derive its expression in the spatial domain, before applying this method to simulated xCT and pCT data. Although the main applications in mind for this method are pCT and TOF PET, the filter is mathematically demonstrated in a conventional xCT setting to avoid technicalities specific to proton/PET imaging. It is then applied to proton CT simulations.

METHODS

Parallel FBP and 1D ramp filter

An image f (x 1 , x 2 ), with (x 1 , x 2 ) the spatial coordinates, can be reconstructed from its parallel projections

p φ (s 1 ) = R f (s 1 θ + s 2 θ ⊥ ) ds 2 , (1)
where s 1 is the distance of the projection line from the center, s 2 is the position along the projection line, θ = (cos φ, sin φ) and θ ⊥ = (-sin φ, cos φ) (see Figure 1). Using a filtered backprojection method, the image is reconstructed by

f (x 1 , x 2 ) = π 0 pφ (x 1 cos φ + x 2 sin φ) dφ, (2)
where

pφ (s 1 ) = R P φ (σ 1 )|σ 1 |e i2πσ 1 s 1 dσ 1 , (3)
with P φ (σ 1 ) defined as the Fourier transform (FT) of p φ (s 1 ).

In practice, the filtered projections pφ (s 1 ) are obtained using a band-limited version of the 1D ramp filter. A common choice proposed by [START_REF] Bracewell | Inversion of fan-beam scans in radio astronomy[END_REF] and [START_REF] Ramachandran | Three-dimensional reconstruction from radiographs and electron micrographs: application of convolutions instead of Fourier transforms[END_REF] is

H(σ 1 ) = |σ 1 | if -W < σ 1 < W 0 otherwise, (4)
with W the Nyquist frequency defined as W = 1/2τ with τ the sampling period of the projections p φ . Directly discretizing this expression to filter the projections in the frequency domain would cause an underestimation of the entire reconstructed values, because of the zero value attributed to the zero frequency bin [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF][START_REF] Crawford | CT filtration aliasing artifacts[END_REF]. In order to avoid this issue, the impulse response of the filter is calculated by taking its inverse Fourier transform

h(s 1 ) = R H(σ 1 )e i2πσ 1 s 1 dσ 1 = 1 2τ 2 sinc 2πs 1 2τ - 1 4τ 2 sinc 2 πs 1 2τ (5)
and subsequently discretized by replacing s 1 = nτ , which yields

h[n] =          1 4τ 2 , if n = 0 0, if n even -1 n 2 π 2 τ 2 , if n odd, (6)
as was proposed by [START_REF] Ramachandran | Three-dimensional reconstruction from radiographs and electron micrographs: application of convolutions instead of Fourier transforms[END_REF]. The discrete Fourier transform (DFT) of this impulse response does not have a zero value at its zero frequency and allows an accurate reconstruction of the image values [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. Equation 6 is used to compute the 1D ramp filter used for FBP.

From the 1D ramp filter to a 2D directional ramp filter

From Equation 2, we see that the projections are filtered then smeared in the image space using the variable change s 1 = x 1 cos φ + x 2 sin φ. We want to apply a 2D version of the ramp filter on the already smeared projections, i.e. do the variable change before filtering.

We recall that (σ 1 , σ 2 ) are the Fourier variables corresponding to (s 1 , s 2 ) in the spatial domain. We introduce the dirac δ function, using the property

R δ(σ 2 )e i2πσ 2 s 2 dσ 2 = 1, in Equation 3 pφ (s 1 ) = R R δ(σ 2 )P φ (σ 1 )|σ 1 |e i2π(σ 1 s 1 +σ 2 s 2 ) dσ 1 dσ 2 (7)
and, given the expression of the Fourier transform P φ (σ 1 ) and the fact that formally

δ(σ 2 ) = R e -i2πσ 2 s 2 ds 2 , note that δ(σ 2 )P φ (σ 1 ) = R R p φ (s 1 )e -i2π(σ 1 s 1 +σ 2 s 2 ) ds 1 ds 2 . (8)
We use the following vector notations for s 1 = x • θ, s 2 = x • θ ⊥ and the corresponding Fourier variables σ 1 = ξ • θ and σ 2 = ξ • θ ⊥ , where x = (x 1 , x 2 ) and ξ = (ξ 1 , ξ 2 ). We also define b φ : R 2 → R, the backprojection of one projection, as

b φ (x) = p φ (x • θ) (9)
and B φ (ξ) its Fourier transform. We call b φ (x) a smeared projection. Using these definitions, we note that Equation 8 is the expression of the Fourier transform of b φ (x), and Equation 7becomes bφ

(x) = R 2 B φ (ξ)|ξ • θ|e i2πξ•x dξ, (10)
and the image can then be reconstructed from

f (x) = π 0 bφ (x) dφ. ( 11
)
This equation allows reconstructing the image in three steps

(1) Smear each projection in the image space, i.e. compute a 2D image b φ (x) for each projection angle;

(2) Filter each 2D image b φ (x) with the corresponding 2D filter |ξ • θ| to obtain bφ (x);

(3) Sum bφ (x) over all angles φ.

Note that this 2D version of the ramp filter depends on φ, hence the name directional ramp (DR) filter. By definition, the backprojection operation consists in smearing each projection in the image space, i.e. interpolating b φ (x) from p φ (s 1 ) using b φ (x) = p φ (x • θ) and summing over all projection angles, i.e. π 0 b φ (x)dφ . In the proposed method, we separate these two steps of the backprojection: steps (1) and (3) correspond to the backprojection and step (2) to the filtering. This method resembles a BPF approach since we start by smearing the projections in the image space, however the angular sum to compute the backprojected image is the last step as in an FBP. The usual 2D cone filter used in BPF methods is ξ , and it is applied after the full backprojection since it is not direction dependent, unlike the DR filter which must be used before the angular sum as it depends on the projection angle. A similar filter to the 2D DR, called transverse ramp, has been used to reconstruct TOF PET data [START_REF] Watson | An evaluation of image noise variance for time-of-flight PET[END_REF]. It was directly sampled in the frequency domain, causing sampling artifacts due to the zero values along the line ξ • θ = 0. We choose to calculate the filter's impulse response as in the 1D case to avoid such artifacts.

Impulse response of the 2D directional ramp filter

A band-limited version of the 2D DR filter is given by

H φ (ξ) = |ξ • θ| if | ξ 1 |, | ξ 2 |< W 0 otherwise. ( 12
)
The corresponding kernel in the direct space is obtained by taking the inverse Fourier transform of this directional filter:

h φ (x) = R 2 H φ (ξ)e i2πξ•x dξ. (13)
For φ ∈ [-π/4, π/4], the integral can be split into two parts,

h φ (x 1 , x 2 ) = W -W -ξ 2 tan φ -W -(ξ 1 cos φ + ξ 2 sin φ)e i2π(ξ 1 x 1 +ξ 2 x 2 ) dξ 1 + W -ξ 2 tan φ (ξ 1 cos φ + ξ 2 sin φ)e i2π(ξ 1 x 1 +ξ 2 x 2 ) dξ 1 dξ 2 , ( 14
)
knowing that ξ • θ is negative below the line given by ξ 1 = -ξ 2 tan φ and positive above it. These integrals lead to

h φ (x 1 , x 2 ) = -x 2 cos φ sinc π(x 1 tan φ -x 2 ) τ -x 1 sin φ cos πx 1 τ cos πx 2 τ + (x 1 sin φ + x 2 cos φ) sinc πx 2 τ cos πx 1 τ 1 2τ (πx 1 ) 2 x 2 + cos φ 2τ 3 sinc πx 1 τ sinc πx 2 τ . (15)
The discretized impulse response, with x 1 = mτ and x 2 = nτ , is given by The expressions for m and/or n equal to zero are found by setting x 1 and/or x 2 to zero in Equation 14. The expression can be extended to angles φ / ∈ [-π/4, π/4] using

h φ [m, n] =                      1 2π 2 τ 3 -cos φ m 2 sinc(π(m tan φ -n)) -cos(πn) cos(πm) sin φ nm if m, n = 0 sin 2 φ 2π 2 n 2 τ 3 cos φ cos(πn) if m = 0, n = 0 cos φ 2π 2 m 2 τ 3 (cos(πm) -sinc(πm tan φ)) if m = 0, n = 0 2 cos 2 φ + 1 12τ 3 cos φ if m = 0, n = 0. (16)
h φ+π [m, n] = -h φ [m, n] and h φ+π/2 [m, n] = h -φ [n, m].
Note that setting φ = 0 in Equation 16gives the expression of the 1D ramp filter multiplied by 1/τ on the line n = 0, and zero elsewhere. It should also be noted that the impulse response of the filter (Equation 15) has infinite support, such that sampling the filter in the frequency domain would always result in aliasing artifacts according to the sampling theorem [START_REF] Crawford | CT filtration aliasing artifacts[END_REF].

Size of the backprojection region

In order for the 2D filtered smeared projections to be equivalent to the 1D filtered projections, we need non-truncated images b φ (x) along θ inside the field of view (FOV) (Figure 2). We do not need the full projections outside the FOV since the DR filter tends to a dirac along θ ⊥ in the direct space (since it is constant in the Fourier space). By taking a backprojection region (region where the projections will be smeared) of the same size as the reconstruction region, part of the data near the borders will be truncated for angles φ = kπ/2, k ∈ Z. In order to keep the full projections inside a circular FOV for all angles, one should set M = √ 2N , where the backprojection region is M × M and the reconstruction region is N × N .

Simulations

We generated parallel ray projections of the 2D Shepp Logan phantom [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF] using RTK [START_REF] Rit | The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK)[END_REF] in order to test the filter. 804 projections were acquired over a 180 degree range, with a detector of 512 pixels. The phantom was reconstructed using grids of 512 × 512 pixels and 1024 × 1024 pixels. Poisson noise was applied to the idealized projections, using an initial fluence of 10 7 protons per detector pixel, and weighting the line integrals with 0.01879 mm -1 , the linear attenuation coefficient of water at 75 keV, to work with physical attenuation values.

A proton CT simulation was performed using the Monte Carlo platform GATE [START_REF] Jan | GATE V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy[END_REF]. The setup is detailed in [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF]. A fan beam proton source was simulated, and the list-mode data was binned into parallel geometry images b φ (x) using the same binning as in [START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF]. The non-linearity of proton trajectories was taken into account by estimating each proton's most likely path [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. We used a spiral phantom, consisting of a disk of water with aluminum inserts (Figure 5). The spatial resolution of each insert was obtained by extracting profiles across the insert edges and determining the frequency corresponding to 10% of the peak value of the modulation transfer functiion (MTF) of these profiles. The reconstruction using the DR filter was compared with the algorithm presented in [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF] which generalizes the FBP to non-linear paths.

RESULTS

Figure 3 shows the "pseudo-continuous" filter, i.e. the oversampled filter calculated from the continuous impulse response (Equation 15) by setting x 1 = mτ and x 2 = nτ with τ τ . It is interesting to observe that, after discretization, for the case φ = 0 (left), the samples outside the central horizontal line lie exactly on the zero crossings. For the case φ = π/3 (right), the filter is not zero outside the filtering line n = m tan φ. The DFT of the filters looks like the usual 1D ramp filter duplicated along the direction θ ⊥ .

The reconstructions of the Shepp-Logan phantom using different parameters are shown in Figure 4. The reconstructed values using the filter sampled in the frequency domain are underestimated. The bias is constant and is equal to -0.17. The rest of the reconstructions are very similar, although some artifacts are visible near the outer edges of the phantom in the case of the non-apodized DR filter for idealized projections. This is confirmed on the difference maps, where the DR reconstructed image is slightly noisier for idealized projections, although this is only visible with a tight color scale ([-0.01, 0.01]). Using an oversampled backprojection (but keeping the same sampling for the projections) reduces the noise, and apodizing the filter with a Hamming window completely suppresses the artifacts. The reconstructions from noisy projections show the DR filter performs as well as the 1D ramp filter, and that the previous artifacts are hidden by noise.

The reconstructions of proton CT data are shown in Figure 5. There is no visible difference between the image reconstructed with an FBP and with the DR filter. The spatial resolution f MTF10% as a function of each bead's distance to the center of the phantom is represented. The difference between the resolution for the FBP and the DR algorithms reaches 0.8 lp/mm for the outermost bead.

DISCUSSION

Our results show that the 2D DR filter, sampled in the spatial domain, allows accurate reconstruction of xCT and pCT data, although some artifacts are apparent in the case of reconstructions from noiseless xCT projections. We know that, in the direct space, the DR filter should look like a dirac in the direction θ ⊥ for the 2D convolution to be equivalent to the 1D one. However, apodizing the filter with a square window in Equation 12 results in a convolution with a 2D sinc function, which explains the shape of the filter in Figure 3. While the discretized filter is a dirac for angles φ = kπ/2, k ∈ Z, since all samples outside the central line are zero, this is not the case for other angles, meaning there will be interferences from other lines during convolution.

In Fourier space, this translates in a discontinuity of the FT at the Nyquist frequency. This is visible in Figure 3 where the DFT of the DR filter for φ = π/3 has very thin black lines on its bottom left corner. This could explain the remaining artifacts visible in Figure 4. A similar issue with the continuity of the FT was reported by [START_REF] Schöndube | Evaluations of a 2D inverse Hilbert Transform[END_REF], where the authors define a 2D inverse Hilbert transform. The image reconstructed using this filter suffered from ring artifacts. One of the proposed solutions is to apodize the filter in the Fourier space using a Hamming window to deal with the discontinuity. Our results show that the apodization effectively suppresses artifacts, although it also causes a slight loss of spatial resolution. An alternative is to increase the sampling of the backprojected image [START_REF] Schöndube | Evaluations of a 2D inverse Hilbert Transform[END_REF], which also reduces the noise as seen in Figure 4. Since the artifacts are related to the sampling of the 2D DR filter, taking more samples i.e. increasing the filter's band limit, will make the reconstruction less noisy. However, this artifact is not visible when using projections with a realistic level of noise.

As was observed for the 1D ramp filter by [START_REF] Crawford | CT filtration aliasing artifacts[END_REF] and [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF], sampling in the frequency domain will cause aliasing artifacts (an underestimation in the spatial domain) due to the infinite support of the impulse response.

In the case of simple line integrals, such as in xCT, using the 2D DR filter should give the same result as a standard FBP but is computationally less effective. The DR filter could rather be used with modalities such as proton CT or TOF PET, where the projections are respectively integrals along non-linear paths due to multiple Coulomb scattering or integrals along weighted lines due to time resolution of the detector. In xCT, the projections are constant along the direction θ ⊥ , therefore using 1D projections and the 1D ramp filter makes sense. Both pCT and TOF PET add a new dimension, the depth of the proton or the time-of-flight, such that the smeared projections vary along θ ⊥ (see Figure 1). Filtering the 2D smeared projections b φ (x) requires either rotating each image so that its samples are aligned with the 1D filter, then perform a line-by-line filtering; or using a directional 2D ramp filter. The DR filter could present some notable advantages in this case: as in BPF methods, the filtering is done in the image space so that we avoid rotations that can cause a loss of spatial resolution; at the same time, unlike other BPF approaches, there is no bias due to the truncation of the backprojection region which does not need to be approximated using very large matrices (2 to 4 times the reconstruction region in [START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojection-then-filtering approach[END_REF] versus √ 2-times). A drawback is that the DR filter is different for each projection angle, but it can be pre-computed for a given image grid and set of projection angles. Our results for proton CT data (Figure 5) have confirmed that the DR filter can improve spatial resolution compared with an FBP algorithm. This improvement was observed only for inserts close to the edge of the spiral phantom because at the center of the object, the image is too blurry due to Coulomb scattering to see any enhancement of spatial resolution.

While this method is theoretically demonstrated for parallel geometries only, the acquisitions can be done in another setup as long as the data is rebinned into parallel geometry, which is straightforward for list-mode pCT data, as was done for our own pCT simulations. It may be possible to derive an expression for the ramp filter in the case of fan beam data, for example, however, this might prove to be more difficult as the transform from the sinogram domain to the image domain is not a simple rotation as for parallel data.

Another possible use of the DR filter is for 3D parallel projections in xCT, when the detector is rotated by an in-plane angle. The projections would need to be realigned with the image grid before filtering. To avoid this interpolation, the 2D DR filter, computed for the in-plane angle direction, could be used on each 2D projection before backprojection.

CONCLUSION

We have proposed a 2D version of the ramp filter and applied it to simulated xCT and pCT data. This method starts by the smearing of each projection in the image space, then the directiondependent filtering, and finally the angular sum. The filter allows an accurate reconstruction and can be used to avoid an interpolation which would degrade spatial resolution.

Figure 1 :

 1 Figure 1: Parallel xCT geometry (left), pCT single tracking setup (middle), and TOF PET setup (right).

Figure 2 :

 2 Figure2: Smeared projection for angle φ = 0 (left) and φ = -π/4 (right). The blue circle is the field of view and the orange square is the reconstruction region. The shaded areas represent the regions where the image is truncated along θ. The green square is the region needed to accurately filter the smeared projection for all angles. For φ = 0, the data inside the orange square is sufficient for accurate reconstruction inside the FOV, while a larger region is needed for oblique angles, with the largest area being for φ = π/4 + kπ/2.

Figure 3 :

 3 Figure 3: Top : Pseudo-continuous impulse response of the 2D DR filter for φ = 0 (left) and φ = π/3 (right). The discrete samples are overlayed as black dots. Bottom : Fourier transform of the discrete impulse response.

FBP

Figure 4 :

 4 Figure 4: Shepp-Logan phantom reconstructed using, from left to right, an FBP, the DR filter sampled in the frequency and spatial domain, the oversampled DR filter (N=1024, in spatial domain), and the filter apodized with a Hamming window. The first two rows correspond to the results (reconstruction and difference with reference) for idealised projections and the two bottom rows to the results for noisy projections.

Figure 5 :

 5 Figure 5: (a) Proton CT reconstructions of the relative stopping power (RSP) map of the spiral phantom using an FBP (left) and the DR filter (right). (b) Spatial resolution of the inserts in the spiral phantom. The shaded areas correspond to ±3σ where σ is the uncertainty on the spatial resolution.

Acknowledgements

The work of Feriel Khellaf was supported by the project DIC20161236452 of the Fondation pour la Recherche Médicale (FRM). The work of Nils Krah was partially funded by the FRM and partially by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 753370. This work was performed within the framework of the SIRIC LYriCAN INCa_INSERM_DGOS_12563 and the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). We thank Rolf Clackdoyle, Laurent Desbat and Michel Defrise for useful discussions.


F Khellaf
email: feriel.khellaf@creatis.insa-lyon.fr

N Krah

J M Létang

Charles-Antoine Collins-Fekete

S Rit

A comparison of direct reconstruction algorithms in proton computed tomography

Keywords: proton CT reconstruction, spatial resolution, direct algorithms

Several direct algorithms have been proposed to take into account the non-linear path of protons in the reconstruction of a proton CT (pCT) image. This paper presents a comparison between five of them, in terms of spatial resolution and relative stopping power (RSP) accuracy. Our comparison includes (1) a distance-driven algorithm extending the filtered backprojection to non-linear trajectories (DD), (2) an algorithm reconstructing a pCT image from optimized projections (ML), (3) a backproject-then-filter approach using a 2D cone filter (BTF), (4) a differentiated backprojection algorithm based on the inversion of the Hilbert transform (DBP), and (5) an algorithm using a 2D directional ramp filter (DR). We have simulated a single tracking pCT set-up using Geant4 through GATE, with a proton source and two position, direction and energy detectors upstream and downstream from the object. Tracker uncertainties were added on the position and direction measurements. A Catphan 528 phantom and a spiral phantom were simulated to measure the spatial resolution and a Gammex 467 phantom was used for the RSP accuracy. Each proton's trajectory was estimated using a most likely path (MLP) formalism. The spatial resolution was evaluated using the frequency corresponding to an MTF of 10% of its peak value and the RSP accuracy using the mean values in the inserts of the Gammex phantom. In terms of spatial resolution, it was shown that, for ideal trackers, the DR and BTF methods offer a slightly better resolution since each proton is directly binned in the image grid according to its MLP. However, all methods but the ML show comparable resolution when using realistic trackers. Regarding the RSP, three algorithms (DR, DD and BTF) show a mean relative error inside the inserts about 0.1%. As the DR and BTF methods are more computationally expensive, the DD -which allows the same spatial resolution in realistic conditions and the same accuracy -and the DBP -which has a fairly good accuracy (< 0.2%) and allows reconstruction from truncated data -can be used for a reduced reconstruction time.

INTRODUCTION

Proton computed tomography (pCT) has been developed to reduce particle range uncertainties in hadron therapy treatment planning [START_REF] Schaffner | The precision of proton range calculations in proton radiotherapy treatment planning: experimental verification of the relation between CT-HU and proton stopping power[END_REF][START_REF] Yang | Comprehensive analysis of proton range uncertainties related to patient stopping-power-ratio estimation using the stoichiometric calibration[END_REF][START_REF] Paganetti | Range uncertainties in proton therapy and the role of Monte Carlo simulations[END_REF]. While current methods are based on a calibration between X-ray CT numbers and proton relative stopping power (RSP), pCT allows for a direct reconstruction of the RSP map and therefore a reduction of range uncertainties [START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF]. Although pCT provides a better estimation of the RSP, it suffers from poor spatial resolution, compared to X-ray CT, due to the stochastic nature of the non-linear proton paths. Disregarding nuclear interactions, a proton trajectory is a random combination of many small angular deflections due to multiple Coulomb scattering (MCS). It has been shown that these trajectories can be estimated using a most likely path (MLP) formalism [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF][START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF], a cubic spline path [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF][START_REF] Fekete | Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging[END_REF] or a polynomial model [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF], and including these estimates in the reconstruction can substantially improve spatial resolution.

The inverse problem in pCT is similar to that in X-ray CT as it consists in reconstructing an image from line integrals. What makes it different are the non-linear integration lines representing the proton path estimates. Several reconstruction methods incorporating a proton path estimate have been proposed for pCT [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF][START_REF] Penfold | Total variation superiorization schemes in proton computed tomography image reconstruction[END_REF][START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF][START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF]Hansen et al., 2014a;Hansen et al., 2014b;[START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF][START_REF] Collins-Fekete | A maximum likelihood method for high resolution proton radiography/proton CT[END_REF][START_REF] Hansen | Fast reconstruction of low dose proton CT by sinogram interpolation[END_REF]. However, since there is no exact analytical solution to invert proton integration paths yet, all proposed algorithms are heuristic solutions. An often used direct reconstruction method is the distance-driven binning algorithm, which generalizes the filtered backprojection algorithm to non-linear paths. Other direct algorithms invert the backprojection and filtering steps as it would be more fitting to list-mode proton data. Alongside direct reconstructions, various iterative algebraic methods were proposed, from the classical ART [START_REF] Gordon | Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[END_REF] to total variation superiorization algorithms. Iterative algorithms have advantages, e.g. including priors to improve image reconstruction (Hansen et al., 2014a), however, their computational cost is generally an order of magnitude larger than direct methods. [START_REF] Hansen | Fast reconstruction of low dose proton CT by sinogram interpolation[END_REF] compared the distance-driven method and three iterative algorithms to reconstruct low dose pCT images. The authors conclude that the distance-driven algorithm allows to reconstruct high-quality images while being at least 3 and up to 22 times faster than iterative algorithms.

However, no comparison between direct reconstruction methods was reported in the literature. The purpose of this paper is to compare five different direct algorithms based on spatial resolution as well as RSP accuracy. Although some of the direct methods are quite close, they use different kinds of approximations which could impact spatial resolution and/or RSP accuracy. For instance, as stated previously, backproject-filter approaches fit the list-mode data more naturally and might give a better spatial resolution, but some of these methods also require the computation of correction terms that might degrade the RSP accuracy. By comparing these methods, this work aims to select the algorithm giving the best spatial resolution while keeping a high RSP accuracy (relative error about 0.1% to 0.5% [START_REF] Penfold | A more accurate reconstruction system matrix for quantitative proton computed tomography[END_REF][START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF]) with a reasonable computational cost.

MATERIALS AND METHODS

Reconstruction in proton CT

The aim of pCT reconstruction is to build an RSP map from the protons' energy loss. For a given proton i, the line integral of the RSP is related to the energy loss using

WEPL i ≡ Γ i RSP(x)dx ≈ E in i E out i dE S water (E) (1)
where Γ i ⊂ R 3 is the proton path, RSP(x) is the stopping power relative to water at position x ∈ R 3 , E in i and E out i are the entrance and exit proton energies, and S water (E) is the stopping power of water for energy E. It is assumed that the RSP is energy independent (less than 0.7% variations in the 80-300 MeV energy range for surrogate human tissues [START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF] for the equality to be valid. The integral corresponds to the water-equivalent path length (WEPL).

One of the proposed set-ups for pCT tracks each proton individually and measures its position and direction before and after traversing the object, as well as its energy loss [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF]. It is possible to get each proton's WEPL from the energy measurement using the right-hand side of Equation 1. Its MLP is computed from the position and angle information. The pCT reconstruction problem consists in using each proton's WEPL and MLP to recover the RSP map.

Figure 1: pCT set-up with a parallel proton beam. The left figure corresponds to the geometry in the frame of reference of the source/detector, i.e. the unit vectors (e u , e w ) define a coordinate system relative to the projection angle, with e w the proton beam direction and e u defining the detector surface. The right figure corresponds to the geometry in the frame of reference of the object, with (e x , e y ) a fixed coordinate system. The MLP is computed for each proton i at each w k giving u i,k and v i,k (the axial coordinate is not represented on this 2D slice). x i,k and y i,k correspond to a rotation of w k and u i,k .

Reconstruction algorithms

The following sections summarize the five algorithms included in this comparison. All algorithms are presented in a parallel beam set-up for simplicity, although they can all be used with other geometries.

Distance-driven algorithm (DD)

It is the first direct reconstruction algorithm incorporating the MLP that was proposed in the literature [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF]. It consists in binning the protons into projections on a series of planes parallel to the detectors (Figure 1) according to the estimated proton path Γi (t) ∈ R 3 at time t ∈ R, rather than binning just in the detector plane. The projection binning is done at several distances, hence the name of the algorithm.

The MLP formalism yields the coordinates u i (w) and v i (w) for each proton i (Figure 1). In practice, we compute this trajectory at a discrete number of positions k ∈ {1, ..., K} such that u i,k = u i (w k ) and v i,k = v i (w k ) with w k the discrete distances. The 3D proton trajectory at the time when proton i crosses the plane parallel to the detectors at distance w k is written Γi (t i,k ) ∈ R 3 . Without loss of generality, we use the 3D path to bin the list-mode data into projections, but only the 2D central slice is considered in the following. The value in the sinogram for pixel j and source position p is given by

g j,p = i∈Ip k ζ j (u i,k , v i,k , w k )WEPL i i∈Ip k ζ j (u i,k , v i,k , w k ) (2)
where ζ j is the indicator function for pixel j defined as

ζ j (y) = 1 if y ∈ R 3 is in pixel j, 0 otherwise, (3)
and I p is the subset of protons emitted from the same source position p. Note that the sinogram for a single slice is in 3D since it depends on the source position and the 2D pixel index.

The filtering and backprojection steps are the same as in a standard filtered backprojection algorithm, except that a voxel-specific backprojection is used to select which distance w to use, instead of the usual pixel-specific backprojection. Specifically, we define the 2D projections g p (u, w) : R 2 → R as g p (u, w) = j g j,p ψ j (u, w), (4) with ψ j (u, w) a given interpolation function. The reconstructed image f DD (x) : R 2 → R is then given by

f DD (x) = p gp (x • e u , x • e w )∆θ. ( 5
)
with x = (x, y) ∈ R 2 , gp the ramp-filtered projection g p in the case of a parallel beam geometry, and ∆θ the constant angular spacing between consecutive projections.

Maximum likelihood method (ML)

This method proposes to maximize the likelihood of the water equivalent thickness (WET) in proton radiographs, i.e. on the projection level, in order to improve their spatial resolution [START_REF] Collins-Fekete | A maximum likelihood method for high resolution proton radiography/proton CT[END_REF]. These optimized radiographs are then used to reconstruct a pCT image. Specifically, the WET's likelihood is maximized for each so-called "channel" c in a radiograph, where a channel is the volume defined by the pixel corners along the source direction (see Figure 1). The WET value for channel c and source position p is obtained given the measured WEPL and the estimated trajectory for each proton using

WET c,p = i∈Ip l 2 c,i L 2 i WEPL i i∈Ip l 2 c,i L 2 i (6)
with l c,i the length of the i-th proton in channel c and L i its total path length. This is similar to the DD algorithm, except that the WEPL is integrated along the w direction to form a single projection per source position while the DD uses a series of depth dependent projections. Also, the weights for the ML correspond to proportions of proton trajectories spent in a given channel

( l 2 c,i L 2 i
) while for the DD method, the weighting function ζ j is binary. To reconstruct the pCT image, the radiographs WET c,p are simply fed to a standard filtered backprojection algorithm [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF]. As in the DD, we define interpolated projections

WET p (u) = c WET c,p ψ c (u), (7)
giving the following reconstruction

f ML (x) = p WET p (x • e u )∆θ, (8)
where WET p (u) is the filtered projection WET p (u).

Backprojection-then-filtering (BTF)

This algorithm inverts the order of the backprojection and the filtering so that we directly backproject the protons along their MLPs, instead of binning them into a sinogram first (Polud-niowski et al., 2014). An intermediate step for calculating the backprojection from the list-mode data is a pixel-wise and direction-wise binning,

b m,n = i,k ζ m ( Γi (t i,k ))ξ n (Θ i,k )WEPL i i,k ζ m ( Γi (t i,k ))ξ n (Θ i,k ) , (9)
with indicator functions

ζ m ( Γi (t i,k )) = 1 if Γi (t i,k ) ∈ R 3 is in pixel m, 0 otherwise, ( 10
)
and

ξ n (Θ i,k ) = 1 if Θ i,k ∈ R is in angular bin n. 0 otherwise. ( 11
)
This binning was first presented by [START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF]. The WEPL for each proton is smeared along the corresponding path in the image grid using ζ m , an indicator function for pixel m, and ξ n , an indicator function for angular bin θ n . The angle Θ i,k is that of proton i at depth w k with respect to e x (Figure 1). This binning was also used for the algorithms described in sections 2.2.4 and 2.2.5. Note that b m,n is not strictly speaking a backprojection but rather an average of the WEPL of protons passing through pixel m with an angle θ n . This direction-wise binning allows to rebin the data into parallel geometry (in case another geometry is used) as it uses the proton's angle Θ i,k with respect to e x rather than the subsets I p as in the previous methods. We interpolate an image b n (x) : R 2 → R defined as

b n (x) = m b m,n ψ m (x), (12)
giving the backprojected image

B(x) = n b n (x)∆θ. ( 13
)
In the BTF method [START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF], the image is reconstructed by convolving the backprojection with a 2D cone filter c(x), yielding

f BTF (x) = (B * c)(x). ( 14
)
The expression for the filter c(x) is given by [START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF] in Equations ( 12)-( 14).

It was derived by taking the inverse Fourier transform of the cone filter ρ = (ν x , ν y ) where (ν x , ν y ) are the 2D Fourier variables.

The backprojection B has an infinite support. In order to mitigate the effects of truncating B during the convolution, the backprojection is computed on a very large matrix (4 to 8 times the size of the reconstructed image). We chose to use a factor of 4 to keep a reasonable computation time. The contribution of distant data of B to the central pixel is approximated by a factor ∆f (Equation 19 [START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojectionthen-filtering approach[END_REF]) obtained by replacing B = f BTF * 1/r in Equation 14when r → ∞, with r = (x, y) . This factor is used to correct for the constant bias (i.e. pixel independent) induced by the truncation of the backprojection matrix, giving the final image

f BTF,∆ = f BTF + ∆f.
(15)

Directional Ramp filtering (DR)

This is an intermediate method between a backproject-first and a filter-first algorithm [START_REF] Khellaf | 2D directional ramp filter[END_REF]. It starts by ( 1) computing the pixel-wise and direction-wise binning (Equation 9), then (2) filtering each direction with the corresponding directional ramp filter and (3) summing all filtered backprojections. The reconstructed image is given by

f DR (x) = n bn (x)∆θ, ( 16
) where bn (x) = (b n * h n )(x). ( 17
)
The expression of h n is given by [START_REF] Khellaf | 2D directional ramp filter[END_REF]. The directional filter is a two-dimensional version of the usual ramp filter. As in backproject-first algorithms, this method starts by smearing the projections into the image space before filtering, but the sum over all projection angles is the last step, as in filter-first methods.

Differentiated backprojection (DBP)

This algorithm is based on the two-step Hilbert transform method [START_REF] Noo | A two-step Hilbert transform method for 2D image reconstruction[END_REF]. It relies on the relationship between the backprojection of the derivative of the projections and the Hilbert filtered image:

f (x, y) = H -1 π 0 ∂ ∂s p(s, θ) s=-x sin θ+y cos θ dθ, (18)
with H -1 the inverse of the Hilbert transform. [START_REF] Noo | A two-step Hilbert transform method for 2D image reconstruction[END_REF] give a formula to invert the finite Hilbert transform, which requires knowledge of the convex hull of the object (an information which is also used in our implementation of the MLP). Simply put, the image can be reconstructed by inverting the finite Hilbert transform of the differentiated backprojection. The advantage of this approach is that, since the derivative is a local operator, truncated data can still be reconstructed in cases when points are lying on at least one line segment having its two end points in the field of view and outside the object to reconstruct. This method was adapted to work with list-mode data by inverting the order of the differentiation and the backprojection steps [START_REF] Zeng | Image reconstruction via the finite Hilbert transform of the derivative of the backprojection[END_REF], giving the following reconstruction formula

f DBP (x, y) = H -1 ∂ ∂x B s (x, y) + ∂ ∂y B c (x, y) . (19)
Due to the partial differentiation, two intermediate weighted backprojections must be computed

B s (x) = - n b n (x) sin θ n ∆θ, ( 20
) B c (x) = n b n (x) cos θ n ∆θ. ( 21
)
This method was applied to pCT data by [START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF].

Simulations

Monte Carlo simulations were conducted using Geant4 [START_REF] Agostinelli | GEANT4 a simulation toolkit[END_REF] through Gate [START_REF] Jan | GATE V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy[END_REF] to acquire pCT list-mode data. A "single tracking" setup was used, consisting of two ideal detectors measuring position, direction and energy, located upstream and downstream from the phantoms. A 200 MeV rotating fan-beam proton source was positioned at 1 m distance from the object center, with a flux of 225 protons•mm -2 •projection -1 at the isocenter, for a total of 720 projections over a 360 degree range. Protons were paired between the two detectors, and each proton's path was estimated using the MLP formalism, after excluding protons with nuclear interactions using 3σ angle and energy cuts [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. Three phantoms were imaged: the CTP528 module of the Catphan phantom (The Phantom Laboratory, NY) with aluminium line pairs of different resolutions in a 20 cm diameter water cylinder; a spiral phantom, consisting of aluminium cylinders placed along a spiral in a 20 cm diameter water cylinder [START_REF] Rit | List-mode proton CT reconstruction using their most likely paths via the finite Hilbert transform of the derivative of the backprojection[END_REF], and the Gammex 467 tissue characterization phantom with different tissues inserts inside a 33 cm diameter water cylinder (Figure 2). The phantoms used for spatial resolution assessment (Catphan and spiral) were reconstructed in images of 1000 × 1000 pixels with a pixel spacing of 0.25 mm. The Gammex 467 was reconstructed on a grid of 800 × 800 pixels of spacing 0.5 mm and required a beam energy of 250 MeV because of its larger diameter.

In order to estimate the impact of detector uncertainties on spatial resolution, an uncertainty on the simulated position and direction data was added using the formalism of [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF]. The parameters were chosen to reflect a typical pCT prototype. In particular, the detector resolution (pitch of strips) was characterized by a Gaussian of σ p = 0.15 mm, the distance between the trackers was set to d T = 10 cm (each real detector consists of a pair of trackers), the distance between the trackers and the object border to 20 cm, and the material budget to x/X 0 = 5 × 10 -3 .

The energy resolution of the detectors and the spread of the incident beam were omitted for simplicity as scintillating detectors [START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF] are capable of achieving a WEPL resolution close to the limit due to proton range straggling (measured WEPL resolution of 3 mm for a range straggling of 2.85 mm, for 200 MeV protons in water), although we note the beam spread can contribute to the WEPL variance up to 20% [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF]. Electromagnetic interactions are taken into account in the simulation and are responsible for the RSP resolution of the reconstructed images.

Image quality metrics

Both spatial resolution and RSP accuracy were evaluated. The spatial resolution was assessed qualitatively using the Catphan phantom and quantitatively with the spiral phantom. The modulation transfer function (MTF) was calculated in the spiral phantom using each aluminium cylinder's radial edge-spread-function (ESF), as was done by [START_REF] Richard | Towards task-based assessment of CT performance: system and object MTF across different reconstruction algorithms[END_REF]). An error function with parameters µ and σ was fitted to sampled points of the ESF, assuming

ESF(x) ∝ 1 2 1 + erf x -µ σ . ( 22
)
The frequency corresponding to an MTF value of 10%, given by

MTF 10% = ln 10 2 1 πσ , (23)
was used to quantify the spatial resolution. The error on MTF 10% was derived from the standard deviation of the error function fit. In addition, we measured the radial and azimuthal resolutions defined as defined by [START_REF] Plautz | An evaluation of spatial resolution of a prototype proton CT scanner[END_REF] in order to assess the (an)isotropy of the spatial resolution. The radial resolution was measured using ESFs along the direction from the isocenter to the center of the bead and the azimuthal resolution was measured from profiles along the direction perpendicular to the radial direction. Profiles were sampled and averaged for angles in a range of ±10 • around these two directions. Each algorithm's RSP accuracy was evaluated by computing the mean value in each insert of the Gammex phantom. A region of interest (ROI) of radius 8 mm was used (the radius of the inserts is 14 mm) so that the measured RSP value would not be impacted by the inserts' spatial resolution. The uncertainty on the measured RSP is computed using the 95% confidence limits of the mean, i.e. ±1.96σ ROI / √ N ROI , with σ ROI the standard deviation of the RSP values in the ROI and N ROI the number of pixels inside the ROI.

RESULTS

The reconstructed images of the Catphan and spiral phantoms are shown in Figure 3. For the Catphan phantom, the five reconstructions are quite similar, as confirmed by the line profiles drawn in Figure 4. All algorithms underestimate the RSP of aluminium in the Catphan phantom due to the impact of axial resolution. The ML reconstruction leads to a lower resolution, causing a larger underestimation of the RSP value of the aluminium inserts. Regarding the spiral phantom reconstructed from data simulated with ideal trackers (middle column of Figure 3), the central beads appear much blurrier than the outer beads for all algorithms. This is expected because the MLP uncertainty is largest for those inserts. A slight artifact is seen for the ML reconstruction where the outer beads appear smeared along the azimuthal direction. This is visible on the zoomed insets, although it is exaggerated by using a color scale with a small range. Figure 5 shows the spatial resolution as a function of each bead's distance to the centre of the spiral phantom, for reconstructions with ideal trackers and realistic trackers. For ideal trackers, resolution is about a factor 4 to 6 higher for beads near the phantom borders compared with the resolution at the centre. The resolution of the ML algorithm is lower than other algorithms: it only reaches 1 lp/mm for the outermost bead and is also smaller for central beads. Differences in spatial resolution among the other four algorithms are much smaller, especially below a distance of about 80 mm from the center. Beyond that distance, i.e. closer to the phantom edge, the DR and BTF algorithms reach similar resolutions (up to 3.5 and 3.3 lp/mm, respectively), while the DBP and DD methods both reach about 2.7 lp/mm. Figure 6 presents a more extensive analysis of the ideal spatial resolution, discriminating between the radial and azimuthal resolutions. The azimuthal resolution is higher than the radial resolution for all algorithms except ML. The right column of Figure 3 shows that the resolution with realistic trackers is considerably lower than with ideal ones. In Figure 5, we observe that the maximum resolution for all algorithms only reaches between 0.39 and 0.42 lp/mm, and that the difference of spatial resolution between the center and the borders is much less important than with ideal detectors. For example, for the DD algorithm, it only increases by a factor of 1.3 between the central bead and the outer bead, versus a factor of 4.6 with ideal detectors. In addition, the differences in spatial resolution observed between the algorithms are mostly insignificant, considering the uncertainty on the MTF. Only the ML algorithm still shows a slightly lower resolution although it is much closer to other algorithms compared with to ideal case.

Table 1 shows the relative error on the RSP value in each insert of the reconstructed Gammex phantom. The most accurate algorithms are the DR and DD, with a mean relative error of 0.07% and 0.08% respectively, followed by the BTF with an error of 0.11%, the DBP with 0.19%, and the ML with a mean error of 0.25%. The uncertainty in Figure 7 shows that the majority of the reconstructed RSP values are very close, with differences within the confidence intervals of the mean. The RSP error in the lung inserts, however, reaches over 1% for the ML method.

DISCUSSION

The aim of this work was to compare the performance of five direct reconstruction algorithms in terms of spatial resolution and RSP accuracy. Spatial resolution in pCT is impacted by various factors. A first factor is MCS which is mitigated to some extent by modeling the MLP [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. This greatly improves spatial resolution compared to naive straight line trajectories, although the effect of scattering is still important especially at the center of the phantoms where the MLP estimate is most uncertain. We have observed that the resolution at the center of a 20 cm wide phantom is about 0.5 -0.6 lp/mm with ideal trackers. For the DR method, the resolution is improved by a factor of 6 between the center and the border of the spiral phantom. Since the center of the object corresponds to the largest uncertainty on the MLP, differences in spatial resolution between the algorithms are not significant near the center. Another consequence of proton scatter is the anisotropy of the spatial resolution. As expected, the azimuthal resolution gets higher than the radial resolution with the distance to the center. Due to scatter, the error envelope of the MLP of protons traversing the same region will depend on the thickness of material traversed, i.e. on their initial direction. For central beads, the Figure 7: Relative error on RSP with confidence limits.

variation of the traversed thickness depending on direction will be smaller, such that resolution is isotropic at the center. The only algorithm for which the azimuthal resolution is poor is the ML algorithm, due to the artifact observed in Figure 3. This artifact can be explained by the uncertainty on the MLP. We mentioned that the ML method basically consists in integrating the DD's binning along the w direction. We know that the projections that are farthest from the detectors along the w direction will be the most blurry due to a large uncertainty on the MLP. For the outer beads of the spiral phantom, the most blurry projection contributing to their reconstruction is the one where protons have had the longest trajectory inside the phantom, i.e. the ones that contribute to the azimuthal resolution. This explains why the outer beads seem smeared on the ML reconstruction. Another factor of resolution loss is the tracker uncertainties which have a significant impact as can be seen in Figure 5. Adding an uncertainty on the position and direction measurements dramatically lowers the spatial resolution compared to ideal simulated data. It should be noted that those results were obtained for a distance of 20 cm between the trackers and the borders of the phantom. Using larger distances would further reduce resolution.

Finally, the choice of the algorithm also has an impact on the resolution as different binnings and interpolations are involved. For example, using ideal detectors, the resolution 96 mm away from the center for the DR method is 250% higher than the ML algorithm's resolution (only 25% at the center), and 26% higher than the DD algorithm's resolution (1% at the center). Our results show that the best performing algorithms in terms of spatial resolution away from the center are the DR and BTF, suggesting that backprojecting the protons before filtering can improve spatial resolution. Indeed, filter-first methods imply two interpolations that can reduce spatial resolution: the projection binning and the backprojection. On the other hand, only one is needed for backproject-first methods. This could preserve the spatial information of proton paths and therefore the resolution. However, we have seen that the DBP method, which starts with the backprojection, has the same spatial resolution as the DD. A possible explanation is that the discrete differentiation computed in the DBP algorithm, using differences of adjacent pixels, may cause a loss of spatial resolution. The ML algorithm shows the worst spatial resolution as it attempts to account for MLPs in the 2D projections, but they can only be described in 3D space. This method might be better suited for proton radiography than CT. Overall, the impact of the choice of the algorithm is negligible compared with MLP and tracker uncertainty.

In terms of RSP, the DR, DD and BTF show good accuracy with mean absolute errors close to 0.1%. The DBP method is slightly less accurate, which can be linked to the computation of a constant during the Hilbert transform inversion [START_REF] Noo | A two-step Hilbert transform method for 2D image reconstruction[END_REF]. However, this is not significant with respect to the 95% confidence interval of the mean (Figure 7). The ML algorithm overestimates the RSP value in the lung inserts compared with the other methods, resulting in a higher mean error. An explanation is that as it uses channels instead of pixels, the WEPL value of a proton that has not traversed a lung insert can still be attributed to it. While this would happen for all inserts, lung inserts would be more impacted because the difference between the RSP of the lung and the rest of the phantom is larger than for other inserts.

The computational cost of each algorithm is an important factor when several algorithms reach similar levels of resolution and RSP accuracy. The binning of the list-mode data is the most costly operation as the MLP of each proton has to be computed. The binning steps to compute projections as in Equation 2 or smeared projections as in Equation 9 are mostly similar in terms of computational cost. This cost increases with the number of protons and the number of evaluations of the MLP. As we have used a high fluence (288×10 6 protons in total) and a small spacing for the MLP evaluation (0.25 mm or 0.5 mm depending on the phantom), our binning is orders of magnitude slower than the other reconstruction steps (for the DD, 8 s/projection × 720 projections for the binning, and 42s for the filtering and backprojection). Both the BPF and DR methods require the backprojection region to be larger than the reconstructed image region, increasing the computation time of the binning. While the DR method requires backprojecting the list-mode data in a region larger by a factor of 2, the BPF method needs a matrix of at least four times the reconstruction matrix. They are therefore the most computationally expensive methods. The cost of the filtering step varies a lot depending on the algorithm. For example, for the DD, a 3D sinogram needs to be filtered (for each slice) compared with the usual 2D sinogram for the ML. However, as we stated above, the cost of the filtering step is quite negligible compared with the cost of the binning step.

CONCLUSION

We have studied the spatial resolution and RSP accuracy of five different algorithms used to reconstruct proton list-mode data. The results show that, using ideal trackers, two of the methods that bin the list-mode data directly in the image space (DR and BTF) offer a better spatial resolution at the borders of the object, with a maximum of 3.5 lp/mm at the borders of a 20 cm thick object. However, considering realistic measurement uncertainties, the impact of the choice of the algorithm on spatial resolution becomes negligible. Furthermore, the DR, DD, and BTF methods show the best RSP accuracy with a mean error close to 0.1%. While the DR and BTF methods have a slight advantage regarding spatial resolution as well as a good accuracy, they are also the most computationally expensive. For faster reconstructions, the DD algorithm offers the same accuracy and an equivalent spatial resolution considering detector uncertainties, and the DBP, while being somewhat less accurate, can be used for reconstruction with truncated data. The ML, as it works in the projection level, is well suited for proton radiography.

Figure 2 :

 2 Figure 2: RSP map of the Gammex 467 tissue characterization phantom.

Figure 3 :

 3 Figure 3: Reconstructions of the Catphan and spiral phantoms for the five algorithms. The Catphan is reconstructed from idealized data. Zoomed insets of a peripheral bead are shown to appreciate the artifact reconstructed by the ML algorithm. The extent of the color scale in the insets is [0.7,1.3] to highlight the artifact. 10

Figure 4 :

 4 Figure 4: Line profiles along first (top) and third (bottom) line pairs in the Catphan phantom. The dotted lines represent the theoretical RSP values of water (RSP=1) and aluminium (RSP=2.1).

Figure 5 :

 5 Figure 5: Spatial resolution measured as the frequency for an MTF value of 10% in function of distance from the isocenter, for ideal trackers (top) and realistic trackers (bottom). The shaded areas represent the uncertainty on the MTF estimation.

Figure 6 :

 6 Figure 6: Spatial resolution for ideal trackers computed on profiles along the normal and tangent directions to the edge. Error bars were omitted for clarity.

Table 1 :

 1 Relative error on the measured RSP values.

	Insert	Reference DD	BTF	DBP	ML	DR
	LN300(1)	0.295	-0.15% 0.02%	-0.44% 1.32%	-0.21%
	IB3	1.094	-0.04% 0.02%	-0.19% -0.02% -0.03%
	LV1 (1)	1.084	-0.03% 0.05%	-0.21% -0.03% -0.01%
	LN450	0.443	-0.27% -0.06%	-0.24% 0.54%	-0.21%
	B200	1.103	-0.04% 0.06%	-0.18% -0.17% <0.01%
	Water Solid (1) 1.005	-0.03% 0.06%	-0.16% 0.05%	-0.01%
	CB2 30	1.286	0.04%	0.10%	-0.14% -0.13% 0.07%
	AP6 (1)	0.940	0.04%	0.08%	-0.15% 0.03%	0.02%
	BR12	0.974	-0.01% -0.17%	-0.15% -0.03% 0.01%
	BRN-SR2	1.073	0.08%	-0.06%	-0.05% 0.01%	0.07%
	CB2 50	1.449	-0.11% -0.18%	-0.25% -0.13% -0.08%
	SB3	1.649	-0.16% -0.19%	-0.20% -0.13% -0.08%
	LN300 (2)	0.295	0.05%	-0.59%	-0.30% 1.16%	-0.07%
	AP6 (2)	0.940	0.14%	<0.01% 0.19%	0.18%	0.18%
	LV1 (2)	1.084	0.04%	-0.08%	-0.06% 0.06%	0.07%
	Water Solid (2) 1.005	0.06%	-0.06%	-0.11% 0.06%	0.07%
	Mean abs error		0.08%	0.11%	0.19%	0.25%	0.07%

Acknowledgements

The work of Feriel Khellaf was partially supported by the project DIC20161236452 of the Fondation pour la Recherche Médicale (FRM). The work of Nils Krah was supported by funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 753370. This work was performed within the framework of the SIRIC LYriCAN INCa INSERM DGOS 12563 and the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).