Novel, full 3D scintillation dosimetry using a static plenoptic camera

Purpose: Patient-speciﬁc quality assurance (QA) of dynamic radiotherapy delivery would gain from being performed using a 3D dosimeter. However, 3D dosimeters, such as gels, have many disadvan-tages limiting to quality assurance, such as tedious read-out procedures and poor reproducibility. absolute dose difference between the reconstructed 3D dose and the expected dose calculated using the treatment planning software Pinnacle average 1.5% of maximum for both integrated IMRT and VMAT deliveries, and below 3% for each individual IMRT incidences. Dose agreement the 3D dose radiochromic ﬁlm acquisition same experimental on average within 2.1% and 1.2% of the maximum recorded dose for the IMRT and VMAT delivery, respectively.


INTRODUCTION
Recent advances in linac technology, on both hardware and software control, have enabled the development and widespread use of dynamic radiotherapy techniques such as intensity-modulated radiation therapy (IMRT), 1, 2 volumetricmodulated arc therapy (VMAT), 3,4 and tomotherapy 5 . The complexity of these delivery techniques calls for a personalized treatment quality assurance (QA), which is usually performed using 2D dosimetry. 6,7 However, as 2D measurements only represent a partial sampling of the planned 3D dose distribution, it is recognized by the medical physics community that only 3D dosimeters offer a truly comprehensive verification of dynamic dose delivery. 7,8 Moreover, it is our hypothesis that the validation of such dynamic treatment techniques would benefit from a real-time dosimeter, i.e., a detector which would record the incident dose distribution as a function of time.
Current dosimetric systems using polyacrylamide, 9,10 Fricke, 11,12 or radiochromic 13,14 gels offer desirable levels of precision and accuracy at a spatial resolution ranging from 0.5 to 5 mm. However, these systems have several disadvantages limiting to quality assurance. Indeed, most gel dosimeters involve laborious and time consuming preparation and processing manipulations, leading to variations between gel batches, as well as delayed dose information. Moreover, gels are, fundamentally, integrating dosimeters and can only provide information on the total dose delivered and do not qualify as real-time detectors. Hence, many groups are researching on the use of scintillator volumes or liquids to perform 3D dosimetry. Many of these dosimeters make use of the concept of tomography to reconstruct a 3D scintillation volume from many different perspectives (or projections). [15][16][17] While these techniques have shown potential to perform accurate dosimetry on a whole 3D volume, their application to dynamic, modulated radiation delivery is not straightforward, since components of the dosimeter need to rotate during acquisition, a fact which limits their real-time capability. Some other approaches use only fixed components to image and reconstruct a scintillating volume in real-time. 18,19 However, these techniques are usually restrained to very simple fields, such as intensitymodulated proton therapy or square photon fields, as a limited number of perspectives are used to image the scintillation volume of interest.
In this work, we investigated the possibility of conducting real-time 3D dosimetry on a volume of plastic scintillator using the light acquisition from a plenoptic camera. The main advantage of plenoptic cameras over conventional ones is that the light rays integrated by each pixel of the camera sensor originate from a known and restrained angle range. Thus, with the constraint of a beam's eye view projection of the deposited dose, the light integrated by each plenoptic camera sensor pixel can be back-projected into the scintillator volume using reconstruction algorithms similar to those developed for tomodosimetry, 20 computed tomodensitometry (CT), or positron emission tomography (PET). This kind of methodology allows for water-equivalent, high resolution, real-time, and full 3D dosimetry using only static detectors and phantoms.
We designed a 3D tomodosimeter based on the plenoptic camera technology to reconstruct the incident dose distribution on a volume of plastic scintillator as a function of delivery time. In this study, we developed the background method of 3D light reconstruction from a plenoptic camera acquisition and evaluated experimentally the detector's performance in the measurement of IMRT and VMAT clinical plans.

2.A. Plenoptic camera
A plenoptic camera, also referred to as a light-field or a single lens 3D camera, consists of a conventional camera having an array of microlenses coupled to its digital image sensor. In photography or microscopy, this additional optical component allows an adjustment of the focus or the point of view of an image after it has been recorded. [21][22][23][24][25] Figure 1 compares the ray tracings for a conventional and a plenoptic camera. The sensor pixels of a conventional camera solely record spatial information of incident photons: each sensor pixel is produced from the sum of all light rays striking the pixel's position. On the other hand, in a plenoptic camera, each microlens covers several sensor pixels, thus separating the light rays passing through the camera's main lens into smaller images recorded by the sensor. In this case, a sensor pixel represents a sum of light rays over a small portion of incident angles, providing spatial and directional information of incident photons. However, a point in the object side of the FIG. 1. Comparison of the ray tracings for a plenoptic (A) and a conventional camera (B). A conventional camera's sensor pixel intensity represents a sum of light rays over all incident angles (dotted lines). A plenoptic camera's sensor pixel intensity is produced by a sum of light rays from a small, known range of incident angles (shaded region). (a) Camera sensor; (b) microlens array; (c) main lens; and (d) scintillator volume. main lens can be imaged in more than one sensor pixel, thus reducing the spatial resolution of the system. 21 Different types of plenoptic cameras differ mainly by the distance between the array of microlenses and their sensor. "Standard" plenoptic cameras set this distance to the focal length of one microlens, which is usually around 500 μm. This configuration has the advantage of providing maximal angular resolution of the captured light-field at the expense of a very coarse spatial resolution. 21 In a "focused" plenoptic camera, this distance is made smaller (for a virtual main lens image) or larger (for a real main lens image) than the microlens focal length. 23 This type of configuration allows for a compromise between angular resolution and spatial resolution, maintaining the latter sufficiently fine for photography standards.
The purpose of a plenoptic camera in 3D dosimetry is based on using the angular information gathered from an imaged scintillator volume to reconstruct the 3D emitted light pattern. Using the spatial and directional information recorded by the plenoptic camera, an iterative, expectation maximization algorithm back-projects the photons received by the sensor pixels into the scintillating volume. The 3D dose distribution inside the scintillator volume can then be obtained from the 3D light pattern, as the light yield of organic scintillator is proportional to the energy deposited in its volume, independently of dose rate or energy in the therapeutic photon energy range. 26

2.B. Dose integral calculation
The relationship between a plenoptic camera pixel intensity and the dose deposited in each region of a scintillator volume is in essence very similar to the relationship between the light yield of a long scintillating optical fiber and the dose profile incident on the fiber. 16,20,27 Mathematically speaking one has: where p = {p j : j ∈ {1, . . . , n}} is the vector of weighted dose integrals, n is the number of plenoptic camera sensor pixels used, D = {D i : i ∈ {1, . . . , m}} is the vector of threedimensional dose distribution, m is the number of dose voxels in the scintillator volume, and A is the n × m projection matrix. Each element A j,i of this projection matrix corresponds to the contribution of the ith dose voxel to the jth weighted dose integral. Details on the calculation of A are provided in Sec. 2.C. The weighted dose integral (p j ) for each plenoptic camera sensor pixel is computed from its collected intensity (I j ) using the formalism developed for 2D tomodosimetry: 20 where the indices "ref" refer to a normalization field of known 3D dose distribution, D ref represents the vector of the threedimensional dose distribution of the normalization field and A j, * represents the jth row of the projection matrix A of Eq. (1a). As was the case for 2D tomodosimetry, D ref should ideally be known to a good precision (e.g., using a square or near-square radiation field) since an error on D ref could bias the computation of the weighted dose integrals, and thus the results of the 3D dosimeter. Additionally, the normalization field should be wide enough as to enable the simultaneous irradiation of the whole scintillator volume. The main difference between the weighted dose integrals calculated for a scintillating optical fiber and those calculated for a volume of scintillator is that, for the latter case, optical artifacts affect the propagation of light inside the scintillator and in the camera lens system. While most of these optical artifacts can be taken into account directly by the projection matrix A (see Sec. 2.C), optical scattering inside the scintillator volume is better modeled using a direct correction of the weighted dose integrals.
Two types of optical scattering were considered in this work. First order optical scattering represents the capacity of the scintillator and lens system to resolve a point light source inside the scintillator. If the magnitude of this effect is large in comparison to the system spatial resolution, its behavior can be modelled by the point-spread function (PSF) of the dosimeter system, and will affect the calculation of the weighted dose integrals as follows: where the PSF is usually experimentally measured for a given scintillator and lens system. Second order optical scattering represents the amount of light emitted from the whole scintillator volume by a point light source as a result of multiple absorption and reemission processes inside the scintillator. Moreover, due to the high refraction indices of scintillator plastic, light rays emitted inside the scintillator volume will be more likely to undergo total internal reflections, hence increasing the amount of light scattered far away from the emission point. Such behavior is modelled in this work as a homogeneous background of scattered scintillation light which is proportional to the total amount of dose deposited inside the scintillator volume. Mathematically speaking, one has where p j,c represents the corrected jth dose integral; β and β ref are proportionality constants; and the summation on each row of the projection matrix ( i A j,i ) represents the geometric area of the scintillator contributing to the jth pixel. Note that since the dose inside the scintillator is usually not known during the dose integral calculation, the total amount of dose deposited inside the scintillator volume can be approximated as proportional to the total amount of light collected from the scintillator volume by all pixels of the CCD camera: where α and α ref are proportionality constants. Note that in the presence of a non-negligible first order optical scattering, D ref should also be convolved using the measured PSF of Eq. (2a).

2.C. Computation of the projection matrix
Each row of the projection matrix A represents the contribution of all dose voxels to the corresponding camera sensor pixel weighted dose projection. In this work, the computation of A was performed using paraxial ray tracing, 28 i.e., light rays were propagated, starting from the camera sensor pixel, through the microlenses and main lens, then finally inside the scintillator volume. The contribution of any given voxel to the weighted dose projection under analysis is then proportional to the path length of the light rays inside this voxel. More specifically, each row of the projection matrix A was calculated taking into account: r The size of the sensor pixel and all possible ray angles between the sensor pixel and its corresponding microlens; r The focal length of each of the microlenses and the main lenses, as well as the optical distance between each optical element; r The total internal reflection at each side of the scintillator volume; r The partial reflection at the back of the scintillator volume, using Fresnel equations; 28 r The optical attenuation occurring inside the scintillator volume, following the exponential decay law. More specifically, the paraxial ray tracing was weighted by an exponential factor e −λ·x sc , where x sc represents the total distance travelled inside the scintillator volume by the considered ray.
The paraxial approximation was used for the ray tracing to optimize computational speed and also because the ray angles encountered were smaller than ±5 • in our configuration. No visible lens aberrations (e.g., distortion, astigmatism, coma) were observed. Lens aberrations were thus deemed negligible, although their effects could also be incorporated into the calculations (e.g., by using nonparaxial ray tracing). Vignetting was also neglected in our calculation of A, since the illumination falloff for pixels far from the central optical axis was the same for the normalization field and the field under analysis, thus cancelling out when applying Eq. (1b) or (2c).

2.D. Reconstruction algorithm
Similarly to previous work on tomodosimetry, Eq. (1a) was solved using an expectation maximization (EM) algorithm. 16,20 In addition to the weighted dose projection acquired using the plenoptic camera, a relative 2D beam's eye view (BEV) projection of the incident dose distribution was measured to guide the 3D dose reconstruction process. Finally, an additional hypothesis on the relationship between the dose and the depth inside the scintillator was incorporated in the reconstruction algorithm. More specifically, it was assumed that the dose at any point in the scintillator was a smooth, decreasing function with respect to the distance to the radiation source (z s ), also taking into account the z s −2 geometric decrease on the dose.
Thus, based on our previously developed method of tomographic dose reconstruction, 16,20 each iteration of the reconstruction algorithm was split into three successive parts: (a) The constraints of Eq. (1a) were satisfied using the EM algorithm: where D (k−1) i represents the value of D i after the three parts of the previous iteration and D (k) i,EM represents the value of D i at the kth iteration of the algorithm after the application of the EM algorithm; (b) The reconstructed dose was constrained by the relative shape of the BEV projection. As seen in Fig. 2, the 2D BEV projection (G) is scaled at the height of each pixel of interest (z i ) and its magnitude is adjusted using D (k) EM from step (a). Mathematically speaking, one has where δ represents the discrete Dirac delta function (which equals 1 only if its argument is zero); x i , y i , and z i represent the spatial coordinate of the ith dose pixel in the radiation source frame of reference; z G represents the perpendicular distance from the source to the BEV projection plane and D (k) i,BEV represents the value of D i at the kth iteration of the algorithm after the application of the BEV projection constraints; (c) The reconstructed dose was constrained by the expected dose vs depth relationship inside the scintillator. This relationship hypothesizes a square law decay of the dose with respect to the radiation source perpendicular distance (z), convolved with a linear relationship between the remaining quantity (analog to the tissuemaximum ratio) and z. To realize this constraint, as seen on Fig. 2, a line passing through the radiation source and each pixel of interest is first plotted. Each value on this line is obtained from the trilinear interpolation of the D (k) i,BEV values of step (b). After applying a z −2 factor on the obtained data to account for the square law decay of the dose, the remainder is then linearly fitted as a function of z. The obtained fit parameters are then used to evaluate the dose at the pixel of interest. Mathematically speaking, one has where D 0,i and a i are the aforementioned linear fit parameters. After an expectation maximization step using the acquired dose projections from the camera, the dose at the plane located at a distance z from the source inside the scintillator volume (large square) is constrained using the relative dose distribution acquired using the beam's eye view projection (G) (left). Then, the dose at the ith dose voxel (small square) is optimized assuming a smooth, decreasing dose with respect to the distance to the radiation source (right).
The stopping criterion for the algorithm was based on the pixel-by-pixel absolute dose variation ( D): where V s represents the total scintillator volume. Namely, the reconstruction process was terminated when D went below a fixed threshold, set in this work at 10 −9 cGy per cubic mm of reconstructed volume. It should be noted that the initial dose distribution in the algorithm ( D (0) ), while arbitrary, should be nonzero. In this work, a homogeneous dose distribution of 0.01 cGy was chosen as D (0) .

2.E. Experimental prototype
The prototype presented in this work is composed of a 10 × 10 × 10 cm 3 green plastic scintillator cube (EJ-260; density, 1.023 g/cc; Eljen technology, Sweetwater, TX) embedded in the center of a 10 cm radius by 17 cm long truncated cylindrical acrylic phantom (see Fig. 3). The scintillator was imaged with a "focused" plenoptic color camera (Raytrix R5, Raytrix GmbH, Kiel, Germany) using a f/2.9, 54 mm focal length lens (JML Optical Industries, Rochester, NY) and a 90 • mirror to enable lead shielding of the camera during irradiation. The scintillator was covered on all sides except the one facing the plenoptic camera by a black cardboard and its center was placed at the isocentre of a Varian Trilogy linear accelerator (Varian Medical Systems, Palo Alto, CA). Light acquisition using the plenoptic camera was performed at a rate of one image per second, and each image recorded by the plenoptic camera was treated for transient noise coming FIG. 3. Experimental prototype used in this work. A 10 × 10 × 10 cm 3 cube of green plastic scintillator (A) was inserted inside a cylindrical acrylic phantom (B) and imaged using a plenoptic camera (C) with the help of a 90 • mirror (D) to enable radiation shielding of the camera. A beam's eye view projection of the incident radiation field was also acquired using the linac portal imager (E). from stray radiation. 29 Only the green pixels of the camera sensor were used in this work.
A beam's eye view projection of the incident radiation field was obtained using the raw images acquired using the linear accelerator portal imaging device. The imager was operated in the continuous acquisition mode at a rate of approximately three acquisitions per second. Each acquisition was corrected for imager sagging and was convolved using a Gaussian kernel to better model the penumbra region of each acquired field, as suggested in earlier publication on IMRT field verification using the EPID. 30 The 3D dose inside the scintillator volume was reconstructed at a resolution of 2 × 2 × 2 mm 3 using the weighted dose projection calculated as in Sec. 2.B and the reconstruction algorithm described in Sec. 2.E. An 11 × 11 cm 2 square radiation field was used as the normalization field in weighted dose projection calculations performed using Eqs. (1b) and (2c). D ref was obtained for this normalization field using calculation performed by the treatment planning system (TPS) (Pinnacle 3 9.2, Philips Healthcare, Andover, MA) on a CT scan of the previously described phantom. An optical index of 1.58 and an attenuation length of λ = 1/150 cm −1 was used in the calculation of the projection matrix (see Sec. 2.C), according to the scintillator manufacturer specifications.
The linac gantry angle, needed in the reconstruction algorithm (see Sec. 2.D), was determined directly from the linac interface for static gantry delivery, and calculated from the camera image acquisition in the case of dynamic gantry delivery. In the latter case, this information was determined using a maximum gradient search method on the intensity image recorded by the plenoptic camera.

2.F. Field acquisition
To evaluate the PSF of the scintillator and imaging system, an orthovoltage x-ray therapy system (Xstrahl200, Xstrahl Ltd., Camberley, UK) was used as a 180 kV photon radiation source to irradiate the cube. Two lead sheets of approximately 2 mm thickness each were used to create a scintillation edge at the center of the cube by shielding half of the top of the scintillator. The sheets were placed perpendicular to the central beam axis and the obtained edge-spread function (ESF) was then numerically differentiated to obtain a linespread function (LSF), which was then fitted with a Gaussian curve to estimate the width of the PSF.
Parameters α and α ref of Eq. (2c) were estimated using a 5 × 5 cm 2 6 MV square photon field produced in the linear accelerator configuration of Sec. 2.E. Measured weighted dose projections calculated using Eq. (1b) were compared to the expected dose projections (p j,e ) calculated using the following formula: where D e represents the expected dose distribution of the 5 × 5 cm 2 field as calculated with the TPS on the phantom CT. Parameters α and α ref were automatically adjusted until the difference between the weighted dose projections calculated from Eqs. (2c) and (7) was minimal.
Finally, two treatment plans from two clinical brain tumor cases were acquired using the dosimeter prototype. The first was a seven incidences step-and-shoot IMRT plan and the second was a 360 • rotation SmartArc VMAT delivery, each delivered with 6MV photons. For each case, the reconstructed dose at a resolution of 2 × 2 × 2 mm 3 was compared to the expected dose calculated using the Pinnacle 3 TPS. Each delivery integrated dose was also compared to a radiochromic film acquisition. To do so for both treatment plans, the scintillator cube was removed from its acrylic phantom and replaced by three 11 × 10 cm 2 sheets of radiochromic film (Gafchromic EBT3, Ashland Inc., Covington, KY) immersed in water in a Lat/AP orientation. Each film was scanned in reflection mode using a flatbed scanner (Epson 10000XL, Epson Canada Ltd., Markham, Canada) before and 40 h after irradiation. 31 Optical densities measured in the red channel were converted to dose with the help of a calibration curve. The latter was obtained using small squares of the same film batch irradiated inside a solid water phantom at the same linac used for irradiation.

3.A. Optical scattering characterization
The full width half maximum (FWHM) of the measured PSF was of (0.5 ± 0.2) mm. As this width is small with respect to the length of each dose voxel (2 mm), the PSF was approximated by a Dirac delta function and Eq. (1a) was used instead of Eq. (2a) in the dose reconstruction process. The optimal value for the parameters α and α ref of Eq. (2c) were found to be 4.0 × 10 −6 and 4.2×10 −6 , respectively. The effect of the second order optical scattering correction on the weighted dose projections is presented in Fig. 4. Without second order optical scattering correction, the weighted dose projections are overestimated in the low-light region of the scintillator, while they are underestimated in the high-light region of the scintillator.

3.B. 3D dose comparison with the TPS
The 3D dose reconstructed inside the scintillator volume for the whole IMRT and VMAT plan are presented in Figs. 5 and 6, respectively. Three orthogonal 2D views of the reconstructed dose are presented for each plan, along with their absolute dose deviation from the expected TPS dose. Dose accuracy for the IMRT field was overall better than 4% of the maximum dose except for regions of high dose gradient. For the VMAT plan the accuracy was overall better than 3% of the maximum dose outside of high dose gradient regions.
The gamma test success rates and mean absolute dose difference for each reconstructed dose distributions are summarized in Table I. Each gamma test was performed using a linearly interpolated TPS dose distribution at a resolution of 0.1×0.1 mm 2 and a dose difference threshold expressed as a percentage of the maximum TPS dose. For all but two IMRT incidences, 3%/3 mm, 2%/2 mm, and 3%/1 mm gamma test success rates were greater or equal to 98.1%, 90.0%, and 93.7% for the pixels receiving over 10% of the maximum dose, respectively. Dose difference statistics were computed over the low gradient (<1% D max /mm) region of each field.
The mean absolute dose differences calculated over this region in the high dose (D > 50% D max ) and low dose (50% > D/ D max > 10%) regions of the field were under 2% of the maximum dose for all but two incidences of the step-andshoot IMRT plan.

3.C. Dose comparison with radiochromic films
Comparison between the reconstructed dose in the scintillator volume and the dose measured using radiochromic films is presented in Figs. 7 and 8 for the whole IMRT and VMAT treatment, respectively. For the IMRT plan, dose differences were overall less than 5% outside of high dose gradi-ent regions. A mean absolute dose difference of 2.1% was obtained in the high dose (D > 50% D max ), low gradient (<1% D max /mm) region of the field and a mean absolute dose difference of 1.6% in the low dose (D < 50% D max ), low gradient (<1% D max /mm) region of the field. For the VMAT plan, dose differences were overall less than 3% outside of high dose gradient regions. A mean absolute dose difference of 1. 14%   TABLE I. Gamma test success rates and mean gamma computed between the reconstructed dose planes and the predicted dose planes from the Pinnacle 3 treatment planning system for a seven incidence step-and-shoot IMRT and an VMAT rotational plan. a Calculated for each field over the region with dose higher than 50% D max and gradient lower than 1% D max /mm. b Calculated for each field over the region with dose lower than 50% D max and gradient lower than 1% D max /mm. was obtained in the high dose (D > 50% D max ), low gradient (<1% D max /mm) region of the field and a mean absolute dose difference of 1.15% in the low dose (D < 50% D max ), low gradient (<1% D max /mm) region of the field.

4.A. Quality of the 3D dosimetry
The presented 3D dosimeter is capable of performing realtime, full 3D dosimetry of intensity-modulated plans at a resolution of 2 mm, which is, to the best of our knowledge, a first in the field of experimental dosimetry. In addition to the novelty of this approach, it also shows good accuracy for both IMRT and VMAT dose distributions. Namely, deviations be-tween the reconstructed 3D dose distribution and the expected dose from the treatment planning system Pinnacle 3 are typically within 3% of the maximum dose in the low dose gradient region of each field, as seen in Table I and Figs. 5 and 6. It should be also be noted that, although not explicitly visible in the presented results, the 3D dose distribution was obtained as a function of delivery time, a desirable feature for dynamic treatment experimental verification.
The benefit of using a plenoptic camera in this methodology comes from the fact that each pixel of the camera sensor is associated with a narrower spatial region of the 3D scintillator volume than in the case of a conventional camera, as seen in Fig. 1. As an example, while a conventional sensor pixel can integrate over a 10 • range of incident light ray angles, the same sensor pixel in our plenoptic configuration covers less than a 1 • range of angles. The loss in spatial resolution of the plenoptic camera image is not of particular importance, since the reduced spatial resolution remains under 0.2 × 0.2 mm 2 per sensor pixel at the point of focus, which is largely adequate for application in millimeter resolution dosimetry.
This level of dose accuracy could not be obtained using only the measured dose projections without some kind of orthogonal dose projection information. Indeed, similarly to what was observed with simulations in cylindrical tomodosimetry, 16 the light impinging on the camera comes from a limited set of relatively small angles. This fact makes it impossible to resolve spatial details in the direction of increasing distance relative to the camera sensor. As the range of ray angles acquired using a single, fixed camera are very limited, the portal imager of the linear accelerator was used as a beam's eye view dose projection to improve the spatial resolution of the reconstructed 3D dose distribution, especially in the direction parallel to the main optical axis of the camera. The choice of the portal imager was made mainly because it is, by construction, always approximately perpendicular to the radiation field incidence, and it is also readily available on most linear accelerators. Without loss of generality and/or applicability of the proposed method, a second camera, plenoptic or conventional, could also have been used to acquire a beam's eye view projection by imaging the scintillator volume from the radiation source's point of view.
The resolution of the presented prototype, 2 mm in each spatial direction, is higher than most dosimetric gel systems 10 but is still considered adequate for IMRT and VMAT QA applications. Indeed, Fourier analysis studies concluded that a dose grid resolution of 2.5 mm or lower would be sufficient to detect dose errors larger than a percent for IMRT cases. 32 While this conclusion was made for 2D dose distributions, it is expected to extend to three-dimensional measurements. It should be noted that the well demonstrated properties of plastic scintillators 26 are also expected to apply to the 3D dosimetry using the present methodology. As such, the presented prototype is expected to show no dose or dose-rate dependence. Moreover, since most camera sensors have an adjustable acquisition time, no maximum usable dose range is foreseen. The time needed for each 3D dose reconstruction was of approximately 3 min per plenoptic camera acquisition on a 2.67 GHz CPU using a MATLAB R implementation, but is expected to require only a few seconds per reconstruction using GPU implementation of the algorithm. 33 It should also be noted that the reconstruction time is expected to be proportional to the number of scintillator volume voxels and the number of camera sensor pixels. Stability of the response over time will mostly depend on the variation over time of the scintillator optical attenuation λ, which usually tends to increase with radiation damage in plastics. Although the optical attenuation is of minor impact in the presented methodology, any variations are expected to be accounted for using detector recalibration. Signal contamination by Cerenkov emission was evaluated to less than 0.5% of the collected scintillation light for all fields investigated, and therefore Cerenkov light filtration was deemed unnecessary.

4.B. Comparison with radiochromic films
To account for possible errors in the TPS, a verification of the delivered dose on the scintillator was made using radiochromic films. Once again, the agreement between the 3D reconstructed dose and the radiochromic films were overall better than 5% for the IMRT delivery, and 3% for the lower gradient VMAT delivery, as seen in Figs. 7 and 8, respectively.
The higher dose differences between the radiochromic films and the 3D reconstructed dose than those of the previous comparison with the Pinnacle 3 TPS are assumed to come from: (a) the high uncertainties of radiochromic films reading and processing, especially at low doses; 31 (b) the fact that the phantom was not 100% identical for the beam delivery on the scintillator volume and the radiochromic film. Indeed, in the latter case, the scintillator phantom cavity section not occupied by radiochromic films was filled with water. Although the scintillator itself is water-equivalent and a subsequent CT scan showed no important difference between the two phantom configurations, some dose inaccuracies could have been noticeable at the border of the scintillator cavity, where a small air gap was present in the scintillator configuration but not in that of the water fill-in; (c) some residual registration error between the scintillator and the 3D reconstructed dose, especially for the IMRT cases where higher dose gradients were observable.

4.C. Optical scattering and other current limitations
Since the spatial magnitude of first order optical scattering (approximately 0.5 mm FWHM) is small with respect to the system's current spatial resolution (2 mm in each direction), we concluded that this effect need not to be considered significant, and thus Eq. (1a) was used instead of Eq. (2a) in the dose reconstruction process. However, as the spatial resolution of the 3D dose distribution is a free parameter in the reconstruction algorithm, this effect may be of interest in future work if the spatial resolution is to be improved.
The impact of second order optical scattering is however more important than first order scattering in the presented methodology. Indeed, without second order scattering corrections, the weighted dose projections can suffer from large relative errors, especially in regions of low light level (i.e., in the out-of-field regions), as seen in Fig. 4. However, such errors have only a small effect on the spatial resolution of the 3D reconstructed dose, as the latter is also constrained by the BEV projection. In our particular case, the BEV projection was free of second order optical scattering consideration, since the portal imager is sensitive to ionizing radiation only. Errors on the weighted dose projections caused by second order optical scattering will however largely affect the reconstructed dose output and should be corrected for optimal accuracy. The correction proposed in this work makes the assumption that the second order optical scattering magnitude is evenly distributed throughout the scintillator volume, and thus corrects each weighted dose projection of a given acquisition by the same amount [see Eq. (2c)]. A 5 × 5 cm 2 field was selected to determine the α and α ref constants of Eq. (2c) mainly because it allowed the clear visualization of both low-light and high-light region in the scintillator (see region A and B in Fig. 4). It should be noted that the value of α ref employed for the reference field correction was 5% higher than the value of α used every other acquired field correction. This difference was mainly attributable to the fact that only the reference field used in this work (11 × 11 cm 2 ) was wide enough as to enable the irradiation of the whole scintillator cube, including its edges and external faces. It is hypothesized that imperfections at these edges and faces may have caused more second order scattering per deposited dose than for smaller radiation field. While second order optical scattering correction works well for most cases, as suggested by the results presented in Sec. 3.A, it may however be considered simplistic; thus, a more detailed model of second order scattering production inside the scintillator may lead to even more accurate results.
Apart from optical scattering, the main source of inaccuracy in the presented 3D dosimeter comes from the BEV projections itself. Indeed, while this projection mainly guides the 3D reconstruction by providing spatial information of the incident 3D dose pattern, it also influences the amount of out-offield dose relative to the in-field dose recorded by the dosimeter. The choice of the portal imager to perform the BEV projection in this work was based on its availability and ease of use; however, it is not considered to be very accurate in the low dose regions without extensive correction factors such as nonuniform backscatter and flatness corrections. 30 An alternative to the portal dosimetry solution would be to use another high resolution 2D detector for the BEV projection, or even to image the scintillator volume from the source point of view using a mirror system or a second camera, either plenoptic or other.
Finally, the acquisition rate of the plenoptic camera is expected to influence the dose reconstruction process of dynamic arc delivery in the proposed method. Indeed, the linac gantry is considered static in the BEV and dose vs depth constraints of the presented algorithm (Eqs. (4) and (5), respectively), which is a valid approximation only if the acquisition time is small with respect to the gantry rotation speed. However, the shorter the camera acquisition time, the noisier each camera pixel data will be and the longer it will take to reconstruct the whole 3D dose for the entire treatment. The acquisition rate of one image per second was considered in this paper as a compromise between the abovementioned effects, but could be reduced in the future using, for example, a high sensitivity CCD sensor. 17,18

CONCLUSION
We developed a 3D dosimeter that allows the iterative reconstruction of full three-dimensional dose distributions in real-time using plenoptic camera technology. Dose reconstructions of both IMRT and VMAT fields showed promising results when compared to the planned dose calculated on the phantom geometry by the treatment planning software, and are comparable to the results obtained using radiochromic film acquisition. The presented prototype achieved waterequivalent, 2 × 2 × 2 mm 3 resolution dosimetry on a 1000 cm 3 volume using only a fixed plenoptic camera and the linac portal imaging device. We expect this kind of dosimeter to be very promising in dynamic treatment quality assurance, especially for rotational treatment such as VMAT and even stereotactic treatment delivery validation. The same approach could also be used for proton and heavy-ion beams as long as correction for scintillation quenching is provided.