Direct reconstruction of the source intensity distribution of a clinical linear accelerator using a maximum likelihood expectation maximization algorithm

Advanced radiotherapy techniques, such as stereotactic radiosurgery (SRS) and intensity modulated radiation therapy (IMRT), commonly utilize small radiation fields (<$2\times 2$ cm²) in order to improve target coverage and reduce the dose to surrounding organs at risk (OAR). However, accurate commissioning of beam models for small fields appears still to be challenging. Most detectors often present significant perturbations in small fields due to volume averaging or the presence of high density materials. During beam model commissioning of linear accelerators, these perturbations may result in mis-configuration of crucial beam parameters in the treatment planning system (TPS). Systematic errors introduced at this step will eventually propagate as dosimetric errors at the patient dose calculation stage.

In order to address this issue, one approach is to apply detector-specific correction factors directly on the measurements (Alfonso et al 2008). Several researchers have shown that corrections for output factors can be derived by Monte Carlo (MC) methods (Cranmer-Sargison et al 2011, Francescon et al 2012). Similar corrections may be needed for dose profiles as well, especially in the tail region (Papaconstadopoulos et al 2014a, Francescon et al 2014). The MC-derived correction factors have been validated by experimental methods (Pantelis et al 2012). However, deriving such corrections is a lengthy process and requires MC expertise as well as detailed knowledge of the design and materials of the detector. Furthermore, each time a new detector is produced, new correction factors are needed. It is clear that alternative pathways for determining the appropriate beam parameters would offer a great service in the dosimetry of small fields.

One of the most important beam parameters is the source size and shape of clinical linear accelerators. In MC beam models, the term source most often refers to the electron spatial distribution incident on the target. In other beam models, it may refer to the bremsstrahlung radiation that is produced in the target, often referred to as the primary x-ray source or focal-spot. In both cases, the source size is most often characterized by the full width at half maximum (FWHM) and the shape is assumed to be a 2D elliptical Gaussian distribution. Sterpin et al (2011) showed that for a nominal 5.5 and 18.0 MeV incident electron energy, the primary x-ray source is located approximately at a depth of 0.13 and 0.25 mm respectively from the top of the target. It was also shown in that study that the photon source FWHM agrees within 0.06 mm to the electron source FWHM for typical source sizes of 0.5–1.5 mm. This important conclusion implies that if direct reconstruction techniques could be developed for the x-ray source, then the electron source could be reproduced as well within the previously stated level of accuracy. Other beam parameters, such as the energy of the incident electron source, have a limited effect in small fields. Furthermore, methods for determining the incident electron energy or even unfolding the full energy spectrum of linear accelerators have been suggested (Ali et al 2012).

As the field becomes small, the primary x-ray source, as seen from the point of measurement, starts to be partially occluded (Aspradakis and IPEM 2010). Photons originating from other scatter sources, such as the flattening filter and primary collimator (extra-focal spot), are almost completely blocked and their dosimetric impact is minimized. The magnitude of the source occlusion effect not only depends on the spatial extent of the primary x-ray source, but also on the exact position of the collimators that define the field (Scott et al 2009). Thus, any model attempting to reproduce the source occlusion effect should reconstruct both parameters.

Several methods have been suggested for measuring the x-ray source directly. These methods may involve a multiple slit γ camera (Lutz et al 1988), determining the inverse Abel transformation of the derivative of fluence profiles (Treuer et al 1993), output ratios (Zhu et al 1995) or a micro-leaf collimator (MLC) (Treuer et al 2003). Perhaps the most widely known and accurate technique is using a rotating slit collimator and a CT reconstruction algorithm (Munro et al 1988, Jaffray et al 1993, Caprile and Hartmann 2009). Other researchers applied a similar method using a moving slit (Loewenthal et al 1992, Sham et al 2008). However, the experimental set-up of a slit collimator technique is difficult to implement and perform routinely in a clinical environment. Furthermore, Chen et al (2011) exhibited that the accuracy of this technique depends strongly on the choice of slit width, height and distance to the source.

An interesting finding by some of the previous investigations was that the x-ray source may not follow a Gaussian functional form (Sham et al (2008) and Chen et al (2011)). Furthermore, the source distribution may vary with time (Jaffray et al 1993) or as currents in the magnets of the electron beam are altered (Munro et al 1988). Despite the above facts, there is currently no requirement for the physicist to directly measure the source as part of the TPS commissioning or of the linac quality assurance procedure.

The purpose of this work is to suggest a method, including a model and an experimental procedure, that would allow the direct reconstruction of the primary x-ray source distribution as generated in the target using clinical measurements. The model is based on iteratively ray-tracing photons from the target to the measurement plane and vice-versa using a maximum-likelihood expectation-maximization (MLEM) algorithm. No prior assumptions on the source distribution or size are needed. The experimental procedure involves the determination of photon fluence profiles in small fields in air using film measurements and a thin lead (Pb) foil as a build-up material. Blurring effects due to the non-zero electron range and photon scatter are taken into account by convolution with pre-calculated point spread functions (PSF). The method also estimates the appropriate collimator settings that reproduce the source occlusion effect for the given accelerator geometry and measurement set. The accuracy of the method in reconstructing the correct source size and shape is first evaluated by performing MC simulations of the experimental set-up for a range of electron sources and then experimentally by performing measurements on a linear accelerator with known MC beam model parameters.

2.1. The inverse problem

Assuming that the photon fluence distribution can be measured at the bottom of the linear accelerator (exit plane, figure 1(a)), the question we are seeking to answer is how the source fluence distribution can be reconstructed at the top of the target (source plane, figure 1(a)). This question defines an imaging problem to derive an object representation from an initial blurred image. This is referred to as the inverse problem in image reconstruction and iterative methods are commonly applied to address it. In the following, the method will be explained as applied in the source reconstruction problem. In this work we will only consider 1D fluence distributions and the reconstruction is performed on the crossplane and inplane orientations separately.

Zoom In Zoom Out Reset image size

2.2. Extracting the system matrix

2.3. The MLEM reconstruction algorithm

2.4. Experimental set-up

2.5. Blurring in the Pb foil

Even though the high density Pb foil reduces the blurring, due to the smaller electron range, it does not eliminate the effect. Furthermore, photon scattering in the phantom will contribute in penumbra broadening. This blurring can be modeled by a point spread function, $\text{PSF}(x)$. The relative dose in water after the Pb foil, ${{\left({{D}_{w}}\right)}_{pb}}(x)$, can then be expressed as a convolution operation between the relative photon fluence, $\Phi(x)$ and the $\text{PSF}(x)$ (equation (7)).

$$\begin{eqnarray}&&{{\left({{D}_{w}}\right)}_{pb}}(x)=\left(\Phi\ast \text{PSF}\right)(x)=\int_{-\infty}^{+\infty}\Phi(\tau )\text{PSF}(x-\tau )\,d\tau \end{eqnarray} \tag{ 7 }$$

To extract the optimum PSF, an optimization procedure based on MC simulations was followed. First, a Pearson VII functional form was chosen for parameterizing the PSF (equation (8)).

$$\begin{eqnarray}&&\text{PSF}\left(x;\gamma,n\right)={{\left(\frac{{{\gamma}^{2}}}{{{\gamma}^{2}}+\left({{2}^{1/n}}-1\right)\centerdot {{x}^{2}}}\right)}^{n}}\end{eqnarray} \tag{ 8 }$$

As a second step, MC simulations using the accurate accelerator model were performed to calculate the relative photon fluence, $\Phi(x)$, and relative dose in water after the Pb foil, ${{\left({{D}_{w}}\right)}_{pb}}(x)$. The values of γ and n varied in a systematic manner following a brute-force approach. For each set of estimated parameter values ($\hat{\gamma},\hat{n}$) the dose profile was deconvolved with the PSF(x; $\hat{\gamma},\hat{n}$), using a maximum likelihood algorithm, to extract an estimation of the photon fluence, $\rm{\hat \Phi}\left(x;\hat{\gamma},\hat{n}\right)$. The optimum set of parameter values ${{\left(\hat{\gamma},\hat{n}\right)}_{\text{opt}}}$ were determined by minimizing the mean square relative error (MSRE) between estimated and MC calculated photon fluence (equation (9)).

$$\begin{eqnarray}&&{{\left(\hat{\gamma},\hat{n}\right)}_{\text{opt}}}=\underset{\gamma,n}{\mathop{\text{arg}~\text{min}}}\,\frac{1}{N}\underset{i}{\overset{N}{\mathop \sum}}\,{{\left(\frac{{{{\rm{\hat \Phi}}}_{i}}}{{{\Phi}_{i}}}-1\right)}^{2}}\end{eqnarray} \tag{ 9 }$$

The above procedure was also performed for other PSF models, including Gaussian and Lorentzian functions, which presented inferior performance than the Pearson VII. The sensitivity of the extracted PSF to the electron source size used in the MC simulation was evaluated by repeating the above procedure for electron source FWHM values of 0.5, 1.0, 1.5 and 2.0 mm and calculating the average PSF and standard deviation (1 σ). The blurring effect was inherently included in the model by convolving each column of the system matrix with the average PSF.

2.6. Determining the collimator jaw position

2.7. Film measurements

2.8. Monte Carlo simulations

2.9. Method evaluation

2.10. Uncertainty analysis

The calculated dose profiles after 2 mm of Pb, before and after deconvolution with the respective PSF, are presented in figure 2 for the crossplane and inplane orientation. In order to evaluate the accuracy of the kernel in debluring the profile from the electron and photon scattering effects, the profiles are directly compared to the expected photon fluence at a SSD of 105 cm. The calculations were performed using the MC accelerator model for the commissioned values of the electron source FWHM and jaw positions. Figure 3 presents the average PSF for electron sources of FWHM equal to 0.5, 1.0, 1.5 and 2.0 mm along with the 1 standard deviation uncertainty level. The PSFs were calculated using the methodology presented in section 2.5. The average PSF had a FWHM (δFWHM) equal to 0.88 (0.05) mm and 0.85 (0.13) mm at the crossplane and inplane orientations respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The reconstructed sources using the MC calculated dose profiles of the experimental set-up as an input are presented in figure 4. The calculations were performed for a range of typical electron Gaussian sources of FWHM equal to 0.5–2 mm. The respective photon sources at a depth of 0.2 mm in the target are also presented. The Gaussian fits to the reconstructed source distributions exhibited a RMSE (%) in the range of 2.3%–5.4% and of 2.2%–4.1% for the crossplane and inplane orientations respectively. The relative intensity differences between the reconstructed sources and the electron and photon sources (figure 5) presented a local agreement within 10% in all cases with the major discrepancies occurring in regions of high intensity gradient and for source sizes of FWHM equal to 0.5 and 1 mm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ability of the MLEM algorithm to reconstruct the FWHM of the expected source and field size is presented in table 1. The reconstructed FWHM agrees within 0.12 mm and 0.10 mm to the electron source FWHM (incident on the target) and photon source FWHM (at 0.2 mm depth in the target) respectively. The photon source appears broader than the electron source by 0.02–0.04 mm. The reconstructed jaw positions reproduce the expected values with an accuracy better or equal to 0.2 mm.

Table 1. Reconstructed FWHM of the source and reconstructed field size (crossplane (X)  ×  inplane (Y)) using the MC calculated dose profiles as input.

e− source (mm2)	γ source (mm2)	MLEM source (mm2)	Field (mm2)
$0.5\times 0.5$	$0.54\times 0.54$	$0.54\times 0.56$	$4.7\times 4.9$
$1.0\times 1.0$	$1.03\times 1.03$	$1.10\times 1.08$	$4.6\times 4.9$
$1.5\times 1.5$	$1.53\times 1.53$	$1.52\times 1.52$	$4.7\times 4.9$
$2.0\times 2.0$	$2.03\times 2.02$	$2.04\times 2.12$	$4.8\times 4.7$

Note: The expected FWHM of the electron source (incident on the target) and photon source (at 0.2 mm depth in the target) are presented. The expected field size was set to the commissioned values ($4.7\times 4.9$ mm²).

To evaluate the accuracy of the method to reproduce the expected dose profile, the reconstructed dose profiles are directly compared to the MC calculated dose profiles that were used as an input (figure 6). Discrepancies are observed mainly for the smallest source size (FWHM  =  0.5 mm) and in the tail region with a local dose difference reaching 17.7% in the 90–10% dose region.

Zoom In Zoom Out Reset image size

Table 2 evaluates the sensitivity of the method to the choice of the exit plane location. The reconstruction was performed for different SSD values of 105, 125 and 150 cm. The MC electron source and field size were set in this case to the commissioned values of $1.25\times 1.10$ mm² and $4.7\times 4.9$ mm² respectively. The reconstructed source FWHM agreed within 0.04 mm and the reconstructed jaw settings within 0.2 mm.

Table 2. Reconstructed FWHM of the source and reconstructed field size (crossplane (X)  ×  inplane (Y)) for SSD values of 105, 125 and 150 cm.

SSD (cm)	MLEM source (mm)	Field (mm2)
105	$1.22\times 1.08$	$4.6\times 4.9$
125	$1.20\times 1.08$	$4.8\times 5.0$
150	$1.24\times 1.12$	$4.8\times 5.0$

Note: The expected source and field size were set to the commissioned values ($1.25\times 1.10$ mm² and $4.7\times 4.9$ mm²).

The reconstructed source using the film dose profile measurements at the Varian Novalis accelerator is presented in figure 7. In the same figure the Gaussian electron source and the respective photon source, as they were determined during commissioning, are also presented. The reconstructed source exhibited a FWHM of 1.22 mm (±0.12) and 1.21 mm (±0.11) at the crossplane and inplane orientation, respectively. The RMSE (%) of a Gaussian fit was found to be 2.4% and 2.7% for the crossplane and inplane orientations, respectively. The relative local intensity differences between reconstructed source and photon and electron sources are presented in figures 8(a) and 8(b) in the crossplane and inplane orientations respectively. The reconstructed field side for this set of measurements was 4.4 mm on both orientations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The propagated effects of the jaw positioning, experimental and PSF uncertainties on the reconstructed source are presented in figure 9. Overall, the jaw positioning uncertainties resulted in the most significant source variations, while PSF and experimental uncertainties mainly affected the tail region. Table 3 summarizes the reconstructed source parameters and uncertainty components.

Table 3. FWHM and FWTM values of the reconstructed source using the film profile measurements, the MC electron source (incident on the target) and the photon source (at 0.2 mm depth in the target) as determined during model commissioning.

	FWHMx (mm)	FWTMx (mm)	FWHMy (mm)	FWTMy (mm)
e− source	1.25	2.26	1.10	2.00
γ source	1.28	2.32	1.13	2.04
rec source	1.22	2.29	1.21	2.32
${{\sigma}_{\text{total}}}$ (rec source)	0.12	0.27	0.11	0.20
${{\sigma}_{\text{jaw}}}$ (rec source)	0.11	0.21	0.10	0.15
${{\sigma}_{\text{exp}}}$ (rec source)	0.03	0.12	0.04	0.11
${{\sigma}_{\text{psf}}}$ (rec source)	0.04	0.16	0.01	0.07

Note: The total and component uncertainties of the reconstruction are presented at the 1 σ level.

Zoom In Zoom Out Reset image size

Finally, the dose profiles were recalculated using the accelerator MC beam model with the electron Gaussian source set to the reconstructed FWHM value of $1.22\times 1.21$ mm² and field size of $4.4\times 4.4$ mm² (figure 10). The dose distributions exhibited an agreement with the film measurements of 1.2% in the 90–10% dose region.

Zoom In Zoom Out Reset image size

Measuring the photon fluence in air using a 2 mm Pb foil as build-up material proved to be challenging and some penumbra broadening still appeared in the dose region below 40%. For that purpose, the use of a PSF function to account for the electron blurring and photon scattering in the build-up material appeared to be important. In fact, performing the MLEM reconstruction without using the PSF, resulted in an over-estimation of the source FWHM up to 0.42 mm.

The PSFs, extracted by the MC calculations and the optimization procedure (section 2.5), presented minimal sensitivity to the source size selected during the simulations (figure 3). This result ensures that the method can be applied to any accelerator of the same geometry without prior knowledge of the source. The observed variations in the PSF distributions resulted in only small variations of the reconstructed source as can be seen in figures 9(e) and (f). Thus, a pre-calculated average PSF can be provided to the user or even included in the system matrix. In the latter case the blurring is inherently included in the model and the end user only needs to input the measured dose profiles. It should be emphasized that different accelerator designs may present different PSFs. Further research will be needed to verify the applicability of a universal PSF for all accelerator designs.

The MLEM reconstruction using the MC calculated dose profiles with known electron Gaussian sources reproduced the expected Gaussian shape with RMSE values up to 5.4% and 4.1% in the crossplane and inplane orientations respectively. The maximum deviations were observed for the smallest source with FWHM of 0.5 mm. The discrepancy is also depicted in the tail and shoulder regions of the reconstructed dose profiles upon convergence of the MLEM algorithm. The difficulty in reconstructing the source shape in this case could potentially be attributed to the steep change of the dose gradient.

The reconstructed source appears to overestimate the FWHM in most of the cases relative to the expected electron Gaussian source incident on the target (figures 4, 5 and table 1). This can be partly explained by the electron blurring occurring in the target. The photon source at 0.2 mm depth in the target indeed appeared broader up to 0.04 mm. This result is in agreement with the work by Sterpin et al (2011) that reported a broadening up to 0.06 mm for a similar incident energy. This finding implies that the reconstruction algorithm actually reconstructs a distribution closer to the photon source and not to the electron source, even though both are in close agreement to each other.

A second reason may be related to the observation that the overestimation of the source was often coupled with an underestimation of the projected field side. In fact in the cases that the correct field side was reproduced (FWHM of 0.5 mm and 1.5 mm in table 1), the MLEM reconstructed intensity distribution presented an exceptional agreement with the photon intensity distribution (figures 4 and 5) and the FWHM was within 0.02 mm to the expected value. The overall accuracy in the jaw position estimation was found to be 0.2 mm, a result which agrees with the initially estimated precision. The above observations emphasize the need for accurate reconstruction of both the source distribution and jaw position in order to properly predict the source occlusion effect. It should be noted that for the largest source size (FWHM  =  2 mm) the uncertainties in jaw positioning estimation increased and different combinations of source size/jaw positions could still provide acceptable agreement between reconstructed and expected dose profiles.

The experimentally reconstructed source reproduced a Gaussian shape with RMSE (%) values up to 2.7%. The results are within the range of the RMSE (%) values reported previously using the MC simulations of known Gaussian electron sources. Relative to the previously commissioned electron source the reconstructed FWHM was  −0.03 mm lower and  +0.11 mm higher in the crossplane and inplane orientations respectively. Relative to the corresponding photon source the reconstructed FWHM was  −0.06 mm lower and  +0.08 mm higher in the crossplane and inplane orientations. In both cases, the reconstructed source agrees within the estimated total uncertainty level (1 σ), even though only marginally for the inplane orientation (figure 7). More importantly, if the reconstructed source FWHM and field size are used as an input to the MC accelerator beam model, an excellent agreement was observed between measured and calculated dose profiles (figure 10).

The uncertainty component analysis, presented in figure 9 and table 3, indicates that the estimation of the jaw position is the major source of uncertainty during reconstruction. A 0.2 mm misestimation of the collimator position would result in about 0.10–0.11 mm misestimation of the source FWHM. This, in turn, would result in output factor variations of about 1.5–3.0% for typical source sizes in the range of 0.8–1.4 mm and penumbra width (80–20%) variations of 3–4%. This level of accuracy still competes with most small field detectors that are currently used in the clinic. The development of independent methods for estimating the correct jaw positions would greatly improve the performance of the source reconstruction algorithm. Uncertainties due to PSF variations with the source and experimental measurements, including mechanical jaw positioning, were less significant in the prediction of the FWHM. However, these effects presented a much more significant impact on the tail regions of the source distribution as illustrated in the reported FWTM uncertainties (table 3).

In this work, we investigated the performance of a MLEM-based source reconstruction technique for clinical linear accelerators using small field photon fluence profiles. The use of a high density build-up material along with an appropriate PSF were found important for extracting accurately the photon fluence. The PSFs exhibited overall a minimal dependence on the source. The model was able to reconstruct the electron source with an overall accuracy of 0.12 mm for typical source sizes with FWHM ranging from 0.5 to 2 mm that were tested in simulations. Experimentally, the method exhibited an overall accuracy level of 0.11 mm with respect to a previously commissioned MC model, with the most significant uncertainty attributed to the estimation of the jaw position at the time of measurements. The results of this study indicate that MLEM-based approaches could be a powerful and practical tool for the direct reconstruction of the source distribution, without any prior assumptions. In the future, such techniques could become part of a quality assurance procedure that would evaluate potential source variations through time and machine usage or between different accelerator types.

PP would like to acknowledge Dr A Ahnesjö, Dr A Reader and S Lee for important conversations on the subject as well as Compute Canada/Calcul Quebec for providing computing resources. PP gratefully acknowledges the financial support by the A S Onassis Public Benefit Foundation in Greece and by the CREATE Medical Physics Research Training Network grant of the Natural Sciences and Engineering Research Council (Grant number: 432290). IL acknowledges research funding from the Research Institute of the McGill University Health Centre, the Montreal General Hospital Foundation, and the Phil Gold Fellowship.