Proton radiography and proton computed tomography based on time-resolved dose measurements

The principle of measuring WEPL based on time resolved dosimetry has been described previously (Lu 2008, Gottschalk et al 2011). Here we will briefly review the method for specific applications to proton radiography.

In most of the proton therapy systems with passive scattering, like the facility at Massachusetts General Hospital, modulator wheels are used to produce the range modulation for the construction of spread out Bragg peak (SOBP) fields (Lu and Kooy 2006). A modulator wheel may contain one or a few circular modulation tracks each consisting of a series of stepwise segments of absorber material (typically two different materials per step, e.g. Lexan or carbon combined with lead) with increasing water equivalent thickness. When a mono-energetic proton beam passes a particular step, its range is reduced, or pulled back, by the water-equivalent thickness of the segment.

For our gantry-mounted nozzle system, the modulation wheel spins at a constant speed of 600 rpm. The beam passes through the wheel one segment at a time along the track, delivering a series of Bragg peaks spread out in depth. The dose rate, as function of time, produced in a medium by such a beam is thus periodic with the period of the wheel rotation, i.e. 100 ms for our system (Lu 2006). The time dependence of such dose rate functions (DRFs) has been measured in previous works inside a water tank with a small ionization chamber (Lu 2008) and with semiconductor diodes (Gottschalk et al 2011). The exploitation of the characteristic patterns of the DRFs, at different depths within the dose plateau of a SOBP, was originally proposed for proton radiography applications in Lu (2008). Indeed, by measuring the dose rate at any point distal to a phantom, the residual proton range can be potentially determined with millimeter accuracy.

In this work we used a prototype, two-dimensional diode-array detector (Sun Nuclear Corporation—Melbourne, FL, USA) consisting of 249 semiconductor diodes with 2 ms sampling-time displaced on a 12 cm side matrix. The diodes are arranged in an octagonal shape layout (see figures 2 and 9(c)) with 7.07 mm diagonal pitch. The linear pitch between diodes on the same row and column is 10 mm and between each row, diodes are shifted by half linear pitch (5 mm). Therefore 12 or 13 diodes, depending on the row and column, are aligned along the longest dimension of the octagon (12 cm) and 5 diodes are aligned along the shortest side of the octagon (4 cm).

Zoom In Zoom Out Reset image size

In order to determine the residual range of protons traversing phantoms with arbitrary thicknesses we first need to acquire a database of DRFs measured at known WEPL. To acquire such a database, we used a water-phantom (DigiPhant PT, IBA, Belgium) to measure DRFs in steps of 1 mm for WEPL from 5.5 cm to 20 cm. To accurately determine the WEPL of protons reaching the diodes array in the detector, the water equivalent thickness of the water-phantom walls was determined through separate depth–dose range-shift measurements. DRFs measurements were performed with a SOBP with 20 cm range and 20 cm modulation width. With such a field and a beam current of 10 nA at the exit of the cyclotron, the dose rate, measured at the center of the SOBP, is 0.72 cGy s⁻¹.

Figure 3 shows measured DRFs at different WEPL in the water-phantom for one of the diodes in the central area of the detector. The entire technique relies on the fact that for each modulator track the DRF has a unique pattern at every point of measurement, and therefore encodes the proton WEPL from the front surface of the phantom to the position of the diode detectors (Lu 2008).

Zoom In Zoom Out Reset image size

Once the DRF-database is acquired, a single DRF-measurement at a point distal to the phantom is sufficient to determine the WEPL through the unknown phantom thickness. Indeed, such DRF will have a characteristics pattern encoding the protons WEPL at the point of measurement. To determine its WEPL, the DRF will need to be compared to the DRFs in the database. Since all the DRFs in the database have been acquired at known WEPL, the one that produces the best pattern-match will provide the WEPL through the phantom. Accordingly, the residual proton range behind the phantom is given by the difference between the nominal range of the beam and the obtained WEPL.

Depending on the characteristic sampling-time of the detection system, two methods have been developed to compare DRFs: the pattern matching technique (Lu 2008) and the rms-width technique (Gottschalk et al 2011). In this work we used the pattern matching technique that is based on the minimization of the least-square difference between the DRF measured for an unknown proton WEPL and all the DRFs of the database. The pattern matching is done individually for each of the 249 diodes thus providing a 2D radiography of the phantom expressed directly in WEPL.

Figure 4(a) shows the mean WEPL of the 249 diodes, determined through the pattern matching technique, for DRFs measured in the water-phantom used for the acquisition of the database, from WEPL of 5.5 cm to 21 cm. Data are plotted versus the nominal WEPL of the points of measurement that was provided by the DigiPhant control system. In the ideal case of perfect matching between newly measured DRFs and DRFs in the database a straight-line correlation is expected. Up to a depth of 194 mm the mean WEPL does match the nominal WEPL, with the exception of only three points for which we can reasonably assume that fluctuations have occurred in the detector electronics or in the proton beam. Distal to 194 mm, the correlation between measured and nominal WEPL is lost. Indeed, diode signals get significantly attenuated because protons range out and the pattern-matching algorithm fails in determining the WEPL. Since the SOBP field has a nominal range of 20 cm, we can then infer that the detector intrinsic build-up is 6 mm. Such a value was used to correct the water equivalent depth for all the DRFs in the database. Figure 4(b) shows the standard deviation of the mean WEPL calculated among all 249 diodes. For each point, fluctuations around the mean WEPL are on the order of a millimeter. Indeed, provided that the field has homogenous proton range in the plane transverse to the beam, the fluctuations in the DRFs measured by diodes at the same WEPL but at different radial position inside the detector are minimal. Therefore we can assume that in a pure water phantom, millimeter accuracy in the determination of the WEPL can be achieved at any radial detector location. As general trend, the standard deviation decreases with the increasing of the WEPL until 194 mm when protons range out and the stdv rises rapidly and goes out of the plotted scale. The initial decrease of stdv is mainly due to the fact that the DRFs pattern changes less significantly for shallow depth and therefore the matching algorithm has higher possibilities of returning a less precise WEPL value. The same pattern in the variation of stdv was observed in a previous study based on ion chamber measurement (Lu 2008).

Zoom In Zoom Out Reset image size

The general experimental set-up for both proton radiography and tomography is shown in figure 2(a). The detector is placed distal to the phantom and aligned to the center of gantry nozzle with the laser alignment system regularly used for patient positioning. For each experiment we kept the air-gap between the phantom and the detector to ∼3–5 cm to minimize proton scattering. A square brass aperture was used to produce a 12 × 12 cm square proton field in order to match the diodes active area and to reduce the irradiation to the detector electronics. A fully modulated SOBP with nominal range of 20 cm was used for all irradiations. Depending on the most suitable set-up arrangement, according to specific phantoms, irradiations were performed with the gantry nozzle at 0° (vertical) or 270° (horizontal).

Proton radiographies were acquired for a series of phantoms with varying complexity: (i) Lucite (density 1.18 g cm⁻³) wedge phantom with square base (side 20 cm) and slope 40° (see figure 6(a)); (ii) Lucite step-like phantom with square base (side 12 cm) and 36 square steps (side 2 cm) with an increasing thickness of 2 mm per step (see figure 6(b)); (iii) hollow plastic ball (diameter 7 cm) filled with water (see figure 8); (iv) human skull (see figure 9).

For the database the shallowest DRF was acquired for a WEPL of 5.5 cm. Thus, solid water shifters were put in front of the detector in order to reach, for each phantom configuration, the minimum required thickness of 5.5 cm.

In order to increase the detector spatial resolution, currently limited by the quite large linear pitch between diodes (1 cm), we placed the detector on a remote controlled motion platform. We then moved the platform step by step, along a 5 × 5 (step-like phantom and plastic ball) and 10 × 10 (human skull) square mesh-grid. For each detector position we acquired one radiography and then merged together the total number of radiographies in order to get images with 'virtual' linear pitch between diodes of 2 mm for the 5 × 5 mesh-grid and 1 mm for the 10 × 10 mesh-grid, respectively.

The minimum time to produce one radiography is given by the revolution period of the modulator wheel (100 ms) that determines the periodicity of the DRFs. To investigate the potential of using the proton radiography technique for real-time imaging of a moving phantom we placed a Lucite cube (side 4 cm) on the remote controlled motion table and acquired serial radiographies, every 100 ms, while the cube was moving above the static detector. Two motion sets were investigated: (A) 2D sinusoidal motion along both x and y-axis with cycle amplitude of 30 mm and frequency of 10 cycles min⁻¹; (B) 1D sinusoidal motion along the y-axis with cycle amplitude of 30 mm and frequency of 30 cycles min⁻¹. We merged all the radiographic images in two movies (available from stacks.iop.org/PMB/58/8215/mmedia), corresponding to motion options A and B, respectively. Each movie frame is constituted by a single radiographic image, acquired in 100 ms. The detection system, acquires data continuously, with an integrated sampling time of 2 ms. Therefore, to create individual movie frames we simply had to re-sample a single data file, recorded over 10 s, into time frames of 100 ms each, according to the internal clock of the detector.

We performed a pCT experiment with the cylindrical phantom sketched in figure 5. The Lucite cylindrical phantom (12 cm diameter, 12 cm height) contained ten rods inserts arranged differently between the upper and lower phantom sections. The inserts were made of bone equivalent material (Gammex 450 cortical bone, density 1.82 g cm⁻³), lung equivalent (Gammex LN-300, density 0.30 g cm⁻³), air and Lucite (same as phantom).

Zoom In Zoom Out Reset image size

The phantom was placed on a remote controlled rotation table as shown in figure 2 and 100 radiographic projections were acquired at angles θ = 0, 1.8°, ...,178.2°.

To increase the spatial resolution, images were reconstructed using interleaved projection measurements of two full rows of the detector, providing 25 radial samples at r = −6, −5.5, ... +6 cm. We denote as f the 2D reconstructed image and as p(θ,r) the interleaved projections. A system matrix A was constructed that describes the geometry of the projection acquisition such that p = Af. The system was modeled as a fan-beam geometry considering the effective source axis distance of the double scattering proton therapy system and the distance from isocenter to the detector plane. The projections itself were modeled as ideal straight lines; we did not include proton MCS in the image reconstruction model.

The reconstruction of the pCT images was performed by iteratively minimizing the following unconstrained objective function:

$$\begin{equation*} {\rm min}\{ {|Af - p|2 + \lambda \;{\rm TV}(f)} \} \end{equation*}$$

in which the first term represents the consistency of the reconstructed image with the measured projections in a least square sense, and the second term is the total variation of the reconstructed image which serves as regularization function.

The reason for using regularization in the image reconstruction arises from the limited number of projections that can be acquired with the current detector, which renders the image reconstruction an ill-posed problem. The parameter λ determined the trade-off between finding a smooth solution and the consistency of the solution with the measurements. As we only performed a single pCT scan at this time, the value of the parameter was determined empirically.

The optimization was implemented employing the accelerated proximal gradient solver NESTA (Becker et al 2011, Nesterov 2004). We previously tested this implementation for the reconstruction of images from incomplete x-ray projections (Verburg and Seco 2012).

Figure 6 shows the proton radiographies of the Lucite wedge phantom and step-like phantom. Radiographies are directly expressed in WEPL and plotted as 3D histograms. Each bar of the histogram represents the WEPL determined by pattern matching of the DRF provided by a single diode; therefore each image contains 249 bars. In figure 6(a), measured WEPL values were fitted with a plane. The slope of the fitting plane (40.3°) agrees well with the physical slope of the wedge phantom (40°). In figure 6(b), the proton radiography of the step-like phantom shows clearly the step-like pattern although individual steps are not distinguishable. This is practically due to the relative low increase in Lucite thickness for each step (2 mm) and, physically, to the effect of protons MCS on the WEPL resolution.

Zoom In Zoom Out Reset image size

In general, the WEPL resolutions of proton radiographies obtained with our time-resolved technique depends, among others, on the data acquisition time. Because of unavoidable fluctuations between DRFs from different wheel modulation cycles, it can be beneficial, to average the DRFs for each diode, over several repetition cycles. This allows producing DRFs that are less affected by fluctuations of the detection electronic or temporary fluctuations in the proton beam current. All the radiographies of static phantoms presented in this section were obtained through the average of ten DRFs from ten successive modulator cycles. The overall integral acquisition time was 1 s.

Figure 7 shows the proton radiographies of the step-like phantom acquired by shifting, step by step, the detector along the 5 × 5 mesh-grid in order to increase the spatial resolution of the image. The resulting 25 images have been merged in the single radiography shown in figure 7(a). In such radiography the number of WEPL points is then 6225 (25 × 249) and the linear pitch between diodes is 2 mm.

Zoom In Zoom Out Reset image size

The MCS of protons traversing the imaged object is the main cause for degradation of the spatial resolution in proton radiography. The many small-angle deflections that protons undergo along their trajectories lead to so-called proton range mixing (PRM) (Bentefour et al 2012). The range mixing increases with the WET of the imaged object and gets intensified by its material and density inhomogeneity. Practically, when imaging a rather complex object like the step-like phantom, it is unavoidable to obtain DRFs with 'distorted' patterns. The pattern distortion arises from the contribution of protons that have travelled different water equivalent paths but that, because of their scattering through different steps of the phantom, end up being detected, altogether, by a single diode. These DRFs will then not be present in the database and their own pattern will depend on the degree of PRM. For such DRFs the pattern-matching algorithm could potentially provide inaccurate WEPL.

The degree of pattern dissimilarity, between DRFs measured behind the phantom and the reference DRFs in the database, can be used to assess the amount of PRM for each point of the radiography. As an example, for all the WEPL points of figure 7(a), we calculated a so-called range-mixing factor defined as the least-square difference between the DRF measured behind the phantom and the DRF in the database provided as 'best-match' by the pattern-matching algorithm. In figure 7(b) we then re-plotted the WEPL points of figure 7(a), but only for those diodes whose range-mixing factor would not exceed the 50% of the maximum value among all the points. This resulted in about ∼30% reduction of the total number of points of figure 7(a). The new radiography presented in figure 7(b) is much sharper than the homologous image presented in figure 7(a), with clearer thickness contrasts and easier identification of the step-like pattern. As expected, the highest range mixing is produced in the phantom area with the highest thickness gradient corresponding to the regions of the 'vertical cliffs' between two nearby ramps of steps.

Figures 8(a)–(b) shows the proton radiographies of the plastic ball filled with water plotted as a 3D histogram and as planar image, respectively. As for the case of step-like phantom of figure 7, 25 images were acquired over a square mesh-grid and merged in a single radiography. Moreover, we performed a bilinear interpolation of the measured WEPL points on a 1 × 1 mm interpolation-grid. The general shape of the ball is well reproduced. In figure 8(b), regularly spaced cold-spots are due to the fact that, out of the 25 merged image-files, one was partially damaged due an electronic malfunction of the detection system.

Zoom In Zoom Out Reset image size

Figure 9 shows a comparison between proton radiographic images ((a), (b)) and a conventional x-rays radiography (c) of the human skull. The x-rays radiography was acquired with the x-ray patient positioning system used in clinical routine. The image shows the skull placed over the sensitive detection area of the Sun Nuclear detector. The octagonal layout of diode detectors is clearly visible. Proton radiographies shown in figure 9 are obtained merging 100 individual images acquired shifting, step by step, the detector over a 10 × 10 square mesh-grid. The planar image of figure 9(b) was further interpolated on an interpolation grid of 0.1 × 0.1 mm. In figure 9(b), regularly spaced cold spots are due to one partially damaged file out of the 100 constituent images. Moreover, in both figures 9(a) and (b), hot and cold square areas in the regions around the positions (−2,−2),(−3,−1) and (−4,−3) are due to three damaged diodes in the detector.

Zoom In Zoom Out Reset image size

The spatial resolution of conventional x-ray images is superior to proton radiographic images presented in figures 9(a)–(b). This is not only due to proton MCS, but in this specific case, also and foremost, to the different intrinsic spatial resolutions of the two different detectors. Sub-millimeter spatial resolution is easily achievable with commercial x-rays panel while the spatial resolution of the detector is currently limited mainly by the linear pitch between diodes (1 cm).

In the supplementary data we present two movies (available from stacks.iop.org/PMB/58/8215/mmedia) obtained merging the serial radiographies acquired in the experiment performed with the Lucite cube moving with 1D and 2D sinusoidal motion, respectively. Each movie frame is constituted by a single radiographic image, acquired for one rotation cycle of the RMW (100 ms) and interpolated on a 5 × 5 mm grid.

Figure 10(a) shows an example of a single movie frame for the case of the Lucite cube moving with the 1D sinusoidal motion described in section 2.3. For this single frame-radiography the dose delivered to the phantom was integrated for only one cycle of the RMW and it is therefore ten times lower compared to the images of static phantoms presented in the previous section. Overall, for both movies (available from stacks.iop.org/PMB/58/8215/mmedia), the shape of the cube is clearly distinguishable and its trajectory distinctly recognizable.

Zoom In Zoom Out Reset image size

Figure 10(b) shows the motion amplitude diagrams for both motion sets along the x- and y-axis. Dashed lines represent the position of the motion table provided by the table control system. Solid lines represent the position of the center of the cube and are obtained applying to each movie frame a contour algorithm used to determine the coordinates of the center of the cube. Table and cube positions are within millimeter agreement for both motion sets. Little mismatches are mainly due to the non-optimal image spatial resolution mainly caused by motion. Indeed, as an example, for motion set B, the maximum cube speed is 15 mm s⁻¹ that results in a 1.5 mm translation during the 100 ms acquisition time.

Figure 11 shows the reconstructed proton tomographic images for the two sections (up and down) of the Lucite phantom shown in figure 5. The images are expressed in proton RSP with respect to water. Lung, air and bone rod inserts are clearly distinguishable. Lucite rod inserts are indistinguishable from the cylindrical phantom that is made of the same material.

Zoom In Zoom Out Reset image size

To assess the accuracy of RSP obtained from the reconstructed pCT images, a comparison with RSP values obtained from standard range-shift measurements is presented in table 1. The reconstructed RSP for Lucite is close to the actual value, while the RSP values of the inserts are somewhat biased toward to RSP of Lucite. The reason for this bias is twofold. First, the spatial resolution of the prototype detector is limited. Second, the present reconstruction model does not consider proton MCS. A detector pixel behind one of the inserts detects not only protons that passed through the inserts, but also protons which passed only through Lucite.

The spatial resolution was estimated as ∼7 mm, defined as the distance between the points at 90% of the RSP for Lucite (1.01) and 110% of the RSP for air (0.21). The image noise defined as the standard deviation of the pixel values in a 3 cm diameter circular region located at the center of the images was ∼0.1%.

In principle, radiographic images obtained by time resolved dosimetry could also be acquired with commercially available 2D ion-chamber/diode-array detectors. Nevertheless, the dose sampling time of such commercial detectors is usually too high (∼ 20 ms) to acquire accurate DRFs. The sampling time of the detector (2 ms) resulted to be good enough to produce fairly accurate proton radiographic images. Nevertheless higher sampling rates are strongly desirable. Indeed the sampling rate defines the accuracy in the determination of the pattern of DRFs, which are the fundamental basics of the pattern matching technique used to infer the WEPL of the measured points. The other fundamental parameter that affects the image quality is the spatial resolution of the detector that is currently limited by the quite large linear pitch (1 cm) between diodes. Nevertheless higher spatial resolution and higher temporal sampling rate can result in larger noise since fewer protons are detected each time. Therefore, a careful analysis of the dose–noise relationship is required.

All the techniques employed to increase the resolution, namely, merging several radiographies and interpolating the resulting image, can only partially bridge the gap with current conventional x-ray imaging technology. Nevertheless, proton radiographies, being directly expressed in WEPL, are quantitative images that can potentially bring valuable information on the residual range of protons traversing the imaged object.

Protons' MCS in the imaged object, leading to PRM, is the main physical cause for degradation of the spatial resolution in proton radiography. We developed a basic method to assess the amount of PRM based on the degree of pattern dissimilarity, between DRFs measured behind the phantom and the reference DRFs in the database.

More advanced techniques for quantitative evaluation of range mixing are currently under investigation. As an example, one can imagine that any DRF can always be constituted by a linear combination of two or more DRF of the database and that the weighting factors of the linear combination can be used to determine the most likelihood WEPL of protons through the imaged object (Bentefour et al 2012). We could expect that the application of such techniques would potentially increase the quality of radiographies. The investigation of such techniques is beyond the scope of this work and will be presented in a separate study.

With the proton beam parameters employed in this work, the dose rate measured at the center of the SOBP is 0.72 cGy s⁻¹. The radiographies of static phantoms presented in the results section were obtained by averaging ten DRFs from ten successive modulator cycles to minimize beam and detection fluctuations. The overall integral acquisition time for each radiography, was then 1 s and the phantom received a dose of ∼0.7 cGy.

We investigated the influence of the acquisition integral time over the WEPL resolution in the case of the step-like phantom. For acquisition time down to 500 ms, corresponding to the averaging of five successive DRFs, radiographies of comparable quality of the one presented in figure 6(b) are easily achievable. For lower integration time a slight deterioration of the image quality is observed. Moreover, the proton beam current also plays a role on the image resolution. A too low beam current can result in low signals from the diodes and therefore distorted pattern of the DRFs. An investigation of the effect of proton beam current on image quality revealed that, with a proton beam current reduced to 3 nA (instead of 10 nA used conventionally) and an integrated acquisition time of 1 s, the resulting radiographies were of comparable quality of the one presented in figure 6(b). Therefore, we can reliably assume that, by reducing both the data integration time and the proton beam current, same quality images can still be produced delivering a dose to the phantom at least 2–3 times lower than the dose used to produce the radiographies presented in the results section. Nevertheless, further investigation for a quantitative assessment of the influence of the integral dose on the image quality is required. Moreover the dose can eventually be further reduced from the currently minimum value of ∼0.07 cGy/radiography (considering an acquisition time of 100 ms) employing more sensitive detectors. This would potentially bring an advantage for the application of the technique to pCT where some dose optimization is probably necessary due to the high number of projections required to obtain high quality reconstructed images.

While proton radiographic images could improve targeting accuracy for soft tissue tumors as shown in figure 1, the WEPL information in these images also offers a unique approach to beam range verification, for example, for the brain fields of pediatric medulloblastoma patients. For such patients, the whole brain is treated by two nearly opposing lateral fields. Ideally, each field should 'range out' right before the skin surface on the opposite side of the head, so that the hair follicles will receive only half of the target dose. However, to take into account range uncertainty, typical safety margins of 3.5%+1 mm have to be added to the prescribed beam range, potentially causing a significant increase in the skin dose. It has been recently shown that permanent alopecia, affecting some of these patients is correlated with high skin dose (Min et al 2013). Proton radiography based on WEPL measurements could be used to perform a pre-treatment 'range tuning' of such brain fields. A low dose 'scout-beam', having somewhat excess range in order to fully traverse the patient head, would be used with the patient already immobilized in treatment position, to produce a proton radiography of the cranium. By comparing the measured WEPL values to those computed by the treatment planning system, one could derive the exact beam range needed to fully cover the cranium target volume and yet spare the hair follicles, potentially reducing the risk of permanent alopecia. The quality of the images presented in figure 9 and the low dose required to produce one radiography is encouraging to envisage the safe clinical application of this technique for pre-treatment range verification and 'tuning' for medulloblastoma patients. Note that this 'pre-treatment range tuning' is required only once for each of the two lateral fields for the entire treatment course, since there should be very little interfractional variations of the WEPL values along the beam paths for this type of treatment if the patient is positioned correctly. An extra dose of ∼1–2 cGy for the 'tuning' would be negligible compared to the prescription dose for the treatment course (20–30 Gy).

The experiment performed with the moving cubic phantom, although realized in an ideal geometry, demonstrates the potential real-time imaging features of the technique. The method is rapid (∼100 ms) and simultaneous over the whole field and therefore it can potentially image mobile tumors without the problem of interplay effects inherently challenging proton radiography techniques based on pencil beams. Nevertheless, further experiments, with, for example, a realistic respiratory lung phantom, are required to assess the clinical applicability of the technique for real-time tumor tracking. One could imagine that lung patient would be eventually imaged for few respiratory cycles prior the treatment, in order to verify, and possibly adjust, the average breathing pattern depicted in the 4D CT used for treatment planning. In such a way the effects of inter-fraction variability of the breathing pattern would be potentially mitigated by adjusting the original plan according to measured variations at the time of treatment. Indeed, 4D inverse planning, incorporating pre-treatment patient-specific respiratory motion information into the treatment plan has the potential of quantifying patient-specific progressive geometrical variations and, more generally, to improve the treatment plan quality (Sonke and Belderbos 2010). Note that uncertainties in tumor localization are more important in particle therapy than in conventional x-ray therapy because charged particles beams are highly sensitive to geometrical and density variations in the beam path length. Therefore, repetitive proton radiographic imaging, expressed directly in protons' WEPL, have the potential of better targeting mobile tumors along with a general possibility of reducing the exposure of organs at risk, facilitating safe dose escalation, and improving local control as well as overall survival (Riboldi et al 2012).

One of the main potential applications of pCT is to provide more accurate proton RSP values for treatment planning. Conversion of x-ray CT numbers to RSP is subject to uncertainty which is estimated to be around ±2% (Paganetti 2012). Since even small discrepancies in the calculated RSP can result in significant changes in range, because they accumulate over the entire beam path (Jiang and Paganetti 2004), a method to directly obtain a 3D distribution of RSP in the patient could result in an increased accuracy of the delivered proton plan.

The images reconstructed from the pCT experiment presented in this work are generally of good quality. The limited spatial resolution of the prototype detector and the relatively low number of projections represent the main limitation to the accuracy of the RSP values obtained from the images. A full-scale detector with higher spatial and time resolution will improve accuracy.

Also, considering MCS in the image reconstruction is anticipated to yield better images. Previous authors have studied algorithms that include the proton path probability (e.g. Schulte et al 2008, Wang et al 2010, Rit et al 2013). Such models are usually applied to pCT systems where individual protons are tracked through multiple planes. Although this information is not available in the case of a single 2D-array detector, it is still possible to incorporate the scattering physics in the reconstruction model. For example, the model may be updated during each iteration, to consider for each pixel the probability of the various paths of detected protons.

The DRFs have a unique pattern, at each water equivalent depth, only within the dose plateau region of the SOBP (Lu 2008). Moreover the pattern of the DRFs in the DRF-database depends on the SOBP range and modulation width. For these reasons, in this work we used a long range, full-modulated SOBP. Only one database was then acquired for such SOBP and this same field was used to image all types of phantoms. For potential future clinical applications, one can imagine to use a full-modulated SOBP with the highest available range provided by the accelerator (in our system: ∼29 cm water equivalent). Therefore only a single database of DRFs will then need to be acquired. This same SOBP-field would then be used, for example, to perform the cranial radiographies for pre-treatment range tuning of any medulloblastoma patients

Because the periodicity of the DRFs is given by the rotation of the RMW, this proton radiography technique is only applicable to passively modulated systems. However, there is no need for producing a flat dose in the plateau region of the SOBP and therefore there would be no need for dedicated optimizations of the RMW geometry and beam current modulation. Potentially, a single RMW could be applied to image any patient site. Therefore one could imagine to develop a movable device for proton pencil beam systems in which a scatterer and a RMW would be sliding in the treatment head nozzle to acquire proton radiography before each irradiation fraction and be successively removed before the treatment delivery.

More generally, the use of proton radiography instead of currently used conventional x-ray imaging could have an advantage also in daily patient positioning. Indeed we can suppose that having a single beam would help suppress the misalignment due to possible discrepancies between the proton beam-eye-view and the current x-ray system beam-eye-view along with speeding up the patient set-up process (Depauw and Seco 2011).