Integrated-mode proton radiography with 2D lateral projections

Integrated-mode proton radiography leading to water equivalent thickness (WET) maps is an avenue of interest for motion management, patient positioning, and in vivo range verification. Radiographs can be obtained using a pencil beam scanning setup with a large 3D monolithic scintillator coupled with optical cameras. Established reconstruction methods either (1) involve a camera at the distal end of the scintillator, or (2) use a lateral view camera as a range telescope. Both approaches lead to limited image quality. The purpose of this work is to propose a third, novel reconstruction framework that exploits the 2D information provided by two lateral view cameras, to improve image quality achievable using lateral views. The three methods are first compared in a simulated Geant4 Monte Carlo framework using an extended cardiac torso (XCAT) phantom and a slanted edge. The proposed method with 2D lateral views is also compared with the range telescope approach using experimental data acquired with a plastic volumetric scintillator. Scanned phantoms include a Las Vegas (contrast), 9 tissue-substitute inserts (WET accuracy), and a paediatric head phantom. Resolution increases from 0.24 (distal) to 0.33 lp mm−1 (proposed method) on the simulated slanted edge phantom, and the mean absolute error on WET maps of the XCAT phantom is reduced from 3.4 to 2.7 mm with the same methods. Experimental data from the proposed 2D lateral views indicate a 36% increase in contrast relative to the range telescope method. High WET accuracy is obtained, with a mean absolute error of 0.4 mm over 9 inserts. Results are presented for various pencil beam spacing ranging from 2 to 6 mm. This work illustrates that high quality proton radiographs can be obtained with clinical beam settings and the proposed reconstruction framework with 2D lateral views, with potential applications in adaptive proton therapy.


Introduction
Proton computed tomography is a developing imaging modality that estimates the spatial distribution of the relative proton stopping power (RSP) using tomographic principles.Alternatively, radiographic images representing the water equivalent thickness (WET) of traversed tissues can be acquired.Proton imaging may be helpful for accurate treatment planning in proton therapy by limiting uncertainties in RSP estimation (Schneider and Pedroni 1995, Schaffner and Pedroni 1998, Schneider et al 2005, Paganetti 2012, Yang et al 2012, Dedes et al 2022).Furthermore, as the treatment and imaging sources are the same, proton imaging creates images that are fully registered with the treatment beam, and enables a direct measurement of the residual proton range in the treatment room, right before delivery.As such, proton imaging may also be useful for patient positioning and adaptive proton therapy (Parodi 2020).Specific use of proton radiographs for adaptive proton therapy include identifying the time for adaptive re-planning in head and neck cancer patients (Fukumitsu et al 2014, Wu et al 2017, Evans et al 2020, Meijers et al 2021), or rapid radiographic imaging of moving tumours towards online adaptive radiotherapy (Steinsberger et al 2022).
Single-event proton imaging (Johnson 2017) is the most promising approach in terms of spatial resolution and WET/RSP accuracy (Krah et al 2018b, Parodi 2020), although practically such systems are hindered by their high costs, slow imaging time and limited compatibility with proton fluxes commissioned for clinical systems (Parodi 2020).Integrated mode imaging, which considers the integrated signal from individual pencil beams (PB) rather than single protons, is a faster alternative generally compatible with clinical systems.While multiple integrated-mode detection systems have been proposed (Zygmanski et al 2000, Testa et al 2013, Rinaldi et al 2014, Farace et al 2016, Darne et al 2017, 2019, Meijers et al 2021, Tendler et al 2021, Darne et al 2022, Schnürle et al 2023), their image quality remains limited compared to single event imaging, mostly due to multiple Coulomb scattering (MCS) blurring the integrated signal, and range mixing effects.Improving image quality for integrated mode imaging therefore necessitates advanced image reconstruction techniques (Parodi 2020).While there have been some efforts to enhance image quality via deep learning models (Heyden et al 2021), this work explores the use of a novel physics-based reconstruction model to refine image quality.
More specifically, the purpose of this work is to present a novel proton radiography reconstruction algorithm for a fast, low-cost integrated-mode proton imaging device.The device is based on a volumetric plastic scintillator that creates a 3D light emission distribution from the dose distribution of individual PBs inside the scintillation detector.2D projections of this signal are captured using optical apparatus such as CCD cameras (Darne et al 2017, 2019, Tendler et al 2021).The conventional setup (Tanaka et al 2016, Darne et al 2017, Tanaka et al 2018, Darne et al 2019) typically uses one CCD camera at the distal end of the scintillator, which integrates the optical signal along the beam's axis.This beam's eye view signal allows fast and straightforward reconstructions, but image quality largely suffers from MCS, and the approach requires an extensive calibration procedure to be quantitative (Darne et al 2019).
The setup can however be used as a range telescope, provided that the CCD camera is oriented to capture a lateral projection perpendicular to the beam's propagation axis; an early investigation has shown promise towards improved image quality (Tendler et al 2021).However, in that approach, the lateral projections did not exploit the beam's spatial position, which could be inferred as lateral projections result in 2D images of the Bragg curve.In addition, there was no consideration of the impact from MCS in the reconstruction process.This limits the reconstruction accuracy and does not fully take advantage of the information present in each 2D image.This study presents a reconstruction algorithm tailored specifically for the problem of WET reconstruction from a set of 2D orthogonal lateral projections which incorporates MCS physics into the reprojection, and considers range mixing effects.
A probabilistic framework for image reconstruction is first presented.Reconstructed radiographs from Monte Carlo simulations on an extended cardiac torso (XCAT) phantom as well as a slanted edge phantom are first reported to illustrate differences in image quality and WET accuracy between the proposed integrated mode method as well as conventional integrated mode distal and lateral views.Integrated mode approaches are also compared with single event proton imaging.Experimental data to reconstruct proton radiographs was also gathered for a contrast phantom (Las Vegas), 9 tissue-equivalent inserts for WET accuracy, and a paediatric head phantom for general image quality.

Methods
2.1.Pencil beam propagation and detection 2.1.1.Imaging setup A pencil beam (PB) scanning setup is considered for proton imaging.Let i be the index noting individual pencil beams PB i , with )and propagates towards the z direction, with initial angle / / q q q q p p = = , , 2, 2, 0 ) .Each PB is assumed Gaussian, and has the same initial spatial and angular spread in the lateral directions x and y, defined in the covariance matrix: In the x direction, the initial spatial spread of the beam is s z x i ( ).The initial angular spread s q z i x ( )is the mean square angular deflection of the beam in x, and the spatial-angular covariance s q z x i , x ( )is typically obtained from a beam's measured emittance (Gottschalk 2012).For proton radiography, the aim is to reconstruct the WET of the object.The reconstruction grid has isotropic voxels of side d 2 k and the position of pixels on the grid are noted with = x y z r , , ,where z k is defined from the distance between the source and the object's entrance, z O .= x y z r , , )is the location of a detection element, with z d defined with respect to the detector entrance, which is at distance z D from z i .The detector considered in this work produces a 2D signal.Figure 1 schematizes the imaging setup for a single PB, while the detection geometry is further detailed in figure 2. A 2D detector system is considered in this work, following the setup of Darne et al (2022).The setup involves a monolithic scintillator scoring 3D energy deposition.A CCD camera or any pixelised optical detector whose field of view covers the entire scintillator is positioned at angle f to gather 2D projections of the 3D signal from the scintillator.In this work, images are acquired in two perpendicular lateral (XZ and YZ) planes, resulting in data acquired respectively at camera angles f = 0 and π/2. Figure 2 shows the imaging setup for the XZ projection, examples of 2D datasets for the XZ and YZ views, and an example of experimental setup.

Probabilistic framework
Each PB i generates a signal in the detector at various locations r d .To construct a WET map, the signal at relevant detector voxels r d is reprojected towards the reconstruction grid, as detailed in section 2.3.Reprojection requires the calculation of f  r r r , ), the likelihood that a PB originating from r i and detected at r d (with a detection setup using camera angle f) crosses voxel r k : The right hand side includes  r r k i ( | ), the probability that the PB crosses voxel r k given that it originates from r i .
( | ) is defined similarly for detector location r d , and camera angle f. f  r r r , )is the probability that the PB is  detected at r d in the detector (using camera angle f) knowing that it originates at r i and passes through the reconstruction voxel located at r k .

Likelihood model based on multiple scattering theory
Relevant results from the Fermi-Eyges theory of multiple scattering (Fermi 1940, Eyges 1948)  ( | ), the likelihood that PB i originating from r i crosses a reconstruction grid voxel at r k , is obtained by marginalizing the beam's 3D localisation probability density function (PDF), p r r i ( | ), over the reconstruction voxel's dimensions: For the integral over z, s z l k ( )is assumed constant within a given voxel (Rescigno et al 2015).The calculation of s z l k ( ), the spatial spread of the beam after propagating from z i to z k , is detailed in appendix B.1.

Transport from the source to the detector
( | ) over relevant detector dimensions, which could be 1D (range telescope), 2D (this work) or 3D (voxelised detector).Considering the 2D XZ projection illustrated in figure 2, the likelihood is 4), but with y d and y i instead of x d and x i .s z l d ( ) is the spread of the beam from z i to z d in the detector, and its calculation is detailed in appendix B.2.

Transport from the object to the detector
), the spatio-angular conditional PDF of protons in the XZ plane, is derived; this PDF is first introduced in equation (A1).Second, this distribution is marginalized over all angles q x to obtain the location conditional PDF p x z r r , , A similar procedure is performed for the y coordinate to obtain p y z r r , , ), and both lateral PDFs are multiplied to obtain, as a third step, the location conditional PDF p r r r , )is interpreted as the spatio-angular PDF of a PB coming from r i and passes through the voxel at location r k .The condition of passing through r k is viewed as the beam being collimated by a virtual 2D collimator at the reconstruction voxel's location between lateral coordinates The calculation of the location PDF for a collimated PB is based on the Fermi-Eyges PB summation method, whose details can be found in the work of Sabbas et al (1987) and Safai et al (2008): Here ¢ x and q ¢ x are integration variables denoting the position and angle of protons at the collimator plane.
is the PDF of a needle-like Gaussian beam propagating from z k to z, which assumes initial values of S = z 0 x k ( ) , and is calculated with equation )is the distribution of protons at the collimator plane, also obtained with equation (A1), with elements of S z x ( ) calculated in the same way as in section 2.2.1.Furthermore, integration over the collimator dimensions in ¢ In equation (6), N x is a normalization constant to obtain probabilities, and the following variables are introduced: Here, the elements of the covariance matrix s * z k ( ) are obtained after propagation from z i to the object voxel's location z k , and  z 2 ( ) is the multiple scattering contribution to spatial spread (equation (A6)) from the object voxel's location z k to an arbitrary depth z.Details to calculate s * z k ( ) and  z 2 ( ) are provided in appendix B.3.The joint PDF p r r r ,  ,d .8 The integral over y equals 1, as it is assumed that the y component of p r r r ) is fully contained within the detector.As there is no analytical solution to equation (8), f=  r r r , 2.3.Reprojection from joint lateral views 2.3.1.Joint likelihood Consider the specific case where two lateral views with camera angles f = 0 and f = π/2 (XZ and YZ planes) are obtained.Assuming that each 2D image is composed of one or more identifiable Bragg peaks that can be associated with various structures traversed by the beam, it is possible to match signals from both 2D views and infer the 3D position of the Bragg peaks.In other words, the x z , ) , where z d is taken as the average position between the two views to account for experimental or algorithmic uncertainties.With the proposed experimental setup, the difference between the two views was found to vary by no more than one pixel.In this situation, the likelihood of equation (2) becomes the intersection of both views: r r r r r r r r r r r r , , , .9 Equation ( 9) can be obtained by considering the derivations of section 2.2 with a 3D detector element instead.

Reprojection equation
For a given detection geometry, each PB i generates a signal S r i d ( )at multiple r d .Projection data can be generated similarly to the distance-driven binning procedure suggested by Rescigno et al (2015), which is an adaptation of Rit et al ʼs distance-driven binning framework (Rit et al 2013).The WET value for reconstruction voxel r k , g r k ( ), is the weighted sum, over all pencil beams and relevant detector locations, of the WET associated with each location r d ,W z d ( ): )is the joint likelihood of equation (9), and W z d ( )is obtained from the depth z d of the detector pixel: with W 0 the WET of the PB and RSP det is the RSP of the detector.Equation (10) includes a summation over all relevant detector locations.In the simple case of figure 2(b), a single pixel is used; in the case of figure 2(c), where multiple structures of different WET are in the way of the beam, the three peaks should be identified and a single pixel reprojected for each.To identify peaks, the 2D signal is integrated over the x or y axis to generate an integral depth dose profile.A peak finding routine (Du et al 2006) is applied to find the z position(s) of candidate peaks, and the corresponding lateral (x d or y d ) central position of each peak is obtained from the pixel with the maximum intensity in the lateral profile in the vicinity of the z position.Then, the depth to reproject, z d , is estimated by taking the pixel in the IDD with the 80% distal intensity fall-off from the integral depth dose profile.The 3D position is then obtained by matching the peaks of both lateral views (XZ and YZ) according to the distance between their coordinates z d .Images were acquired at an energy of 200 MeV using a spot size of s s = = z z x i y i ( ) ( ) = 3 mm, and angular divergence of s s = = q q z z i i x y ( ) ( ) 3 mrad.Raw data is acquired with pencil beams spaced by 1 mm.This creates a large dataset which can be sub-sampled to study the impact of beam spacing on image quality.
Two phantoms were imaged: one phase of the extended cardiac phantom (XCAT) (Segars 2010), and an aluminum cube (5 cm of side) with a slanted edge at 2.5°placed inside a 10 cm water tank.The XCAT was acquired with a 30 × 30 cm 2 field of view (FOV), and is reconstructed using spacings of 2-6 mm in increments of 1 mm, to evaluate the impact of beam spacing.The slanted edge was acquired using a 15 × 15 cm 2 FOV and reconstructed with a 3 mm spacing to report resolution with a spacing that limits imaging dose.As the FOV of the XCAT matches the size of the detector, partial signal loss was observed in some PBs at the edge of the scintillator.All algorithms where however robust to partly losing data from the Bragg curve, as illustrated in figure 3.

Experimental datasets
Proton radiographs were acquired at Mayo Clinic Arizona using a 15 × 15 × 15 cm 3 plastic volumetric scintillator (EJ-260, Eljen Technologies, Sweetwater, TX) and a CMOS camera (PGE-23S3M-C, Teledyne FLIR, Wilsonville, OR) with a 25 mm lens (M118FM25, Tamron, Commack, NY), housed in a light-tight box made from black acrylic plastic.Projected 2D lateral views of the 3D energy deposition in the scintillator were acquired; the two lateral views (XZ and YZ planes) were obtained sequentially by rotating the couch by 90 • between acquisitions.The imaging field of view was set to 13 × 13 cm 2 .The following phantoms were scanned: Las Vegas (Tendler et al 2021) for contrast, 9 tissue-substitute cylindrical inserts (CIRS 062M electron density phantom, Sun Nuclear, Norfolk, VA) for WET accuracy, and a paediatric head phantom (Sun Nuclear ATOM Phantom, Model 704) for general image quality with a thicker object (>12 cm WET).The 9 tissue-substitute inserts were: lung (inhale), lung (exhale), adipose, plastic water, muscle, liver, bone 200 mg/cc, bone core 800 mg/cc, and bone 1500 mg/cc.Data was acquired using clinical settings at 135.6 MeV for all phantoms but the paediatric head, which was acquired at 189 MeV.Beam widths s z l ( ) at isocenter are respectively 3.1 and 2.5 mm for 135 and 189 MeV.Data acquisition was done using pencil beam spacings of 2, 3, 4 and 5 mm.A sampling of 6 mm is also generated by sub-sampling the 3 mm acquisition.The experimental setup is shown in figure 2(d); the detector was set on the treatment couch, and an additional couch was inserted between the nozzle and detector to hold phantoms.The distance from the nozzle to the top of the phantom was set to ≈20 cm, and the phantom to scintillator distance was 11 cm.Finally, the imaging dose as well as dose considerations are reported in appendix C.
Raw data from the CCD camera is corrected for multiple optical artifacts according to Robertson et al (2013).Briefly, images are corrected from perspective and refraction by projecting pixel intensities to the front of the detector according to geometrical and refraction principles.In addition, as the proposed image reconstruction framework assumes parallel beams, the small divergence angle of each beam is corrected by rotating the image by the negative of the divergence angle, assuming the center of rotation is the central position of the beam at = z 0

D
. Finally, all WET values are corrected by subtracting the WET of objects in the beam's path (5 mm for the entrance window of scintillator, 8.3 mm for the couch).

Comparison with other methods
Four reconstruction methods are compared in this work, and are referred to as distal, lateral 1D, lateral 2D, and single event.The distal method is an integrated-mode imaging method that uses the beam-eye view (XZ) with the setup of figure 2 for image reconstruction.The signal is integrated along the beam's axis, and does not provide a direct measurement of the WET.However, the intensity can be mapped back to the WET by producing an energy-specific calibration curve that relates scintillator intensity to WET; the approach presented in Darne et al (2022) is used.To limit multiple scattering induced blur, only the pixels whose intensity is  to 70% of the maximum intensity for each pencil beam image are reprojected.
The lateral 2D method is the one proposed in this study, which makes use of the two lateral views.To illustrate the benefits of the proposed lateral approach, a simpler lateral 1D method based on the work of Rescigno et al (2015) is also implemented.Briefly, the signal from one of the lateral views (either XZ or YZ) is summed over the lateral coordinate to produce a percent depth light curve, which reproduces a range telescope measurement.For reprojection, the approach of Rescigno et al (2015) is used, which is a specific case of the proposed algorithm; the reprojection equation is the same as equation ( 10), but the likelihood is directly  r r k i ( | ), as the detection event is assumed to have a probability of 1.The last method is a single event tracking method based on a front-tracker-binning method (Volz et al 2020), and is only used on the Monte Carlo data to compare the three integrated mode methods with a reference method in proton imaging.Plane tracking detectors are assumed ideal and only score primary generated particles (in Geant4, CreatorProcess 0) as an idealised particle imaging filter.

Data analysis
The XCAT is characterised in terms of WET accuracy.Absolute error maps D WET , calculated between the ground truth XCAT WET and the reconstructed WET as well as the mean absolute error (MAE) calculated over the entire image are reported.The slanted edge phantom is used to estimate the resolution by following the methodology of Fujita et al (1992); the resolution is defined as the spatial frequency at which the modulation transfer function (MTF) drops to 10%.The Las Vegas phantom is composed of 28 holes with variable diameter and depth.The contrast is reported for all visible holes of the phantom, in a manually defined region of interest in the inner 50% of each hole; missing data in figure 5(c) should be interpreted as not visible holes in the resulting radiographs.
The contrast C is defined as ˜, where W ref is a reference WET value in the close vicinity of the hole, and W W , min ˜are respectively the lowest and mean WET values in the ROI.W ref was estimated manually for each hole by considering profile lines around the center of each hole and averaging the intensities on each side of the hole.
For WET accuracy, the ground truth WET of each cylindrical insert was measured using a scanning water tank and a multi-layer ionisation chamber.Estimated values from radiographs are taken as the average in a region of interest representing the inner 70% of each insert to limit edge effects.The mean absolute and relative errors over all inserts are reported, along with the standard errors on the mean (SEM).For N inserts, each with standard deviation s n , the SEM is defined as 2 .Experiments described in sections 2.4-2.6 are summarised in table 1.

Results
Monte Carlo results are first presented in figures 3 and 4, while figures 5-7 illustrate experimental results.Finally, the impact of beam spacing is reported in figure 8, for both Monte Carlo and experimental datasets.
Figure 3 compares the general image quality and WET accuracy of the proposed reconstruction method with other reference methods introduced in section 2.5 for the XCAT phantom, using a beam spacing of 3 mm to limit imaging dose and time.
The Las Vegas phantom is shown in figure 5, where general image quality and a quantitative analysis of contrast are reported for the lateral 1D and lateral 2D methods.
Averaged over all hole diameters and hole depths, shown in figure 5(c), the lateral 2D method provides an average relative increase in contrast by 36% compared to the lateral 1D method.Figure 6 illustrates the WET accuracy of each tissue substitute insert estimated with the lateral 2D approach proposed in this work.Averaged over the 9 inserts, the mean relative WET error is 1.2 ± 0.3%, or 0.4 ± 0.1 mm in absolute.While not shown, the WET accuracy of the lateral 1D method was found to be similar, at 0.5 ± 0.1 mm.
WET maps from the proposed reconstruction method for the paediatric head phantom are shown in figure 7 for multiple beam spacings.This qualitatively illustrates the image quality that can be achieved with the current setup on a realistic geometry.Figure 8 shows, for the proposed reconstruction approach, the impact of beam spacing on reconstruction for (a) the simulated XCAT phantom and (b) the experimentally acquired Las Vegas phantom.

Discussions and conclusions
The results obtained in this work highlight the potential of the proposed reconstruction framework with 2D lateral projections to maximise image quality for integrated mode proton radiographs.While the distal view approach has been used for scintillation-based integrated mode detectors (Darne et al 2022), Monte Carlo simulations demonstrate that the distal view provides the overall lowest proton radiograph image quality compared to other views.Quantitative accuracy for the XCAT phantom as well as image resolution estimated from the slanted edge phantom are the lowest compared to all other methods evaluated.The blurriness visible in the distal WET maps of figures 3 and 4 is mainly attributed to MCS, as the integrated signal reprojected towards the imaging plane in the distal view is largely contaminated by MCS.The impact of MCS is limited for lateral views-while uncertainty due to MCS is considered to generate  r r r , ), the reprojected data itself is not  largely contaminated by MCS.In figure 3, it is noticeable that integrated mode methods produce more biased images than single event methods.This can be attributed to the thick, heterogeneous geometry (XCAT phantom) that result in highly complex dose distributions.While range mixing can be deconvolved in relatively simple heterogenous geometries (for instance, geometries that would produce dose distributions such as figure 2(c)), it is more difficult to manage with thick scatterers that produce a broad dose distribution.Indeed, the peak finding routine can fail in the situation where individual peaks associated with different tissues or interfaces cannot be resolved from one another.Results highlight the clear improvements obtained with the proposed method (lateral 2D) compared to a conventional treatment for lateral data (lateral 1D).The 1D method is equivalent to a range telescope that reprojects a single WET value per PB towards the imaging plane with weights  r r k i ( | ) (equation ( 3)), per Rescigno et al (2015).The  proposed approach can be seen as a generalisation of Rescigno et al to a 2D detector.In terms of WET accuracy, while it was found that the experimental WET accuracy on 9 uniform tissue inserts is similar between the lateral 1D and 2D methods (0.5 versus 0.4 mm mean error), it was also reported that the mean absolute error in WET is reduced by 15% in the simulated XCAT phantom (3.2-2.7 mm) with the lateral 2D approach.As the XCAT accuracy is calculated over the entire image, this suggests that the lateral 2D approach is better at handling heterogeneous geometries.Image resolution, calculated on the simulated slanted edge, also improves relatively by 22%.Finally, based on experimental radiographs of the Las Vegas phantom, contrast improves on average by 36%.The improvement from lateral 1D to lateral 2D is mainly attributed to two factors.First, the difference in the reprojection point spread function-namely equation (9) in the lateral 2D case, and equation (3) in the lateral 1D case.The point spread function of equation (9) has the advantage of modelling scattering inside the detector and provides a more realistic model of transport from the object to the detector.In other words, the current approach accounts for the full 3D path of the beam, whereas Rescigno et alʼs work only considers the initial position of the beam.Second, the lateral 2D approach with both views also has the advantage of spatially localising the position of multiple Bragg peaks in a single dataset, acting similarly to a pixel detector (Krah et al 2018b), limiting the impact of range mixing effects.
Experimental image quality is found to be maximised with the proposed reconstruction framework.Excellent WET accuracy of 0.4 mm, calculated over 9 tissue-equivalent inserts, is obtained.This value is consistent with Krah et al (2018a), who have reported 0.2-0.5 mm WET error on tissue-equivalent inserts using a multi-layer ionisation chamber (Qube by De.tec.tor,Turin, Italy), and an improvement upon (Deffet et al 2020), which have reported 1.5mm WET accuracy with a deconvolution approach.Results also appear at least as accurate as the carbon-ion radiographs of Magallanes et al (2019), which reported WET relative errors below 1.3% for most of 6 tissue-equivalent slabs, although adipose and lung exceeded 3%.
An advantage of the proposed imaging framework is that high quality radiographs can be acquired and reconstructed rapidly.Oria et al (2023) have used a flat panel based integrated mode system combined with energy  scanning.Using 10 energies, which results in an imaging time of approximately 40 min for a 27 × 27 cm 2 field of view, a mean absolute WET error of approximately 5 mm is obtained in simulated datasets.In the present study, assuming an acquisition of 3 ms per pencil beam (Langner et al 2017), which could be obtained with the low dose measurements reported in table 2, the full 30 × 30 cm 2 field of view of the XCAT phantom with a 3 mm spacing can be acquired in 30 s, with a mean absolute WET error of 2.7 mm.The head phantom shown in figure 7 could be obtained in 5.5 s with a beam spacing of 3 mm.This opens applications for adaptive proton therapy in terms of patient positioning, in vivo range verification, and adaptive re-planning.Rapid radiographic imaging of moving tumours can also be considered-for instance, the 10 × 10 cm 2 square ROI identified in figure 8 with the 6 mm spacing could be acquired in 0.9 s.However, further investigations are required to explore usability for tracking purposes, as the coarse 6 mm sampling can result in distorted image features, as highlighted in figure 8(b).
This study comes with some limitations.First, the proposed approach does not simulate nor correct optical blurring that arises with scintillation-based detectors (Tendler et al 2021).Imaging quality obtained in silico may therefore be overestimated compared to current experimental results.Incorporating optical blurring into the proposed formalism could however lead to improved experimental image quality.Furthermore, fast imaging is conditional on (1) the use of a high-speed camera (frame rate > 330 Hz (Langner et al 2017), which is commercially available), (2) a triggering system (either hardware based (Alsanea et al 2019) or software based) to synchronise image acquisition with pencil beams, and (3) a fast reconstruction framework.Due to limitations in frame rate of the current camera, experimental datasets were acquired at relatively high dose per PB to increase acquisition time per PB, which is not clinically practical.However, we have tested, with single pencil beams acquired at the doses presented in table 2, that the estimation of z d is robust to noisy images, and we do not anticipate any problems with lower dose images (see section C).The current reconstruction time is in the tens of seconds for a 3 mm spacing, and should be improved via GPU acceleration.Future work therefore include optimising both hardware and software components to explore adaptive radiotherapy applications.Furthermore, the current detector setup (figure 2(d)) constrains the gantry angle to be at either 0 or 90 degrees (or 180 and 270), which may limit applications in treatment adaptation.
In this work, it was found that the spatial resolution estimated on the slanted edge phantom varies minimally with beam spacing, for instance ranging from 0.35 (lateral 2D, 1 mm spacing) to 0.33 lp mm −1 (lateral 2D, 3 mm spacing, shown in figure 4).The main improvements in resolution are obtained by modifying the image reconstruction strategy.The relationship between spatial resolution and beam spacing requires further investigation-for small objects, larger beam spacing leads to distorted images (figure 8), which can have unexpected impact on phantoms such as the slanted edge.We recommend that further investigations are performed with adequate line pair modules.
This work introduces a novel reconstruction framework for integrated-mode proton radiographs using 2D lateral views of individual pencil beams, which provides improved image quality compared to using a range telescope.High quality proton radiographs can be obtained with clinical beam settings using a fast, low-cost scintillation-based system with two cameras capturing both lateral views.

Figure 1 .
Figure 1.XZ projection of the PB scanning setup.The PB propagates z O cm in air before traversing a WET W z d( ) within the object.An additional distance (calculable using z D ) is traveled in air before energy is deposited at depth z d in the detector.Details on how the hull of the object is considered are given in appendix B.

Figure 2 .
Figure 2. (a) Scintillation detector setup illustrated without an object resulting in a 2D detection system.The camera captures a XZ projection (at camera angle f = 0) of the 3D energy deposition profile, resulting in the detection element shown in the bottom right corner of the scintillator.(b) Examples of projections with the top (YZ) and lateral (XZ) views, where the central pixel of the 2D Bragg curve representing the 80% distal fall-off is identified in green in each plot.(c) Additional example illustrating a case where the pencil beam crosses structures (not shown in the figure) with different WET and result in multiple peaks.(d) photograph of a proton radiography experimental setup; see section 2.4.2 for details.
is obtained by integrating over the detector dimensions.Using the example of the 2D XZ projection shown in figure 2, one obtains the YZ views shown in figure 2 are combined to obtain = x y z r

2 D
2.4.Data production 2.4.1.Monte Carlo simulations Monte Carlo simulations were performed with the Geant4 software toolkit (version 10.6) (Agostinelli 2003), using the same geometry as presented in figure 1. Detailed simulation parameters following AAPM TG-268 guidelines (Sechopoulos et al 2018) are the same as in table 1 of Tendler et al (2021), with the following differences: the source geometry consisted of two 2D Gaussian functions, the first Gaussian, containing 75% of the beam intensity, and having a spread corresponding to the spot size and the second Gaussian having a spread twice the spot size with the remaining 25% intensity.Angular spread was sampled from a Gaussian distribution characterized by a defined beam divergence.There were 4 10 6 histories per pencil beam position, which were centered on evenly distributed points across the field of view.A distance of = z cm 40 D was used and the object was centered at / z from the beam.A 30 × 30 × 30 cm 3 detector object with the density and composition of the EJ-260 plastic scintillator (Darne et al 2022) was used to simulate the volumetric scintillation detector and score energy, which was transformed to quenched light emission following Birks' law (Tendler et al 2021).

Figure 3 .
Figure 3. WET maps (top) and absolute error D WET (bottom) in mm of the simulated XCAT phantom using a 3 mm beam spacing, for the four reconstruction methods (left to right).Lateral 2D is the proposed approach of this work.

Figure 4 .
Figure4.Top: water equivalent thickness maps (mm) of the simulated slanted edge phantom using a 3 mm beam spacing, for the four reconstruction methods (left to right).Bottom left: resulting edge spread functions (colour) and associated fit (dashed lines).Bottom right: modulation transfer functions for each method.

Figure 5 .
Figure 5. (a) Experimental WET maps of the Las Vegas phantom for a 3 mm beam spacing using the lateral 1D and 2D methods.The dynamic range is from 53 to 59 mm.The observed intensity dip in the lateral 1D figure (top right of profile line p 2 ) is due to a noisy signal from a single PB.(b) Line profiles corresponding to the p 1 and p 2 lines identified in sub-figure (a). (c) Corresponding calculated contrast for each visible inserts of the Las Vegas phantom for each reconstruction approach.Data is shown as a function of hole diameters for various hole depths d.

Figure 6 .
Figure 6.Experimental WET accuracy for 9 tissue substitute inserts, estimated with the proposed lateral 2D approach.Relative differences on the WET are shown along error bars representing the standard deviation inside the selected region of interest.

Figure 7 .
Figure 7. (Left) photograph of the paediatric head phantom, (right) experimental WET maps (mm) of the phantom reconstructed using the lateral 2D method for beam spacings of 2, 3 and 4 mm.

Figure 8 .
Figure 8. Water equivalent thickness maps reported as a function of beam spacing for (a) the simulated XCAT phantom and (b) the experimental Las Vegas phantom, reconstructed with the lateral 2D method.For scale, the red square in (a) represents a 10 × 10 cm 2 region of interest.
are first summarized in appendix A. Expressions for  r r

Table 1 .
List of scanned/simulated phantoms with the reconstruction methods and performance metrics used.