# A tomographic reconstruction algorithm for cross-sectional imaging of IMRT beams from six projections

Odran Pivot, P. Pittet, Rolf Clackdoyle, Laurent Desbat, Simon Rit

► To cite this version:

Odran Pivot, P. Pittet, Rolf Clackdoyle, Laurent Desbat, Simon Rit. A tomographic reconstruction algorithm for cross-sectional imaging of IMRT beams from six projections. Physics in Medicine and Biology, 2023, 68 (9), pp.095015. 10.1088/1361-6560/acc925 hal-04061147v2

HAL Id: hal-04061147

https://hal.science/hal-04061147v2

Submitted on 27 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A tomographic reconstruction algorithm for
cross-sectional imaging of IMRT beams from six
projections

Odran Pivot<fnref n="1" />, Patrick Pittet<fnref n="2" />, Rolf Clackdoyle<fnref n="1" />, Laurent Desbat<fnref n="1" />, and Simon Rit<fnref n="3" />

E-mail: odran.pivot@protonmail.com

February 2023

**Abstract.** Objective: Patient-specific Quality Assurance (QA) measurements are of key importance in radiotherapy for safe and efficient treatment delivery and allow early detection of clinically relevant errors. Such QA processes remain challenging to implement for complex IMRT radiotherapy fields delivered using a multileaf collimator (MLC) which often feature small open segments and raise QA issues similar to those encountered in small field dosimetry. Recently, detectors based on long scintillating fibers have been proposed to measure a few parallel projections of the irradiation field with good performance for small field dosimetry. The purpose of this work is to develop and validate a novel approach to reconstruct MLC-shaped small irradiation fields from six projections.

**Approach:** The proposed field reconstruction method uses a limited number of geometric parameters to model the irradiation field. These parameters are iteratively estimated with a steepest descent algorithm. The reconstruction method was first validated on simulated data. Real data were measured with a water-equivalent slab phantom equipped with a detector made of 6 scintillating-fiber ribbons placed at 1 m from the source. A radiochromic film was used to acquire a reference measurement of a first dose distribution in the slab phantom at the same source-to-detector distance and the treatment planning system (TPS) provided the reference for another dose distribution. In addition, simulated errors introduced on the delivered dose, field location and field shape were used to evaluate the ability of the proposed method to efficiently identify a deviation between the planned and delivered treatments.

**Main results:** For a first small IMRT segment, 3%/3 mm, 2%/2 mm and 2%/1 mm gamma analysis conducted between the reconstructed dose distribution and the dose measured with radiochromic film exhibited pass rates of 100 %, 99.9 % and 95.7 %, respectively. For a second and smaller IMRT segment, the same gamma analysis performed between the reconstructed dose distribution and the reference provided by the TPS showed pass rates of 100 %, 99.4 % and 92.6 % for the 3%/3 mm, 2%/2 mm and 2%/1 mm gamma criteria, respectively. Gamma analysis of the simulated treatment

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France.</footnote>

<footnote n="2">² Institut des Nanotechnologies de Lyon (INL), CNRS UMR5270, Université de Lyon, Université Claude Bernard Lyon 1, F-69100, France.</footnote>

<footnote n="3">³ Univ. Lyon, INSA-Lyon, UCB Lyon 1, UJM Saint-Etienne, CNRS, Inserm, CREATIS, Centre Lyon Bérard, 69373 Lyon, France.</footnote>
delivery errors showed the ability of the reconstruction algorithm to detect a 3 %
deviation between the planned and delivered doses, as well as shifts lower than 7 mm
and 3 mm when considering an individual leaf and a whole field shift, respectively.

Significance: The proposed method allows accurate tomographic reconstruction of IMRT segments by processing projections measured with six scintillating-fiber ribbons and is suitable for water-equivalent real-time small IMRT segments QA.

*Keywords*: tomographic reconstruction, radiotherapy quality assurance, scintillating-fiber detector, multileaf collimator

Submitted to: *Phys. Med. Biol.*
# 1. Introduction

Modern radiotherapy techniques implementing Intensity Modulated Radiation Therapy (IMRT) use multileaf collimators (MLC) to shape complex irradiation fields [1, 2, 3]. IMRT Quality Assurance (QA) is an important process for safe and efficient patient treatment because it allows early detection of clinically relevant errors in the treatment delivery. In this context, QA measurements ideally use tissue-equivalent detectors with a negligible volume effect. The radiochromic film characteristics fit these requirements and films are widely adopted for small field QA measurements [4]. However, they must be manually scanned and are not suited to automatic QA before treatment fractions (i.e. single patient-specific QA test prior to plan delivery). In the past few decades, many strategies have been investigated for achieving accurate and real-time measurement-based IMRT QA [5]. Simple QA strategies use the log files obtained from the radiotherapy device for real-time identification of the MLC leaf positions [6], but it is not suitable for an independent assessment of the dose distribution because it is based on information provided by the LINAC manufacturer and not on independent measurements. Electronic portal imaging devices (EPIDs) can be considered for IMRT field QA since they provide high-resolution images. However, they show a non-water-equivalent dosimetric response caused by the high-Z material in some layer of the detector, which results in a higher sensitivity to lower-energy photons [7], and thus require a conversion procedure to obtain 2D dose distributions. High spatial resolution diode arrays have been proposed for stereotactic QA [8] and can also be used for QA of MLC-shaped irradiation fields. However, they have a millimeter detector spacing (e.g. 2.47 mm for the QA system studied in [8]), which could be insufficient in case of small open segments. Moreover, diodes are made of silicon, which is not a tissue-equivalent material. Organic scintillating fibers which have the intrinsic advantage to be made of tissue-equivalent material can be considered for real-time QA measurements of small fields. However, they measure the integrated dose along the fiber axis and are therefore intrinsically 1D detectors. With a fiber ribbon, one can measure the parallel projection of the field along the ribbon axis, but not the 2D dose distribution directly.

By employing tomographic reconstruction, 2D dosimetry based on scintillating fibers has been reported from 18 sequentially acquired parallel projections of the irradiation field [9]. More recently, a system with 6 stacked SciFi detectors has also been developed to provide six projections acquired simultaneously (i.e. on the same beam irradiation) [10, 11]. Tomographic reconstruction methods adapted to simple irradiation field shapes (e.g. an arbitrary disk [12], or an arbitrary rectangle [13, 11]) have been proposed and validated on this QA system. In the case of complex MLC-shaped fields, the small number of projections measured with such a QA system presents a challenge for robust and accurate irradiation field reconstruction. In this article, we describe the development, implementation, and preliminary validation of a suitable algorithm for reconstruction of MLC-shaped irradiation fields from six parallel projections.

Over the last twenty years, there have been various approaches to performing
tomographic reconstruction from few projections. Today, in medical imaging, "few projections" would mean anything less than about 50 projections. For general imaging purposes, the number of projections needed to recover a function with specified resolution features has been analyzed in detail (see chapter 3 of [14]). When fewer projections are available, reliable reconstruction is not possible unless the class of objects under study is restricted. In 2002, an iterative method was proposed to minimize the $L^p$ norm of the image with $p = 1.1$ (rather than the usual Euclidean $p = 2$) subject to matching the projection data [15]. In an underdetermined system with few projections, this choice of $L^p$ norm heavily favours solutions that have the most possible zero values. Simulations of sparse blood-vessel trees were shown to be reconstructible with as few as 8 projections. However, the irradiation fields under consideration cannot be considered to be sparse objects. Extending the class of objects to include large constant areas (rather than only zero areas), an alternating total variation (TV) minimization with ART reconstruction [16] (Algebraic Reconstruction Technique) was proposed in 2006 [17], with illustrations of high-quality reconstructions of the Shepp Logan phantom using 20 fan-beam projections. In 2011, using a variation of this idea where the TV minimization was performed with a different algorithm, and alternated with SART (Simultaneous Algebraic Reconstruction Technique) instead of ART, simulations of a high-contrast Shepp Logan phantom were satisfactory with 21 projections but poor when only using 15 projections [18]. The more challenging Forbild phantom, even though it also has large constant regions, required 61 projections but did not perfectly resolve the small high-frequency features. For the situation of irradiation fields imaged using scintillating fibers, a similar alternating TV minimization with an EM algorithm gave good reconstructions from only 18 projections [9]. The "objects" in this case were irradiation fields which were nearly binary images, and without high frequency structure. If further restrictions are placed on the object functions, satisfactory reconstructions can be obtained with fewer projections. In 2011, the DART reconstruction algorithm (Discrete Algebraic Reconstruction Technique) was introduced, which is capable of reconstructing objects consisting of a small set of known densities [19]. The algorithm proceeds by alternating SART reconstructions with identifying and removing from the reconstruction problem the areas of known constant density. Over a range of object functions, it was observed that 8 to 10 projections were generally sufficient for highly satisfactory reconstructions. Further in this direction, if the object is a binary function (i.e. just one known density on a zero background), then the tools of geometric tomography [20] can be consulted. It is known that if the binary object is convex, then it can be unambiguously identified from any seven projection directions, or from four well-chosen projection directions. If the convexity requirement is dropped then the situation is more complex and depends on intermediate constraints, such as h-convexity ("horizontal convexity") whereby all horizontal cross-sections of a suitable rotation of the object are convex (in a one-dimensional sense).

Reliable reconstruction of an IMRT segment would require uniqueness of the solution (i.e. the measurements correspond to one unique MLC irradiation field). It
has been demonstrated that even for binary objects, uniqueness cannot be ensured but most of the non-unique situations can be avoided by using a detector with enough projection angles [21]. For six projections, non-uniqueness situations can arise with 32 or more leaf-pairs. This would be a limitation of using six-projection detectors for reconstructing large fields, but the assumption that the solution is unique is reasonable when handling less than 32 pairs of leaves. In this study, we only consider situations with less than 32 leaf-pairs and we assume for the rest of the article that the solution is unique.

All the previous methods listed above use a pixel-based model for the object function, but none of those approaches appear suitable for reconstruction from six irradiation field projections. The TV-based methods appear to break down at about 15 projections. The DART method requires knowledge of the binary image intensity, which might be possible to estimate but even that approach seems to require 8 to 10 projections. Satisfactory results with 6 projections could doubtless be achieved if the irradiation field were always a convex shape but this requirement is too restrictive for our purposes. Rather than using a pixel-based image, we suggest taking advantage of the highly structured irradiation field geometry. For example, two opposite leaves of an MLC depict a rectangle and an MLC irradiation field can therefore be decomposed as contiguous rectangles with the same orientation, height and vertical distribution as the leaves. We propose to use this prior information to model the forward problem as a function of a few spatial variables describing the rectangles formed by the leaf positions and, following [13], an isotropic gaussian kernel modeling the blurring penumbra. This description greatly reduces the number of unknown parameters, compared to a pixel model. In this work, we develop and test a custom algorithm for reconstructing MLC irradiation fields from six scintillating-fiber projections. We present simulation results for typical QA situations to show feasibility of this approach, although further experimental testing would be required to validate application to clinical data.

## 2. Materials and methods

### 2.1. The forward problem

#### 2.1.1. The irradiation field model
The $N$ pairs of leaves of an MLC describe a shape composed of $N$ rectangles, where the $l$-th rectangle is characterized by four parameters: its width $w_l \in \mathbb{R}$, its height $h \in \mathbb{R}$ and the coordinates of its center $(x_l, y_l) \in \mathbb{R}^2$ (see figure 1). We assume that the collimator plane is parallel to the plane of the scintillating-fiber detector (i.e., the gantry angle is 0), which allows for a simpler implementation of the proposed approach (with no additional transformation needed to recover the projected geometry of the MLC in the detector plane). Therefore, the rectangles all have the same height $h$, corresponding to the width of the MLC leaves projected at the isocenter (i.e. in the plane of the detector). In addition, they are regularly distributed and their vertical position can therefore be deduced from the first one with
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic of the shape described by an MLC. $(x_l, y_l)$ is the center of the $l^{th}$ rectangle, $w_l$ its width and $h$ is the height of all rectangles.</figcaption></figure>

$y_l = y_1 + (l - 1)h$, for $l = 2, 3, \dots, N$. Therefore, the shape described by an MLC can be characterized by the $2N + 2$ geometric variables $(x_1, w_1), (x_2, w_2), \dots, (x_N, w_N), y_1, h$. In the rest of this article, we assume that the rectangles height $h$ is known, and we denote $\mu = (w_1, w_2, \dots, w_N, x_1, x_2, \dots, x_N, y_1)^T \in \mathbb{R}^{2N+1}$, the vector constructed by concatenating the $2N + 1$ unknown geometric variables.

Physical phenomena occurring near the leaf borders induce a penumbra that we model as an isotropic two-dimensional gaussian smoothing function [13]. In the plane of the detector, the MLC irradiation field $f_\mu$ is then given by

$$f_{\mu} = d \sum_{l=1}^{N} \chi_{\mu,l} *_{2D} g, \qquad (1)$$

where $d \in \mathbb{R}$ is a constant dose, $*_{2D}$ the two-dimensional convolution operator, $g$ the standard two-dimensional gaussian function

$$g(x\vec{e}_x + y\vec{e}_y) = \frac{1}{2\pi\sigma^2} \exp(-\frac{x^2 + y^2}{2\sigma^2}), \qquad (2)$$

with a known standard deviation $\sigma \in \mathbb{R}$ and $\chi_{\mu,l}$ the indicator function of the $l$-th rectangle:

$$\chi_{\mu,l}(x\vec{e}_x + y\vec{e}_y) = \begin{cases} 1, & \text{if } x \in [x_l - \frac{w_l}{2}, x_l + \frac{w_l}{2}] \text{ and } y \in [y_l - \frac{h}{2}, y_l + \frac{h}{2}] \\ 0, & \text{otherwise.} \end{cases} \qquad (3)$$

In equations (2) and (3), $(\vec{e}_x, \vec{e}_y)$ is the orthonormal basis of the plane such that $\vec{e}_x$ is colinear with the direction of the displacement of the collimator leaves.
For a full description of the irradiation field model of equation (1), the 2N + 1 geometric variables in $\mu$ and the additional parameter $d$ are needed.

### 2.1.2. Measurement model
The measured signal provided by the scintillating-fiber detector is the integrated dose along the fibers, or, equivalently, a finite number of samples of the Radon transform of the irradiation field $f_\mu$, yielding the following measurement model:

$$p_{\mu}(\phi, s) = \int_{-\infty}^{\infty} f_{\mu}(r\vec{u}_{\phi} + s\vec{v}_{\phi}) dr = R_{\phi}f_{\mu}(s), \quad (4)$$

where $p_\mu(\phi, s)$ is the measurement at projection angle $\phi \in \Phi = \{\phi_1, \phi_2, \dots, \phi_P\}$ and pixel position $s$, $(\vec{u}_\phi, \vec{v}_\phi)$ is the orthonormal basis of the plane such that $\vec{u}_\phi$ is colinear with the direction of the fibers at projection angle $\phi$ and $R_\phi$ is the Radon transform operator at projection angle $\phi$. The set $\Phi$ of projection angles has a small number of elements $P$ (in this work, $P = 6$) and is predetermined by the detector geometry and the known orientation of the collimator. By incorporating equation (1) into equation (4) and taking advantage of the properties of the Radon transform (linearity and convolution), one can derive the following forward model:

$$p_{\mu}(\phi, s) = d \sum_{l=1}^{N} (R_{\phi}\chi_{\mu,l} *_{1D} R_{\phi}g)(s), \quad (5)$$

where $*_{1D}$ denotes the one-dimensional convolution operator.

The Radon transform of a unit intensity rectangle $R_\phi\chi_{\mu,l}$ appearing in equation (5) has the shape of a trapezoid depending on four variables $(m_l, s_l, \delta_l^-, \delta_l^+)^\in \mathbb{R}^4$:

$$R_{\phi}\chi_{\mu,l}(s) = \begin{cases} m_l, & \text{if } s \in [s_l - \delta_l^-, s_l + \delta_l^-] \\ \frac{m}{\delta_l^+ - \delta_l^-}(s - (s_l - \delta_l^+)), & \text{if } s \in [s_l - \delta_l^+, s_l - \delta_l^-] \\ -\frac{m}{\delta_l^+ - \delta_l^-}(s - (s_l + \delta_l^+)), & \text{if } s \in [s_l + \delta_l^-, s_l + \delta_l^+] \\ 0, & \text{otherwise} \end{cases}, \quad (6)$$

where $m_l = \min\left(\frac{h}{|\sin\phi|}, \frac{w_l}{|\cos\phi|}\right)$ is the maximum value of the projection of the unit intensity rectangle, $s_l = -x_l \sin\phi + y_l \cos\phi = (x_l, y_l) \cdot \vec{v}_\phi$ is the signed distance between the origin and the projection of $(x_l, y_l)$ onto the line in the direction $\vec{v}_\phi$ passing through the origin, and $\delta_l^- = \frac{1}{2}|w_l|\sin\phi - h|\cos\phi|$ and $\delta_l^+ = \frac{1}{2}|w_l|\sin\phi + h|\cos\phi|$ the absolute distances between the projection of $(x_l, y_l)$ and respectively the closest and the furthest rectangle corner projections, all on the line in the direction $\vec{v}_\phi$ passing through the origin. The notations are illustrated in figure 2.

The Radon transform $R_\phi g$ of the two-dimensional gaussian is a one-dimensional gaussian with the same standard deviation:

$$R_{\phi}g(s) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{s^2}{2\sigma^2}\right). \quad (7)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Diagram of the Radon transform $R_{\phi\chi_{\mu,l}}$ of the l-th unit intensity rectangle for $\phi = -40^{\circ}$.</figcaption></figure>

## 2.2. Irradiation field reconstruction

### 2.2.1. Dose from zero-order moment and geometric variables
We let $M_{\phi}^0$ denote the zero-order moment defined as the integral of the projection at angle $\phi$: $M_{\phi}^0 = \int_{-\infty}^{\infty} p_{\mu}(\phi, s)ds$. The moment $M_{\phi}^0$ is then the integral of the irradiation field along the whole plane and can therefore be expressed according to the dose and the geometric variables:

$$M_{\phi}^{0} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f_{\mu}(r\vec{u}_{\phi} + s\vec{v}_{\phi}) \, dr \, ds = \sum_{l=1}^{N} dw_l h = dh \sum_{l=1}^{N} w_l. \quad (8)$$

If the projections are noiseless, the zero-order moment does not depend on $\phi$ ($M_{\phi}^0 = M^0, \forall \phi \in \Phi$). In practice, when several projections are measured, $M^0$ can be robustly estimated from the average of the individual estimates $M_{\phi}^0$.

From equation (8), we deduce that the dose can be retrieved from the geometric variables:

$$d = \frac{M^0}{h \sum_{l=1}^{N} w_l}. \quad (9)$$

Finally, by incorporating equation (9) into equation (5), we obtain the following modified forward model:

$$p_{\mu}(\phi, s) = \frac{M^0}{h \sum_{l=1}^{N} w_l} \sum_{l=1}^{N} (R_{\phi\chi_{\mu,l}} *_{1D} R_{\phi g})(s), \quad (10)$$

where $R_{\phi\chi_{\mu,l}}$ and $R_{\phi g}$ are given by equations (6) and (7), respectively.

### 2.2.2. Least-squares estimation of geometric variables
Having expressed the dose as a function of the measurements and the geometric variables, only the $2N+1$ variables in $\mu$
remain unknown and we estimate them by solving the following least squares problem:

$$ \tilde{\mu} \in \arg \min_{\mu \in \mathbb{R}^{2N+1}} C(\mu), \qquad (11) $$

where

$$ C(\mu) = \sum_{\phi \in \Phi} \int_{-\infty}^{\infty} (p_{\mu}(\phi, s) - q(\phi, s))^2 ds \qquad (12) $$

is the cost function to be minimized ($p_{\mu}$ and $q$ being the parameterized forward model and the measurements, respectively).

Starting from an initial guess $\mu^{(0)}$ (which may be defined from the treatment plan), we iteratively estimate $\tilde{\mu}$ with the following steepest descent algorithm:

$$ \mu^{(n+1)} = \mu^{(n)} - \lambda^{(n)} \nabla C(\mu^{(n)}), \qquad (13) $$

where $\lambda^{(n)}$ is the (positive) step length at iterate $n$ computed using a backtracking line-search algorithm [22] and $\nabla C(\mu^{(n)}) \in \mathbb{R}^{2N+1}$ the cost function gradient at $\mu^{(n)}$ (details of its computation can be found in the Appendix). It is worth noting that the cost function of equation (12) is non-convex and may admit some local minima, and we assume that the initial guess is close enough to the solution (assumed to be unique) so that the descent algorithm converges to the global minimum.

The algorithm stops when the norm of the gradient has been reduced by a preset factor $\epsilon$, such that $||\nabla C(\mu^{(n)})||/||\nabla C(\mu^{(0)})|| < \epsilon$, and we set $\tilde{\mu} = \mu^{(n)}$.

## 2.3. Experiments

### 2.3.1. Testing the reconstruction algorithm with simulated data

**Simulated fields** We simulated three typical IMRT segments to evaluate the proposed reconstruction method in different situations. They were all implemented with $N = 10$ pairs of leaves (which is a typical number of leaf-pairs for shaping small IMRT segments) with $h = 2.5$ mm (mimicking the commercial Varian HD120 multileaf collimator), illustrating situations of varying difficulty. The dimension of the optimization problem to reconstruct these fields was $2N + 1 = 21$. The first irradiation field had a convex shape while the other two were non-convex (see second column of figure 4). The second irradiation field had a narrow aperture between the leaves while the first and the third had a much wider one. The orientation angles of the MLC were arbitrarily chosen at 20°, 35° and 15° respectively.

**Simulated data** We simulated the measurements using equation (5) with the penumbra parameter set to $\sigma = 2.1$ mm and considering $P = 6$ projections (equally distributed over 180 degrees) with 128 pixels per projection (with a 0.4 mm pitch). These parameters were taken to mimic our prototype detector (see section 2.3.2). In order to simulate a realistic noise level in the measurements, we added 1 % multiplicative gaussian noise
to the simulated projections. A precise estimation of the noise magnitude in the real measurements would require a specific study that was out of the scope of this study. We therefore selected 1 % for the noise factor to obtain simulated noise magnitudes that qualitatively resembled those observed in the real data (see left column of figure 11).

*Reconstruction from a rectangular field initial estimate* The first simulation experiment aimed at testing the reconstruction algorithm for initial images far from the true solution (i.e. without prior knowledge about the field to reconstruct) to assess its capability and limitations, free from any real delivery errors. The initial estimate $\mu^{(0)}$ was a centered rectangular irradiation field of width 10 mm (see left column of figure 4) and the stopping criterion parameter was set to $\epsilon = 10^{-3}$.

*Detection of treatment delivery errors* We then performed three more simulation experiments mimicking several treatment delivery errors of the linear accelerator that could arise in the clinical context, in order to test the ability of the algorithm to detect them. The magnitudes of the introduced deviations for each of the three types of errors were chosen to be representative of deviations that have been reported in clinics (see for example [23] and [24]). For these three simulation experiments, the initial estimate $\mu^{(0)}$ was set equal to the treatment plan and the stopping criterion parameter was set to $\epsilon = 10^{-3}$.

We first simulated dose deviation errors: the delivered dose was either lower or greater than the planned one but the location and shape of the irradiation field corresponded to the treatment plan. The delivered dose was 5 % greater, 5 % lower, and 3 % greater for the first, second, and third treatment plans, respectively. For this particular simulation, we assume that our system has been previously calibrated against an absolute dosimeter under reference irradiation conditions, and allows absolute dose evaluation.

The second type of error that was simulated was a global spatial shift of the irradiation field: the delivered dose and field shape corresponded to the planned ones, but the whole delivered field was translated with respect to the planned field. The considered translations were (–4 mm, –7 mm) for the first treatment plan, and (–4 mm, +1 mm) and (+2 mm, +2 mm) for the second and third ones, respectively.

We finally considered individual leaf shifts: one or two leaves had an erroneous position but the dose *d* and field location were correct. For the first treatment plan, the right leaf of the sixth leaf-pair was translated by 7 mm to the right while the left leaf of the same leaf-pair was translated by 14 mm to the left for the second treatment plan. For the last treatment plan, the left leaves of the third and sixth leaf-pairs were both translated by 10 mm to the left.

### 2.3.2. Preliminary testing on real data
**System description** Our SciFi prototype detector simultaneously measures $P = 6$ projections (equally distributed over 180 degrees), with 128 pixels per projection and with 0.4 mm pixel pitch [12, 11]. Further description of the detector can be found in [11]. The detector was placed in an RW3 solid water phantom at a source-to-detector distance of 100 cm and at the depth of maximum dose (i.e. 1.4 cm for 6 MV beams). Irradiation was performed by a Novalis TrueBeam STx equipped with a HD120 high definition multileaf collimator. The width of the innermost 32 pairs of leaves was 2.5 mm and 5 mm for the 28 outer ones, as projected at isocenter (i.e. onto the detector plane). In this work, we only used the innermost pairs of leaves and we set $h = 2.5$ mm. Figure 3 shows a schematic of our system with the 6 measured projections.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Schematic of our scintillating-fiber detector. Left: global view ; right: measurement for each ribbon orientation.</figcaption></figure>

*Measured fields* Two reference dose distributions were measured with EBT3 Gafchromic films. The first one was a 10 × 20 mm² rectangular field (see top-left image of figure 8) shaped by $N = 8$ open pairs of leaves and the second one was a typical small IMRT segment (see top-left image of figure 10) shaped by $N = 9$ open pairs of leaves. In addition, we implemented another typical small IMRT segment in the treatment planning system (TPS) (see top-left image of figure 13) which has not been measured with a radiochromic film. The latter was shaped by $N = 12$ pairs of leaves but was much thinner than the one measured with a Gafchromic film. In the rest of this article, the segment shaped with $N = 9$ leaf-pairs is referred as the first one, while the one shaped with $N = 12$ pairs of leaves is referred as the second one.
**Measured data** Using our prototype detector, we performed acquisitions of both typical small IMRT segments (top-left images of figures 10 and 13) with the collimator oriented at 45°. Since the ribbons of scintillating fibers of our detector were stacked on top of each other, the zero-order moments of the projections exhibited small variations that made the projections inconsistent with each other. We compensated for this inconsistency by dividing each projection by the corresponding zero-order moment, and multiplying by the average of their six estimates.

**Determination of the penumbra parameter** Aiming at determining the penumbra parameter $\sigma$, we performed an experiment consisting of comparing the fields constructed with the model of equation (1) (so-called *modeled fields*) with the $10 \times 20 \text{ mm}^2$ reference rectangular field. We considered values of the gaussian standard deviation $\sigma$ varying from 1.9 mm to 2.3 mm, with a 0.1 mm step. The determination of the penumbra parameter should be performed once and for all on a given radiotherapy device.

**Validation of the field model** We validated the field model of equation (1) with a comparison between the reference dose distribution of the first typical IMRT segment and the corresponding field simulated using the proposed model with a penumbra parameter set to $\sigma = 2.1 \text{ mm}$, as determined by our search using the rectangular field. The geometric parameters used to construct the modeled field were directly obtained from the TPS.

**Reconstruction from measured data** The reconstruction algorithm was then tested using the measured data. The initial estimate $\mu^{(0)}$ for the reconstruction procedures was a centered rectangular irradiation field of width 10 mm, the penumbra parameter was set to $\sigma = 2.1 \text{ mm}$ (again, as determined by our search using the rectangular field) and the stopping criterion parameter was set to $\epsilon = 10^{-3}$.

**2.3.3. Reconstruction quality evaluation** The accuracy of the proposed reconstruction method was evaluated using the global gamma index [25], computed as follows:

$$ \gamma(x, y) = \min_{x', y'} \sqrt{\frac{(x - x')^2 + (y - y')^2}{\text{DTA}^2} + \frac{(f(x', y') - \hat{f}(x, y))^2}{\text{DD}^2}}, \quad (14) $$

where $f$ and $\hat{f}$ are respectively the reconstructed and reference irradiation fields, DTA is the distance-to-agreement criterion (in mm) and DD is the dose-difference criterion defined as follows:

$$ \text{DD} = \rho \max_{x,y} \hat{f}(x, y), \quad (15) $$

where $\rho \in [0, 1]$. For example, a gamma index criterion of 2%/1 mm means that $\rho = 0.02$ and DTA = 1 mm. The gamma index test is said to be satisfied for given DD and DTA and for a given point $(x, y)$ if $\gamma(x, y) \le 1$. In the simulation experiments, we
used the widely adopted 2%/2 mm tight criterion [5]. For determining the penumbra
parameter and validating the field model we used, apart from the 2%/2 mm criterion,
the tighter 2%/1 mm criterion for a more precise evaluation. In order to provide
complete information about the performances of the proposed system and reconstruction
method, we used the 3%/3 mm, 2%/2 mm and 2%/1 mm gamma criteria to evaluate
the reconstructions from real data. The mean, maximum and pass rate (fraction of
pixels with a satisfied gamma test) of the gamma index images were all computed by
ignoring all pixels of the reconstructed irradiation field image whose dose value was less
than 10 % of the maximum dose value.

The reference fields for the simulation experiments were the planned and delivered
(planned with introduced error) fields. For the preliminary testing on real data, the
reference fields were obtained from scanned pictures of EBT3 Gafchromic film, except
for the second field for which the reference field was the field image provided by the TPS.
To obtain the reference fields from the fields measurements, we first applied a median
filter to the film images, using a 0.42 × 0.42 mm² square kernel. We then applied a
gaussian filter with a standard deviation of 0.42 mm. The dose response was then
linearized [26] and normalized by the average value of the 3 % pixels having the greatest
dose. A rigid and scale-invariant image registration algorithm [27] was used to align
the reference with the image to be evaluated. Finally, we down-sampled the images to
500 × 500 pixel images with square pixels of width 0.08 mm.

According to IAEA TRS-398 [28], "There are three basic methods currently used for the absolute determination of absorbed dose to water: calorimetry, chemical dosimetry and ionization dosimetry". Neither radiochromic films nor plastic scintillators can be considered for absolute dosimetry. Both technologies can only be used for relative dosimetry: each technology must be cross-calibrated against a reference ionization chamber under reference conditions as explained for instance in [29] and [30] for radiochromic films and plastic scintillators, respectively. In this article, we therefore only performed relative dose evaluation between films and the proposed system.

3. Results

3.1. *Simulated data*

3.1.1. *Initialization with a rectangular field* The results of the first simulation experiment, where the delivered irradiation fields corresponded to the three chosen treatment plans and were reconstructed from an initial rectangular field, are shown in figure 4. The gamma index images evaluate the reconstructed fields with the planned ones taken as reference, and with the 2%/2 mm criterion. The number of iterations for the reconstructions and the corresponding computation times for the same experiment are shown in table 1, along with the 2%/2 mm gamma mean, maximum and pass rates.

3.1.2. *Detection of treatment delivery errors*
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Results of testing the reconstruction algorithm on simulated data (initialized with a rectangular field). Rows one, two and three correspond to the first, second and third treatment plans, respectively. Columns one, two and three correspond to the initial estimates, planned and reconstructed irradiation fields, respectively; column four corresponds to the 2%/2 mm gamma index maps evaluating the reconstructed fields against the planned ones.</figcaption></figure>

<table><caption>**Table 1.** Quantitative results for testing the reconstruction algorithm on simulated data. The gamma analysis refers to the 2%/2 mm gamma criterion. The computation times correspond to a single CPU operating at 3.6GHz.</caption><thead><tr><th>Treatment plan #</th><th>Treatment plan 1</th><th>Treatment plan 2</th><th>Treatment plan 3</th></tr></thead><tbody><tr><td>Number of iterations</td><td>272</td><td>117</td><td>211</td></tr><tr><td>Computation time (s)</td><td>74</td><td>34</td><td>61</td></tr><tr><td>Mean gamma</td><td>0.29</td><td>0.25</td><td>0.16</td></tr><tr><td>Max gamma</td><td>1.26</td><td>0.98</td><td>1.16</td></tr><tr><td>Pass rate (%)</td><td>99.7</td><td>100</td><td>99.8</td></tr></tbody></table>

*Dose deviation* The first type of simulated treatment delivery errors evaluated the ability of the algorithm to measure dose deviations of +5 % for the first treatment plan and −5 % and +3 % for the other two. The results are shown in figure 5 for each of the three treatment plans, where columns four and five show the 2%/2 mm gamma index maps evaluating the reconstructed fields against the planned and the simulated erroneous delivered fields, respectively. The quantitative results for this experiment are shown in table 2.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Results of the first type of simulated treatment delivery errors (dose deviation). Rows one, two and three correspond to the first, second and third treatment plans, respectively. Columns one, two and three correspond to the planned, simulated delivered (planned with error) and reconstructed irradiation fields, respectively; columns four and five correspond to the 2%/2 mm gamma index maps evaluating the reconstructed field against the planned and the simulated delivered fields, respectively.</figcaption></figure>

<table><caption>**Table 2.** Quantitative results of the first type of simulated treatment delivery errors (dose deviation). The gamma analysis refers to the 2%/2 mm gamma criterion. The computation times correspond to a single CPU operating at 3.6GHz.</caption><thead><tr><th>Treatment plan #</th><th>Treatment plan 1</th><th>Treatment plan 2</th><th>Treatment plan 3</th></tr></thead><tbody><tr><td>Number of iterations</td><td>285</td><td>12</td><td>69</td></tr><tr><td>Computation time (s)</td><td>82</td><td>3</td><td>20</td></tr><tr><th>Reference field</th><th>plan.</th><th>del.</th><th>plan.</th><th>del.</th><th>plan.</th><th>del.</th></tr><tr><td>Mean gamma</td><td>0.66</td><td>0.08</td><td>0.53</td><td>0.02</td><td>0.52</td><td>0.06</td></tr><tr><td>Max gamma</td><td>2.52</td><td>0.87</td><td>2.53</td><td>0.03</td><td>1.44</td><td>0.87</td></tr><tr><td>Pass rate (%)</td><td>81.8</td><td>100</td><td>90.0</td><td>100</td><td>88.6</td><td>100</td></tr></tbody></table>
**Global spatial shift** Figure 6 shows the results of the second type of treatment delivery errors, where the simulated errors were global spatial shifts by 8.1 mm, 4.1 mm and 2.8 mm for the first, second and third treatment plans, respectively. The quantitative results for the same experiment are shown in table 3. Again, the 2%/2 mm gamma criterion was considered.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Results of the second type of simulated treatment delivery errors (global spatial shift). Rows one, two and three correspond to the first, second and third treatment plans, respectively. Columns one, two and three correspond to the planned, simulated delivered (planned with error) and reconstructed irradiation fields, respectively; columns four and five correspond to the 2%/2 mm gamma index maps evaluating the reconstructed field against the planned and the simulated delivered fields, respectively.</figcaption></figure>

<table><caption>**Table 3.** Quantitative results of the second type of simulated treatment delivery errors (global spatial shift). The gamma analysis refers to the 2%/2 mm gamma criterion. The computation times correspond to a single CPU operating at 3.6GHz.</caption><thead><tr><th>Treatment plan #</th><th>Treatment plan 1</th><th>Treatment plan 2</th><th>Treatment plan 3</th></tr></thead><tbody><tr><td>Number of iterations</td><td>75</td><td>99</td><td>61</td></tr><tr><td>Computation time (s)</td><td>20</td><td>27</td><td>18</td></tr><tr><td>Reference field</td><td>plan.</td><td>del.</td><td>plan.</td><td>del.</td><td>plan.</td><td>del.</td></tr><tr><td>Mean gamma</td><td>2.07</td><td>0.26</td><td>1.40</td><td>0.32</td><td>0.86</td><td>0.09</td></tr><tr><td>Max gamma</td><td>4.06</td><td>0.86</td><td>2.26</td><td>0.94</td><td>1.66</td><td>1.25</td></tr><tr><td>Pass rate (%)</td><td>19.4</td><td>100</td><td>22.1</td><td>100</td><td>46.5</td><td>99.6</td></tr></tbody></table>
**Individual leaf shifts** The last type of simulated treatment delivery errors dealt with individual leaf shifts. For the first treatment plan, one right-hand side leaf was shifted by 7 mm to the right. For the second one, one left-hand side leaf was shifted by 14 mm to the left. For the third treatment plan, two left-hand side leaves were shifted to the left, both by 10 mm. The results are shown in table 4 and figure 7, where the gamma analysis refers, again, to the 2%/2 mm gamma criterion.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Results of the third type of simulated treatment delivery errors (individual leaf shifts). Rows one, two and three correspond to the first, second and third treatment plans, respectively. Columns one, two and three correspond to the planned, simulated delivered (planned with error) and reconstructed irradiation fields, respectively; columns four and five correspond to the 2%/2 mm gamma index maps evaluating the reconstructed field against the planned and the simulated delivered fields, respectively.</figcaption></figure>

## 3.2. Preliminary testing on real data

### 3.2.1. Determination of the penumbra parameter
The first experiment using real data aimed at determining the penumbra parameter. The quantitative results (mean, maximum and pass rate of the 2%/2 mm and 2%/1 mm gamma index images) are shown in table 5, for $\sigma$ varying from 1.9 mm to 2.3 mm. The results for the 3%/3 mm criterion are not shown since the pass rates for the tighter 2%/2 mm criterion all are 100%.

We show the modeled planned field constructed with a penumbra parameter $\sigma = 2.1$ mm in figure 8, along with the reference field image and the gamma index images. The horizontal and vertical central profiles of the reference and modeled fields are shown in figure 9.
<table><caption>**Table 4.** Quantitative results of the third type of simulated treatment delivery errors (individual leaf shifts). The gamma analysis refers to the 2%/2 mm gamma criterion. The computation times correspond to a single CPU operating at 3.6GHz.</caption><thead><tr><th>Treatment plan #</th><th>Treatment plan 1</th><th>Treatment plan 2</th><th>Treatment plan 3</th></tr></thead><tbody><tr><td>Number of iterations</td><td>594</td><td>375</td><td>338</td></tr><tr><td>Computation time (s)</td><td>168</td><td>104</td><td>99</td></tr><tr><th>Reference field</th><th>plan.</th><th>del.</th><th>plan.</th><th>del.</th><th>plan.</th><th>del.</th></tr><tr><td>Mean gamma</td><td>0.30</td><td>0.20</td><td>0.52</td><td>0.13</td><td>0.52</td><td>0.08</td></tr><tr><td>Max gamma</td><td>2.61</td><td>1.13</td><td>2.89</td><td>0.80</td><td>4.22</td><td>0.86</td></tr><tr><td>Pass rate (%)</td><td>94.0</td><td>99.8</td><td>81.2</td><td>100</td><td>83.5</td><td>100</td></tr></tbody></table>

<table><caption>**Table 5.** Quantitative results for determining the penumbra parameter.</caption><thead><tr><th rowspan="2">Gamma index criterion</th><th colspan="4">2%/2 mm</th><th colspan="4">2%/1 mm</th></tr><tr><th>1.9</th><th>2.0</th><th>2.1</th><th>2.2</th><th>2.3</th><th>1.9</th><th>2.0</th><th>2.1</th><th>2.2</th><th>2.3</th></tr></thead><tbody><tr><td>σ (mm)</td><td>1.9</td><td>2.0</td><td>2.1</td><td>2.2</td><td>2.3</td><td>1.9</td><td>2.0</td><td>2.1</td><td>2.2</td><td>2.3</td></tr><tr><td>Mean gamma</td><td>0.20</td><td>0.21</td><td>0.23</td><td>0.25</td><td>0.27</td><td>0.33</td><td>0.35</td><td>0.39</td><td>0.43</td><td>0.48</td></tr><tr><td>Max gamma</td><td>0.74</td><td>0.73</td><td>0.73</td><td>0.73</td><td>0.73</td><td>1.35</td><td>1.18</td><td>1.02</td><td>1.01</td><td>1.11</td></tr><tr><td>Pass rate (%)</td><td>100</td><td>100</td><td>100</td><td>100</td><td>100</td><td>99.5</td><td>99.8</td><td>100.</td><td>100.</td><td>98.8</td></tr></tbody></table>

### 3.2.2. Validation of the field model
The objective of the second experiment was to validate the irradiation field model of equation (1) with the penumbra parameter set to $\sigma = 2.1$ mm by performing gamma analysis of the modeled first field, taking the EBT3 film as reference field. The quantitative results for the 2%/2 mm and 2%/1 mm gamma criteria are shown in table 6. For the same experiment, the reference and modeled fields, along with the gamma index images are shown in figure 10.

<table><caption>**Table 6.** Quantitative results for validating the irradiation field model.</caption><thead><tr><th>Gamma index criterion</th><th>2%/2 mm</th><th>2%/1 mm</th></tr></thead><tbody><tr><td>Mean gamma</td><td>0.25</td><td>0.44</td></tr><tr><td>Max gamma</td><td>1.07</td><td>1.24</td></tr><tr><td>Pass rate (%)</td><td>99.9</td><td>98.6</td></tr></tbody></table>

### 3.2.3. Reconstructions from measured data
For the experimental cases, the algorithm stopped after 248 iterations for the first field, and 109 iterations for the second one, corresponding to computation times of 82 s and 33 s on a single CPU operating at 3.6 GHz, respectively. The quantitative results of these physical experiments are shown in table 7 (first field) and table 8 (second field).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Results of determining the penumbra parameter. Top-left: reference field (EBT3 film) ; top-right: modeled planned field with $\sigma = 2.1\text{mm}$ ; bottom: gamma index maps.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Results of determining the penumbra parameter. Profiles of the reference (blue) and modeled (orange) rectangular fields. Left: central horizontal profiles ($y = 0$) ; right: central vertical profiles ($x = 0$).</figcaption></figure>

The measured sinograms and the reconstructed ones (the sinograms obtained from the reconstructed fields) of these experiments are shown in figure 11 along with their absolute differences. Figures 12 and 13 show the irradiation fields and gamma index image for the first and second fields, respectively.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Results of validating the irradiation field model. Top-left: reference field (EBT3 film); top-right: modeled planned field with $\sigma = 2.1\text{mm}$; bottom: gamma index maps evaluating the modeled field against the reference.</figcaption></figure>

<table><caption>**Table 7.** Quantitative results for testing the reconstruction algorithm on real data (first field): gamma analysis evaluating the reconstructed field against the reference and planned fields.</caption><thead><tr><th></th><th>Gamma index criterion</th><th>ref.</th><th>plan.</th></tr></thead><tbody><tr><td>Mean gamma</td><td rowspan="3">3%/3 mm</td><td>0.18</td><td>0.09</td></tr><tr><td>Max gamma</td><td>0.72</td><td>0.33</td></tr><tr><td>Pass rate (%)</td><td>100</td><td>100</td></tr><tr><td>Mean gamma</td><td rowspan="3">2%/2 mm</td><td>0.27</td><td>0.14</td></tr><tr><td>Max gamma</td><td>1.08</td><td>0.49</td></tr><tr><td>Pass rate (%)</td><td>99.9</td><td>100</td></tr><tr><td>Mean gamma</td><td rowspan="3">3%/3 mm</td><td>0.48</td><td>0.24</td></tr><tr><td>Max gamma</td><td>1.53</td><td>0.96</td></tr><tr><td>Pass rate (%)</td><td>95.7</td><td>100</td></tr></tbody></table>

## 4. Discussion

This paper presents a reconstruction algorithm for quality assurance of IMRT segments, using a detector made of scintillating-fiber ribbons with 6 projection angles. For the purpose of this proof-of-concept, the algorithm has been tested on both simulated and
<table><caption>**Table 8.** Quantitative results for testing the reconstruction algorithm on real data (second field): gamma analysis evaluating the reconstructed field against the planned field provided by the TPS.</caption><thead><tr><th>Gamma index criterion</th><th>3%/3 mm</th><th>2%/2 mm</th><th>2%/1 mm</th></tr></thead><tbody><tr><td>Mean gamma</td><td>0.19</td><td>0.28</td><td>0.49</td></tr><tr><td>Max gamma</td><td>0.86</td><td>1.28</td><td>2.34</td></tr><tr><td>Pass rate (%)</td><td>100</td><td>99.4</td><td>92.6</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Physical experiments sinograms. Top: first field ; bottom: second field ; left: measurements ; middle: sinograms obtained from the reconstructed fields ; right: absolute differences.</figcaption></figure>

experimental data and has shown its capability to efficiently detect some treatment delivery errors that could arise in the clinical context.

### 4.1. Irradiation fields modeling

A key challenge for the reconstruction of IMRT segments is to accurately model the penumbra. In the experiment for determining the penumbra parameter, rectangular fields constructed with the proposed model and with different values of the standard deviation of the gaussian were compared to a reference EBT3 film. For the 2%/1 mm gamma criterion, the pass rate reached a maximum for $\sigma = 2.1$ mm and $\sigma = 2.2$ mm (see table 5). For the same criterion, the maximum gamma value for these penumbra parameters was very close (1.02 and 1.01, respectively), and the mean gamma value was lower for $\sigma = 2.1$ mm (0.39) than for $\sigma = 2.2$ mm (0.43) so we selected the value
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Results of testing the reconstruction algorithm on real data (initialized with a rectangular field) for the first field. Top-left: reference field (EBT3 film) ; top-middle: planned field (from equation (1)) ; top-right: reconstructed field ; row two: gamma index maps evaluating the reconstructed field against the reference ; row three: gamma index maps evaluating the reconstructed field against the planned one.</figcaption></figure>

$\sigma$ = 2.1 mm for the rest of the article. The MLC and LINAC manufacturer specified a 20-80 % penumbra width lower than 3.5 mm [31], which was consistent with the gaussian model with $\sigma$ = 2.1 mm, for which the 20-80 % penumbra width was 2.4 mm. The choice of $\sigma$ = 2.1 mm was also consistent with the penumbra measured at isocenter and at the depth of 1.5 cm with the same MLC but with an older LINAC model (Novalis TX), which was of 2.5 mm (as defined between 20-80 %) for a 20-by-20 mm$^2$ irradiation field [32].

For the IMRT segment, the penumbra region was large for such a small field with respect to the maximum dose region, making the gamma test difficult to fully satisfy. For example, most of the unsatisfied 2%/1 mm gamma pixels were located in the penumbra region (see figures 10 and 12). Other strategies for the penumbra modeling (e.g. using
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Results of testing the reconstruction algorithm on real data (initialized with a rectangular field) for the second field. Top-left: planned field (from TPS); top-right: reconstructed field; bottom: gamma index maps evaluating the reconstructed field against the planned one.</figcaption></figure>

an anisotropic gaussian function instead of an isotropic one in equation (2) or using a combination of several gaussian functions instead of a single one) might improve the accuracy of the proposed algorithm. This would require further investigations and was out of the scope of this work. Furthermore, a method for experimental determination of the penumbra parameters (a couple of parameters instead of a single one if an anisotropic gaussian blurring or a combination of two isotropic gaussian blurrings is used instead of our single isotropic gaussian blurring) and a more comprehensive experimental validation of the proposed field model would be required to bring our system to the clinic. The penumbra parameters determination could be performed, for example, by a search similar to the one presented in this article, but with a variety of field shapes and sizes.

## 4.2. Reconstruction algorithm

The reconstruction of the irradiation fields starting from a rectangular irradiation field as initial estimate (i.e. without prior knowledge about the leaf positions) has shown satisfactory results for the three simulated treatment plans with 2%/2 mm gamma pass rates greater than 99.7% (see figure 4 and table 1).

From the simulations of treatment delivery errors, the proposed reconstruction algorithm has shown its potential to perform quality assurance of IMRT segments in a clinical context, based on the tight 2%/2 mm gamma index criterion.
The results of the first simulated errors, where only the delivered dose was different from the treatment plan, indicate that these situations were very well managed. For each of the three treatment plans, the 2%/2 mm gamma pass rates achieved a maximum of 90 % when comparing the reconstructed field with the planned field (even for a dose deviation as low as 3 %), and equal to 100 % when comparing it with the delivered field (see figure 5 and table 2). As the dose estimation was implicitly found from the zero-order moment $M^0$ and depended only on the field geometry (see equation (9)), the proposed approach was very well suited to these dose deviation situations. Without noise in the measurements, the projection of the initial estimate (i.e. the planned field) would fit perfectly and the cost function value and its gradient would be zero, resulting in no iterations.

The second type of treatment delivery errors that were simulated demonstrated the ability of the algorithm to detect a translation of the whole irradiation field compared to the treatment plan, even by less than 3 mm (see figure 6) with, again, low gamma pass rates when the reconstructed irradiation field was evaluated against the planned one (lower than 50 %) and satisfactory pass rates when the reconstructed field was evaluated against the delivered field (greater than 99.6 %).

The algorithm was also able to identify errors in individual leaf positions. In figure 7, we note that the delivered irradiation fields have been accurately reconstructed. The simulated error for the first treatment plan was a shift by 7 mm of an individual leaf and the gamma pass rate evaluating the reconstructed irradiation field against the planned one was lower than 95 %. Again, the gamma pass rate when the delivered field was taken as reference was satisfactory (99.8 %). It is worth mentioning that a 95 % threshold is widely adopted for treatment QA in clinical routine [5]. For the second and third treatment plans, the detection of the leaf shifts was easier (i.e. the gamma pass rates evaluating the reconstructed fields against the planned fields were much lower) since we simulated a greater shift for the second treatment plan, and the shift of two leaves for the third one.

The proposed reconstruction method has also shown satisfactory results on real data. The difference between the measurements and the projections of the reconstructed irradiation fields remained lower than 5 % of the greatest projection values for both fields (see figure 11). In addition, starting from an initial estimate with different doses, locations and shapes, the reconstructed fields matched very well the reference irradiation fields acquired with an EBT3 Gafchromic film for the first field (see figure 12) and the TPS field image for the second field (see figure 13). For the first field, all the gamma index tests exhibited pass rates greater than 95 %, resulting in the same QA conclusion when evaluating the reconstructed field against the reference or the planned one. Observing figure 12, we note that, similar to the testing of the field model, most of the unsatisfied pixels of the 2%/1 mm gamma index map comparing the reconstructed field with the reference lie in the penumbra region. For the second field, the gamma pass rates for the 3%/3 mm and the 2%/2 mm gamma criteria were greater than 95 %, but was only 92.6 % for the 2%/1 mm criterion. It is interesting to note that the low
gamma pass rate for the 2%/1 mm criterion comes from the substantial role played
by the penumbra for such a thin field. This suggests that a more accurate penumbra
modeling could yield submillimetric quality assurance for small MLC fields.

The observed computation times were between 6 s and 3 min as measured with
our Python implementation on a single CPU. A low-level programming language
implementation (e.g. C/C++) coupled with the use of GPUs can lower by several
orders of magnitude the reconstruction time. The computation time also depended
on the number of iterations to reach convergence and the computation time for a
single iteration. The latter depends on the number of leaves and the number of
iterations in the line-search algorithm. An optimization of the trade-off between the
number of iterations and the number of line-search iterations was out of the scope of
this study, but could probably substantially lower the computation time of the whole
algorithm. It is worth mentioning that IMRT plans typically have ten beams, each
containing a hundred segments but the reconstruction of these segments can be fully
parallelized. The proposed reconstruction algorithm is therefore appropriate for real-
time QA measurements.

A fundamental limitation of the use of scintillating-fiber detectors for reconstructing MLC irradiation fields is that for a large number of pairs of leaves *N*, there can be many non-uniqueness configurations (i.e. multiple global minima). For 6 projections, non-uniqueness has been demonstrated for as few as *N* = 32 leaf-pairs [21]. Increasing the number of measured projections would be required for handling more pairs of leaves.

The optimized cost function is non-convex so the convergence to the global
minimum cannot be ensured using the steepest descent algorithm. To ensure reliable
tomographic reconstructions of irradiation fields, modification of the optimization
algorithm to escape any potential local minimum would be required by using, for
example, a basin-hopping approach [33] that alternates local minimization (as presented
in this article) and random perturbations. This was beyond the scope of this study.

The reconstruction algorithm uses the analytical projections of the irradiation field
model of equation (1). The latter is defined under the assumption that the gantry angle
is 0 which is not always the case in practice. In order to deal with non-zero gantry
angle situations, a geometric linear transformation should be applied to the field model
of equation (1) and the analytical formulas for the measured projections (equation (4))
and the cost function gradient (Appendix) should be modified accordingly.

4.3. Toward clinical applications

This study illustrated the potential of tomographic reconstruction of IMRT segments
using scintillating-fiber ribbons to perform high-resolution 2D beam characterization.
We have shown that this technique can verify the accuracy of the MLC position and
detect individual leaf position errors. Furthermore, the success rate determined on
real data (with EBT3 film measurements used as a reference) was greater than 95 %,
which is consistent with the AAPM TG-218 recommendations for quality assurance of
verification based on IMRT measurements, but using even tighter criteria of 2%/1 mm
instead of the recommended 3%/2 mm [5]. If these preliminary results are confirmed
by a more comprehensive study, the proposed technique could be considered for IMRT
QA measurements using the Perpendicular Field by Field (PFF) delivery method [5].
This is an appropriate QA method for fixed-gantry IMRT delivery employing step-
and-shoot method. For rotational IMRT delivery such as tomotherapy and Volumetric
Modulated Arc Therapy (VMAT), according to AAPM TG-218 recommendations, QA
measurements should be performed using the True Composite (TC) delivery method.
TC delivery methods require a QA device which has negligible angular dependence.
This is the case for our technology based on plastic scintillating-fiber ribbons which
have been shown to be tissue-equivalent detectors with low angular dependence [34, 35].
However, the technique proposed for tomographic reconstruction of IMRT segments
would need to be further developed since it currently assumes that the radiation beam
is perpendicular to the measurement device detector plane, which would no longer be
the case for TC delivery method.

**5. Conclusion**

In this work, we have proposed an iterative reconstruction method for quality assurance suitable for radiotherapy irradiation fields shaped with multileaf collimators, which uses a detector made of scintillating fibers with six projections. The method has been tested on real data acquired with our prototype detector, and has shown satisfactory gamma pass rates down to the 2%/2 mm criterion. QA success with the 2%/1 mm criterion remains challenging for fields with a large penumbra region compared to the field size. In addition, based on the simulation of treatment delivery errors that could arise in the clinical context, the algorithm has shown its ability to identify a deviation in the delivered dose and in the irradiation field location or shape. These promising results pave the way towards the use of scintillating-fiber detectors for real-time quality assurance of IMRT segments in clinical routine.
# Appendix

## General expression of the cost function gradient

The elements $\frac{\partial C(\mu)}{\partial \mu_i}$ of the cost function gradient $\nabla C(\mu) \in \mathbb{R}^{2N+1}$ are computed as follows:

$$ \frac{\partial C(\mu)}{\partial \mu_i} = \sum_{\phi \in \Phi} \int_{-\infty}^{\infty} 2 \frac{\partial p_{\mu}(\phi, s)}{\partial \mu_i} (p_{\mu}(\phi, s) - q(\phi, s)) ds, \quad (16) $$

where $\frac{\partial p_{\mu}(\phi, s)}{\partial \mu_i}$ is the forward model partial derivative according to the i-th variable given by the following expression:

$$ \frac{\partial p_{\mu}}{\partial \mu_i} = \frac{\partial}{\partial \mu_i} \left( \frac{M^0}{h \sum_{l=1}^{N} w_l} \right) \times \sum_{l=1}^{N} R_{\phi} \chi_{\mu,l} *_{1D} R_{\phi} g + \frac{M^0}{h \sum_{l=1}^{N} w_l} \times \sum_{l=1}^{N} \frac{\partial}{\partial \mu_i} (R_{\phi} \chi_{\mu,l} *_{1D} R_{\phi} g). \quad (17) $$

We recall that $\mu = (w_1, w_2, \dots, w_N, x_1, x_2, \dots, x_N, y_1)^T$ (see section 2.1.1).

## Partial derivatives of the dose

The dose partial derivatives are the following:

$$ \frac{\partial}{\partial w_i} \left( \frac{M^0}{h \sum_{l=1}^{N} w_l} \right) = - \frac{M^0}{h \left( \sum_{l=1}^{N} w_l \right)^2} \quad (18) $$

$$ \frac{\partial}{\partial x_i} \left( \frac{M^0}{h \sum_{l=1}^{N} w_l} \right) = 0 \quad ; \quad \frac{\partial}{\partial y_i} \left( \frac{M^0}{h \sum_{l=1}^{N} w_l} \right) = 0. \quad (19) $$

## Partial derivatives of the projection of the unit intensity rectangle

The partial derivatives of the unit intensity rectangular irradiation field projection are given by:

$$ \frac{\partial}{\partial \mu_i} (R_{\phi} \chi_{\mu,l} *_{1D} R_{\phi} g) = \begin{cases} B, & \text{if } \cos \phi \sin \phi = 0 \\ A + C, & \text{if } w_l |\sin \phi| = h |\cos \phi|, \\ A + B + C & \text{otherwise} \end{cases} \quad (20) $$

where the conditions $\cos \phi \sin \phi = 0$ and $w_l |\sin \phi| = h |\cos \phi|$ are fulfilled when $R_{\phi} \chi_{\mu,l}$ is a rectangular and a triangular function, respectively. The expressions of A, B and C
are obtained using the Leibniz integral rule:

$$
\begin{align}
A ={}& \frac{\partial m_l}{\partial \mu_i} \times \left( \frac{1}{\delta_l^+ - \delta_l^-} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} t R_{\phi} g(s-t) dt + \frac{\delta_l^+ - s_l}{\delta_l^+ - \delta_l^-} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} R_{\phi} g(s-t) dt \right) \tag{21} \\
& + \frac{\partial s_l}{\partial \mu_i} \times \left( m_l \times R_{\phi} g(s - (s_l - \delta_l^-)) - \frac{m_l}{\delta_l^+ - \delta_l^-} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} R_{\phi} g(s-t) dt \right) \nonumber \\
& + \frac{\partial \delta_l^-}{\partial \mu_i} \times \left( -m_l \times R_{\phi} g(s - (s_l - \delta_l^-)) + \frac{m_l}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} t R_{\phi} g(s-t) dt \right) \nonumber \\
& + m_l \frac{\delta_l^+ - s_l}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} R_{\phi} g(s-t) dt \nonumber \\
& + \frac{\partial \delta_l^+}{\partial \mu_i} \times \left( \frac{-m_l}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} t R_{\phi} g(s-t) dt + \frac{m_l(s_l - \delta_l^-)}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l - \delta_l^+}^{s_l - \delta_l^-} R_{\phi} g(s-t) dt \right), \nonumber
\end{align}
$$

$$
\begin{align}
B ={}& \frac{\partial m_l}{\partial \mu_i} \times \int_{s_l - \delta_l^-}^{s_l + \delta_l^-} R_{\phi} g(s-t) dt \tag{22} \\
& + \frac{\partial s_l}{\partial \mu_i} \times m_l \times \left( R_{\phi} g(s - (s_l + \delta_l^-)) - R_{\phi} g(s - (s_l - \delta_l^+)) \right) \nonumber \\
& + \frac{\partial \delta_l^-}{\partial \mu_i} \times m_l \times \left( R_{\phi} g(s - (s_l + \delta_l^-)) + R_{\phi} g(s - (s_l - \delta_l^-)) \right), \nonumber
\end{align}
$$

$$
\begin{align}
C ={}& \frac{\partial m_l}{\partial \mu_i} \times \left( - \frac{1}{\delta_l^+ - \delta_l^-} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} t R_{\phi} g(s-t) dt + \frac{\delta_l^+ + s_l}{\delta_l^+ - \delta_l^-} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} R_{\phi} g(s-t) dt \right) \tag{23} \\
& + \frac{\partial s_l}{\partial \mu_i} \times \left( - m_l \times R_{\phi} g(s - (s_l + \delta_l^-)) + \frac{m_l}{\delta_l^+ - \delta_l^-} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} R_{\phi} g(s-t) dt \right) \nonumber \\
& + \frac{\partial \delta_l^-}{\partial \mu_i} \times \left( - m_l \times R_{\phi} g(s - (s_l + \delta_l^-)) - \frac{m_l}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} t R_{\phi} g(s-t) dt \right) \nonumber \\
& + m_l \frac{\delta_l^+ + s_l}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} R_{\phi} g(s-t) dt \nonumber \\
& + \frac{\partial \delta_l^+}{\partial \mu_i} \times \left( \frac{m_l}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} t R_{\phi} g(s-t) dt - \frac{m_l(\delta_l^- + s_l)}{(\delta_l^+ - \delta_l^-)^2} \int_{s_l + \delta_l^-}^{s_l + \delta_l^+} R_{\phi} g(s-t) dt \right). \nonumber
\end{align}
$$

Note that unlike $R_{\phi}\chi_{\mu,l}$, its gaussian convolution $R_{\phi}\chi_{\mu,l} *_{1D} R_{\phi}g$ is differentiable with respect to any $\mu_i$ for any $s \in \mathbb{R}$.

Integrals involving the gaussian function Equations (21), (22) and (23) contain two integrals involving the gaussian function. They can be computed as follows:

$$
\int_{t_1}^{t_2} R_{\phi} g(s-t) dt = \left[ -\frac{1}{2} \operatorname{erf} \left( \frac{s-t}{\sigma\sqrt{2}} \right) \right]_{t_1}^{t_2}, \qquad (24)
$$

$$
\int_{t_1}^{t_2} t R_{\phi} g(s-t) dt = \left[ \frac{s}{2} \operatorname{erf} \left( \frac{t-s}{\sigma\sqrt{2}} \right) - \sigma^2 R_{\phi} g(s-t) \right]_{t_1}^{t_2}, \quad (25)
$$
where erf is the error function defined as:

$$ \text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt \qquad (26) $$

Partial derivatives of the trapezoid parameters

$$ \frac{\partial m_l}{\partial w_l} = \begin{cases} \frac{1}{|\cos \phi|}, & \text{if } h |\cos \phi| > w_l |\sin \phi| \\ 0, & \text{if } w_l |\sin \phi| > h |\cos \phi| \end{cases} \qquad (27) $$

$$ \frac{\partial \delta_l^{-}}{\partial w_l} = \begin{cases} \frac{|\sin \phi|}{2}, & \text{if } w_l |\sin \phi| > h |\cos \phi| \\ -\frac{|\sin \phi|}{2}, & \text{if } h |\cos \phi| > w_l |\sin \phi| \end{cases} \qquad (28) $$

$$ \frac{\partial s_l}{\partial x_l} = -\sin \phi \quad ; \quad \frac{\partial s_l}{\partial y_l} = \cos \phi \quad ; \quad \frac{\partial \delta_l^{+}}{\partial w_l} = \frac{|\sin \phi|}{2} \qquad (29) $$

$$ \frac{\partial m_l}{\partial x_l} = \frac{\partial m_l}{\partial y_l} = \frac{\partial s_l}{\partial w_l} = \frac{\partial \delta_l^{-}}{\partial x_l} = \frac{\partial \delta_l^{-}}{\partial y_l} = \frac{\partial \delta_l^{+}}{\partial x_l} = \frac{\partial \delta_l^{+}}{\partial y_l} = 0 \qquad (30) $$

References

[1] S. V. Spirou and C. S. Chui. Generation of arbitrary intensity profiles by dynamic jaws or multileaf collimators. *Medical Physics*, 21(7):1031–1041, 1994.

[2] J. Stein, T. Bortfeld, B. Dörschel, and W. Schlegel. Dynamic x-ray compensation for conformal radiotherapy by means of multi-leaf collimation. *Radiotherapy and Oncology*, 32(2):163–173, 1994.

[3] A. S. Shiu, H. M. Kooy, J. R. Ewton, S. S. Tung, J. Wong, K. Antes, and M. H. Maor. Comparison of miniature multileaf collimation (MMLC) with circular collimation for stereotactic treatment. *International Journal of Radiation Oncology-Biology-Physics*, 37(3):679–688, 1997.

[4] S. Devic, N. Tomic, and D. Lewis. Reference radiochromic film dosimetry: Review of technical aspects. *Physica Medica*, 32(4):541–556, 2016.

[5] M. Miften, A. Olch, D. Mihailidis, J. Moran, T. Pawlicki, A. Molineu, H. Li, K. Wijesooriya, J. Shi, P. Xia, et al. Tolerance limits and methodologies for IMRT measurement-based verification QA: recommendations of AAPM Task Group No. 218. *Medical Physics*, 45(4):e53–e83, 2018.

[6] A. Agnew, C. E. Agnew, M. W. D. Grattan, A. R. Hounsell, and C. K. McGarry. Monitoring daily MLC positional errors using trajectory log files and EPID measurements for IMRT and VMAT deliveries. *Physics in Medicine and Biology*, 59(9):N49–N63, 2014.

[7] B. J. Zwan, B. W. King, D. J. O’Connor, and P. B. Greer. Dose-to-water conversion for the backscatter-shielded EPID: A frame-based method to correct for EPID energy response to MLC transmitted radiation. *Medical Physics*, 41(8Part1):081716, 2014.

[8] M. S. Rose, L. Tirpak, K. Van Casteren, J. Zack, T. Simon, A. Schoenfeld, and W. Simon. Multi-institution validation of a new high spatial resolution diode array for srs and sbrt plan pretreatment quality assurance. *Medical Physics*, 47(7):3153–3164, 2020.

[9] M. Goulet, L. Archambault, L. Beaulieu, and L. Gingras. High resolution 2D dose measurement device based on a few long scintillating fibers and tomographic reconstruction). *Medical Physics*, 39(8):4840–4849, 2012.
[10] P. Pittet, J. Esteves, J.-M. Galvan, G.-N. Lu, F. Blanc, G. Haefeli, P. Hopchev, S. Rit, L. Desbat, J. Ribouton, and P. Jalade. SciFi detector and associated method for real-time determination of profile and output factor for small fields in stereotactic radiotherapy. *Medical Physics*, 47(4):1930–1939, 2020.

[11] J. Esteves, O. Pivot, J. Ribouton, P. Jalade, A. Zouaoui, L. Desbat, S. Rit, F. Blanc, G. Haefeli, P. Hopchev, J. M. Galvan, G. N. Lu, and P. Pittet. A novel QA phantom based on scintillating fiber ribbons with implementation of 2D dose tomography for small-field radiotherapy. *Medical Physics*, 50(1):619–632, 2023.

[12] P. Pittet, J. Ribouton, P. Jalade, J. Esteves, F. Blanc, G. Haefeli, P. Hopchev, J. M. Galvan, and G. N. Lu. Highly spatially resolved 2D plastic scintillation detector system for small field real-time QA. In *ESTRO 2020 Congress*, 2020.

[13] L. Desbat, S. Rit, R. Clackdoyle, P. Jalade, J. Ribouton, and P. Pittet. Geometric tomography for measuring rectangular radiotherapy fields from six projections. In *2019 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)*, pages 1–4, 2019.

[14] F. Natterer. *The Mathematics of Computerized Tomography*. Wiley, 1986.

[15] M. Li, H. Yang, and H. Kudo. An accurate iterative reconstruction algorithm for sparse objects: application to 3d blood vessel reconstruction from a limited number of projections. *Physics in Medicine & Biology*, 47(15):2599–2602, 2012.

[16] R. Gordon, R. Bender, and G. T. Herman. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography. *Journal of Theoretical Biology*, 29(3):471–481, 1970.

[17] E. Y. Sidky, C. M. Kao, and X. Pan. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT. *Journal of X-ray Science and Technology*, 14(2):119–139, 2006.

[18] H. Yu and G. Wang. A soft-threshold filtering approach for reconstruction from a limited number of projections. *Physics in Medicine & Biology*, 55(13):3905, 2010.

[19] K. J. Batenburg and J. Sijbers. DART: A practical reconstruction algorithm for discrete tomography. *IEEE Transactions on Image Processing*, 20(9):2542–2553, 2011.

[20] R. J. Gardner. *Geometric tomography*, volume 58. Cambridge University Press Cambridge, 1995.

[21] O. Pivot, R. Clackdoyle, S. Rit, and L. Desbat. Estimation of radiotherapy dose fields from a few projections: How many projections will ensure uniqueness? In *2020 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)*, pages 1–4, 2020.

[22] S. Wright and J. Nocedal. Numerical optimization. Springer Science, 35(67-68):7, 1999.

[23] A. Mans, M. Wendling, L. N. McDermott, J.-J. Sonke, R. Tielenburg, R. Vijlbrief, B. Mijnheer, M. Van Herk, and J. C. Stroom. Catching errors with in vivo epid dosimetry. *Medical physics*, 37(6Part2):2638–2644, 2010.

[24] M. Pasler, J. Kaas, T. Perik, J. Geuze, R. Dreindl, T. Künzler, F. Wittkamper, and D. Georg. Linking log files with dosimetric accuracy–A multi-institutional study on quality assurance of volumetric modulated arc therapy. *Radiotherapy and Oncology*, 117(3):407–411, 2015.

[25] D. A. Low, W. B. Harms, S. Mutic, and J. A. Purdy. A technique for the quantitative evaluation of dose distributions. *Medical Physics*, 25(5):656–661, 1998.

[26] S. Devic, N. Tomic, S. Aldelaijan, F. DeBlois, J. Seuntjens, M. F. Chan, and D. Lewis. Linearization of dose-response curve of the radiochromic film dosimetry system. *Medical Physics*, 39(8):4850–4857, 2012.

[27] M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup. Efficient subpixel image registration algorithms. *Optics letters*, 33(2):156–158, 2008.

[28] *Absorbed Dose Determination in External Beam Radiotherapy*. Number 398 in Technical Reports Series. INTERNATIONAL ATOMIC ENERGY AGENCY, Vienna, 2001.

[29] T. Santos, T. Ventura, and M. do Carmo Lopes. A review on radiochromic film dosimetry for dose verification in high energy photon beams. *Radiation Physics and Chemistry*, 179:109217, 2021.
[30] P. Carrasco, N. Jornet, O. Jordi, M. Lizondo, A. Latorre-Musoll, T. Eudaldo, A. Ruiz, and M. Ribas. Characterization of the exradin w1 scintillator for use in radiotherapy. *Medical physics*, 42(1):297-304, 2015.

[31] VARIAN medical system. TrueBeam STx System Specifications, 2015.

[32] Z. Wang, Z. Chang, Q. Wu, S. Zhou, C. Huntzinger, and F. Yin. SU-GG-T-446: Dosimetric characteristics of high definition multi-leaf collimator. *Medical Physics*, 35(6Part16):2827-2827, 2008.

[33] D. J. Wales and J. P. K. Doye. Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. *The Journal of Physical Chemistry A*, 101(28):5111-5116, 1997.

[34] C. Devic, J. Plagnard, and M. Munier. Characterization of an innovative detector based on scintillating fiber for personalized computed tomography dosimetry. *Sensors*, 22(1):90, 2021.

[35] L. L. W. Wang, D. Klein, and A. S. Beddar. Monte carlo study of the energy and angular dependence of the response of plastic scintillation detectors in photon beams. *Medical Physics*, 37(10):5279-5286, 2010.

## Acknowledgments

With financial support from ITMO Cancer AVIESAN (Alliance Nationale pour les Sciences de la Vie et de la Santé, National Alliance for Life Sciences & Health) within the framework of the Cancer Plan (QAsys project).

This work is supported by the University Grenoble Alpes IDEX grant CQFD, the "Fonds unique interministériel" and the European Union FEDER in Auvergne Rhône Alpes (3D4Carm project) and by the ANR (ROIdoré ANR-17-CE19-0006-01 ; labex CAMI, ANR-11-LABX-0004_01 ; labex PRIMES, ANR-11-LABX-0063 and SPECT-Motion-eDCC, ANR-21-CE45-0026).

The authors would like to thank Josué Esteves, Patrice Jalade and Julien Ribouton for having performed the physical experiments.