A dose error evaluation study for 4D dose calculations

Previous studies have shown that respiration induced motion is not negligible for Stereotactic Body Radiation Therapy. The intrafractional breathing induced motion influences the delivered dose distribution on the underlying patient geometry such as the lung or the abdomen. If a static geometry is used, a planning process for these indications does not represent the entire dynamic process. The quality of a full 4D dose calculation approach depends on the dose coordinate transformation process between deformable geometries. This article provides an evaluation study that introduces an advanced method to verify the quality of numerical dose transformation generated by four different algorithms. The used transformation metric value is based on the deviation of the dose mass histogram (DMH) and the mean dose throughout dose transformation. The study compares the results of four algorithms. In general, two elementary approaches are used: dose mapping and energy transformation. Dose interpolation (DIM) and an advanced concept, so called divergent dose mapping model (dDMM), are used for dose mapping. The algorithms are compared to the basic energy transformation model (bETM) and the energy mass congruent mapping (EMCM). For evaluation 900 small sample regions of interest (ROI) are generated inside an exemplary lung geometry (4DCT). A homogeneous fluence distribution is assumed for dose calculation inside the ROIs. The dose transformations are performed with the four different algorithms. The study investigates the DMH-metric and the mean dose metric for different scenarios (voxel sizes: 8 mm, 4 mm, 2 mm, 1 mm; 9 different breathing phases). dDMM achieves the best transformation accuracy in all measured test cases with 3–5% lower errors than the other models. The results of dDMM are reasonable and most efficient in this study, although the model is simple and easy to implement. The EMCM model also achieved suitable results, but the approach requires a more complex programming structure. The study discloses disadvantages for the bETM and for the DIM. DIM yielded insufficient results for large voxel sizes, while bETM is prone to errors for small voxel sizes.


Introduction
The technical opportunities in the field of precise dose calculation and accumulation of dose inside the dynamic anatomy of the human body are not yet fully exploited in clinical practice Zhong 2008, Jaffray et al 2010). Dose accumulation is applicable in several tasks of Stereotactic Body Radiation Therapy (SBRT). Dose summation is well known for different beam angles or the cumulative effect of beamlets in the case of a dynamic IMRT. However, all these applications are based on a static geometry. For lung cancer treatments, dose accumulation has to consider the full deformation process inside the lung and the abdomen. Due to breathing induced motion, the model has to generate several dose distributions which are all based on separate single geometries (e.g. consecutive CT series in a 4DCT scan). A deformation grid (Söhn et al 2008) links these geometries in an mathematical way. The result is a group of vector fields which represents the spatial relationships. This allows transformations of dose onto an alternative geometry. The definition of a reference set enables the realization of a total accumulated static dose distribution. The basic 4D accumulation principle is well established and verified with different applications by Rosu et al (2006), Keall et al (2004), Janssens et al (2009) and Söhn et al (2009).
In recent years a number of approaches have been established that implement single dose transformations. A milestone for precise dose transformations is the energy model and its adaptions invented by Siebers and Zhong (2008) and Zhong and Siebers (2009). A more elegant approach is the direct voxel tracking developed by Heath and Seuntjens (2006) and Heath et al (2011). However, due to the complexity of the accumulation task, a list of fundamentally different solutions were published (Rosu and Hugo 2012) with more or less accurate comparative results. Heath et al (2008) and Yan et al (2012) offered methods to calculate the dose error for different dose transformation models. Li et al (2014) and Heath et al (2009) published studies that compare the fundamentals of dose mapping and energy transformation. Based on these evaluations, this work tries to elaborate on more differences to highlight advantages and disadvantages of the respective models. Therefore, this paper introduces an advanced dose transformation algorithm, called divergent dose mapping model (dDMM). The study compares a selection of established approaches. Furthermore, it introduces an accurate dose validation metric to verify the quality of the performed dose transformation. The used validation metric is a refinement of the metric used by Yan et al (2012) and depends on the conservation of the dose mass histogram (DMH) to measure the dose transformation quality.
The main objective of this manuscript is the comparison of dose transformation algorithms by using several comparison metrics. The results during the study contributed to the development of dDMM. Its introduction as an advanced solution of DIM is a secondary objective.

Description of inhomogeneous dose distributions
Doing the accumulation of dose in a entirely numerical procedure requires the discretization of inhomogeneous dose distributions during calculation. The effect is well described in the dose definition for a non-uniformly and partially irradiated object (Brahme 1984, ICRP 2007. The calculation of the absorbed dose for a macroscopic object is accompanied with the dose approximation of several point doses to a total dose value D . The result is the mass weighted average of the involved dose points (Brahme 1984) ⃗ with m being the mass. ρ ⎯→ x ( ) and ⎯→ D x ( ) represent the mass density and the dose for a coordinate ⃗ x inside the investigated volume V. A discrete formula is necessary for dose calculations. Such an approach considers unique voxels x with related masses m(x) and energy values E(x) that are deposited in the respective voxel: (2)

Dose accumulation
The goal of 4D dose accumulation Zhong 2008, Rosu and is the summation of different dose distributions. These distributions are based on time dependent geometries. A series of different breathing phases j that cover the full breathing cycle has to be generated. One distinguishes between the reference set i and the following breathing phases j. Phase i is is the master geometry, or reference set and it collects the dose distributions arising from all other phases j. i is connected with the deformation grid (Siebers and Zhong 2008) v ji to every breathing phase j. This enables coordinate transformations to each time set in the breathing cycle. Dose accumulation can be written as: D i∑ is the summed dose in the reference set. T j is the transformation operator representing the special transformation model. D i∑ is directly defined by a summation of single dose transformations related to all breathing phases n (Söhn et al 2009).

Models of dose transformation
Previous studies regarding 4D dose accumulations established a set of different transformation models (Heath and Seuntjens 2006, Siebers and Zhong 2008, Zhong and Siebers 2009, Yan et al 2012, Rosu and Hugo 2012. The investigation of the used algorithms reveals two elementary branches (Li et al 2014). Almost every popular method is associated with one of these two groups. As a matter of fact, there exist two representations of the total dose defined by (2). Here, D is expressible either with the mass weighted average of the dose or with the ratio of absorbed energies per mass. Therefore, one can derive two groups for dose accumulation algorithms: The Dose Mapping Model (DMM), which deals with the transformation of the applied dose distributions and the Energy Transfer Model (ETM), which computes the dose mapping with the transfer of deposited energy values. Figure 1 presents an overview of algorithms related to one of the basic approaches.
2.3.1. Dose interpolation method (DIM). DIM follows the conventional rules of trilinear interpolation (see figure 1(a)) and has been well described in the field of 4D dose accumulation (Rosu et al 2005, Rosu and). (dDMM). This paper introduces dDMM to fulfill the requirements of the dose definition (1), which are not considered by DIM, i.e. it considers the mass fractions of single voxels for the mapping process. While DIM calculates every single voxel x i starting inside the reference set i with projection v ji (x i ) on the breathing state j, dDMM goes the opposite way v ij (x j ). For this reason, all dose voxels x j in j are separated into a sufficient number of sub-voxels x j sub . Every sub-voxel x j sub receives the dose D j (x j ) of its parent and is composed by the fractional mass of its parent The dose D x ( ) j j sub of a sub-voxel is weighted by its mass m x l ( ) j . The calculation is performed discretely using a high resolution deformation grid v x ( ) ij j sub . The algorithm determines one target x i inside the reference set for every sub-voxel x j sub . The dose D i (x i ) in a target voxel is normalized to the set of all sub-masses that migrate to x i . A single dose transformation can be written as:

Divergent dose mapping model
with W xi being a mathematical set which unites all sub-voxels x j sub that have a connection to x i . This results in the mass weighted average of all involved dose values D x ( ) j j sub for D i (x i ). dDMM compensates the mass effect (neglected in DIM) and avoids this loss of already calculated dose information. The schematic work flow of dDMM is illustrated in figure 1(b). (2008) and Rosu and Hugo (2012). The fundamental idea is to transform energy events with real coordinates ⃗ x j instead of the numerical dose voxels based on voxel indices x j (see figure 1(c)). Therefore, all necessary energy events →

( )
W E x i in j are collected. They are connected to the destination voxel x i . Its dose is given by: 2.3.4. Energy mass congruent mapping (EMCM). The EMCM model was introduced by Zhong and Siebers (2009). The model splits doses in energy and mass, but uses the same components for energy deposition and for dose calculation. For this purpose, the basic model (bETM) is used, but the final dose calculation inside the reference geometry i is performed with a mass mapped structure of j (see figure 1(d)). This should prevent possible dose peaks that are generated by bETM with defective deformation grids inside mass heterogeneous structures. For example (using bETM), a wrong transformation vector connecting the target voxel x i and the source voxel x j results in a large dose error, if their masses differ a lot. The result could be an extreme dose overestimation, if the source voxel receives a large energy based on a heavy mass value and the target mass voxel value is comparatively low. In practice, a numerical deformation grid for real data will never be perfect. Hence, energy mass displacements can not be completely prevented. However, the use of a mass mapped structure of j in i for dose calculation should reduce these errors.

Deformation grid and patient data
This evaluation study uses the data and the deformation grid of the POPI-model, a Pointvalidated Pixel-based Breathing Thorax Model. It is an established scientific data set with lung cancer patient images containing a full 4DCT. It is provided by the Léon Bérard Cancer Center and the CREATIS Laboratory (Lyon, France). The model is fully described by Vandemeulebroucke (2007) (www.creatis.insa-lyon.fr/rio/popi-model). This study uses the 4DCT series of the preprocessed images containing ten static CT series. The transformation is based on the deformation data (v ij ) processed by the parametric deformation method of the POPI-model. The various mapping methods used in the study require deformation grids which operate in opposite directions. Inverse vector fields i , which are necessary for various transformation algorithms, are approximated with the Newton-Raphson optimization implemented in the iPlan™ RT framework. A validation of the inverse fields revealed a mean deviation of 0.003 mm and a maximum deviation of 0.31 mm. Barring DIM all methods use the inverse vector fields. Hence, almost all transformations are based on the same deformable image registration (DIR) data. Thus, errors based on the DIR algorithm are negligible for this study.

Theory of evaluation
To measure the quality of the explained mapping procedures the mapping process is repeated for a set of arbitrary regions of interest (ROI). The analyzed volumes are small cuboid ROIs (Δx = 32 mm, Δz = 24 mm, Δy = 16 mm, inspired by Yan et al (2012)) distributed over the full 4D geometry. The large number of test cases m = 900 (value of considered ROIs) enable a qualitative assessment. The dose calculation (see section 2.8) and the dose transformation process is performed for every ROI. The results of the different transformation models are compared by error values. The investigation is focused on the transformation itself. Effects of the deformation grid, the dose algorithm and any other sources of error are neglected. Figure 2 illustrates the main idea: The first step is the selection of a breathing phase j. The random ROI in j provides the volume for the dose calculation D j (x, y, z). The subsequent dose transformation is performed with all transformation models: DIM, dDMM, bETM, EMCM. This step generates four different target distributions D i (x, y, z). The determination of the mapping error η enables comparative studies. The procedure is carried out for four different voxel sizes. Then, the same calculation is repeated for the next breathing phase until the full breathing cycle is reached. A total of 900 arbitrary ROIs is evaluated.

An error metric based on the dose mass histogram
The quality of a treatment plan cannot be determined by an integral dose or a single dose value. It is defined by the entire volume based dose distribution. Dose gradients and mass inhomogeneities lead to different mass dose fractions inside the volume. The transformation process should not alter this dose information, otherwise the accuracy of the basic dose calculation is obsolete. In clinical practice, one of the most common tools for plan assessment is the dose volume histogram (DVH). However, the DVH assumes a constant mass inside all voxels. This is sufficient for homogeneous organs, but not for heterogeneous objects as considered in this work. The dose mass histogram (DMH) is an extension of the DVH. It is a more accurate approach. The DMH concept uses the fractionated mass that absorbs a certain dose. Wei (2005), Mavroidis et al (2008) and Nioutsikou et al (2005) investigated the DMH approach for lung complications in comparison to the DVH. The main conclusion was a dose overestimation associated to the DVH concept. In contrast, the effectiveness of the DMH model was closer related to the radiation effects. Furthermore, deformation and expansion of anatomical volumes hinder the use of the DVH concept in this task. Therefore, a promising tool for the transformation assessment is the analysis of the integral ROI based DMH. The DMH value should not change during dose transformation, because this guarantees the conservation of each dose mass fraction. Due to numerical effects, an exact DMH conservation is not feasible. A good approach should approximate the unmapped DMH as close as possible.
It is not suitable to compare the differential DMH. Voxel unions could cause small dose shifts for internal sub-volumes. These shifts would cause the same error as model based dose displacements. In contrast, the cumulative DMH is resistant to this effect. Hence, it is necessary to analyze the cumulative DMH: The DMH error for a specific dose D is the difference between the mapped (i) and the unmapped cumulative DMH (j). It is defined by: It describes the mass discrepancy of both dose distribution with respect to a specific dose D. A cumulative error parameter (over the full dose range 0 ... D max ) can be defined to compare the full dose distributions: It is useful to normalize the integral dose difference for a qualitative assessment and to distinguish the outcome of different transformation models and ROIs with different mass fractions. Hence, the total normalized DMH error η DMH is defined by: For a better understanding of this approach figure 3 graphically illustrates the cumulative DMH and the DMH error for one sample.
2.6.1. Technical remarks. Numerical issues are responsible for total mass deviations at D = 0% for DMH(D) in figure 3 (left). The deformation process leads to the fact that more The final doses D (dose-to-medium, after dose transformation) are calculated either with m i (voxel masses of the reference set i) for bETM, DIM, dDMM or with m j → i (mapped masses of j inside the reference set i) for EMCM. Hence, the model related DMHs have to be calculated differently. EMCM uses m j → i instead of m i for DMH determination. This leads to small deviations of the total masses at D = 0% shown in the DMH graph of figure 3.

An error metric based on the mean dose
In some cases it is also useful to analyze the mean dose error η < D > . This value is based on the ROI related average dose. Only voxels x j are considered that count to the total number l of exposed voxels inside the ROI. The dose average is written as: η <D> compares the mean dose before ( j) and after (i) dose transformation. Again, the normalization with the untransformed value has to be considered. η < D > is defined by: η < D > is less meaningful than η DMH . It serves not for global conclusions. The parameter was mainly introduced to filter a dose peak effect. The effect is generated by faulty energy mass displacements (see section 2.3.4). The error is not easily observable in the DMH graph, because it is mainly caused by very small mass values.

Dose algorithm
The study uses a simplified Monte Carlo approach (dose-to-medium) to calculate the dose inside the arbitrary ROIs. To simulate a beam that would deliver a homogeneous dose distribution in water, a spatially homogeneous and isotropic photon particle fluence Φ is assumed in the entire ROI. Furthermore, the energy spectrum of the investigated photons is mono-energetic containing exclusively 6 MeV photons. In general, the fluence is defined as: ROI (16) with N being the number of particles crossing an area A and ∑dl being the sum of all photon path lengths that traverse the defined volume dV ROI . The dose algorithm generates a determined number of photon particles N. Therefore, specific coordinates ⎯→ x of their traversing paths ∑dl are drawn randomly inside the volume of the ROI. These photon coordinates are randomly distributed based on an equal distribution to simulate the homogeneous fluence. Furthermore, the algorithm calculates the probability for an interaction with matter in this coordinate. Derived from the probability density function of the attenuation law, the length of the free photon path length s is inversely proportional to the mass attenuation coefficient μ of the investigated photon energy. In reverse, a higher mass attenuation coefficient leads to more interaction due to shorter free path lengths. This dose algorithm does not distinguish between different interaction types. Also secondary particles are neglected. Hence, the mass energy absorption coefficient μ en given by Hubbell and Seltzer (1995) is used to estimate the deposited energy. The particle interaction and the deposited energy (dose-to-medium) depends on the local absorption coefficient μ ⎯→ x ( ) en . For small coefficients, the particle traverse the volume without interaction with a higher probability. The resulting dose distributions show dose gradients that are merely based on mass density variations inside the ROI (see figure 4). The conversion of Hounsfield Units (HU) to mass density values and to the atomic composition is done using the XVMC code (Fippel and Nüsslin 2001). All energy events occur point wise in single coordinates. A point energy deposition event is not realistic, but sufficient for this application. Furthermore, it simplifies the analysis of various ETM models. For each dose calculation N = 2· 10 5 uniformly distributed photon particles are considered within the ROI. Due the small size of the ROI, the approach yields a relatively realistic dose distribution illustrated by figure 4. However, the aim is not to generate particularly accurate dose distributions, but to investigate the dose transformation process. For example, figure 4 illustrates the difference of a sample dose distributions before (c) and after (d) dose transformation (EMCM) with a specified vector field.

Full implementation and parameter variations
The source code for the different transformation methods was written using the software package IDL (Interactive Data Language from ITT, Version 7.1). In detail, the simulation generates 900 ROIs. The position of the volumes is equally distributed inside the patient geometry. When a sample ROI is selected inside the patient, the evaluation algorithm chooses a coarse voxel size of 8 × 8 × 8 mm 3 . The procedure generates a dose calculation for every breathing phase (in total 9 phases) except the reference set. Thereafter, the actual dose transformation is performed using the mentioned transformation models (DIM, dDMM, bETM, EMCM). The transformation to the reference geometry allows to determinate the mapping error (η DMH , η < D > ). After evaluation, the algorithm repeats the full procedure with finer voxel sizes with side lengths of 4 mm, 2 mm and 1 mm. This procedure is repeated for all ROIs (see figure 2).

DMH error regarding different voxel sizes
The comparison and the assessment of the quality regarding different transformation models (DIM, dDMM, bETM, EMCM) is the main goal of this evaluation. Quality limitation for dose transformations could lead to numerical or model specific inaccuracies. Hence, the analysis is performed with different voxel sizes. The results are illustrated in figure 5. The graph displays the mapping error η DMH as a statistical box-plot. The results of different models reveal significant differences. dDMM generates the best outcome (exception: DIM for 1 mm) visible for all statistical properties (median, mean, all percentiles). The median is at least 3-4% better than the median of the other models. The same holds for the average and the upper quartile. For the lower quartile and the minimum values exist even greater differences. The lower quartile error of dDMM is up to 4-8% lower than the relative value of the energy models (bETM, EMCM). Hence, the most suitable method for DMH conservation is the dDMM approach. EMCM achieves also accurate results. The statistical properties (median, mean, all percentiles) are up to 2-3% better than the values of bETM. Every model shows a continuous improvement for all statistical parameters (median, mean, percentiles, variance) with decreasing voxel size, which is generally expected. Low resolutions typically create a worse outcome due to numerical effects. However, the variance for large voxel sizes is particularly large for the DIM approach. While it achieves even the best results for 1 mm voxel size, the DMH error goes up to 80% for other sizes. The reason for this behavior is that DIM neglects the mass of the dose voxel during dose transformation. Figure 6 illustrates the breathing phase related results for η DMH . Every sample group is quite small (158 ROIs), because the data derives from a fixed voxel size (2 mm) and a special destination area inside the lung. Here, only ROIs are considered which are located in an area of homogeneous density inside the lung. The homogeneous area ensures only moving ROIs for testing. The supposed dependency is observable for bETM. The large displacement is associated with a large error. The maximum error for bETM is located at the maximal amplitude 60% . The model inherits a high deformation grid dependency. A large displacement leads to more defective grids in the vector field. This increases the probability of energy mass displacements (see section 2.3.4). EMCM, DIM and dDMM achieve constant error results for the full breathing cycle. They are resistant to large energy mass displacements. Indeed there are small fluctuations visible for dDMM, but these results may be explained by the small sample group of this test.

Investigation of the mean dose error
In previous tests, the effect of energy mass displacement (see section 2.3.4) was barely demonstrable. This could be explained with very small masses that were affected. Those parts have a low influence on the DMH error due to the low mass weighting. In contrast, these parts are normally processed with η < D > . The results regarding η < D > are illustrated in figure 7 for different voxel sizes. The main effect is especially apparent for the voxel size of 1 mm. The mean value of bETM is quite large. The maximum values are also larger for bETM in comparison to any other model (not visible in figure 7). The test confirms the effect of energy mass displacements. This explains the limited efficiency of bETM and EMCM compared to DIM and dDMM regarding small voxel sizes. In summary, dDMM achieves the best score for all voxel sizes.

Discussion
A very important goal associated with dose accumulation is the numerical unification of inhomogeneous dose distributions. The mass weighted average derived from the definition in (1) is a good approximation to describe inhomogeneous doses in the context of dose accumulation. Different implementation strategies and algorithms are presented to realize the discrete dose transformations. DIM neglects the mass effect, which was demonstrated by large errors for large voxel sizes. The energy models (bETM, EMCM) fulfill the requirements in an indirect way due to the dissection of the dose in masses and energy shares, but this can lead to undesirable effects especially for small voxel sizes. Here, we introduced the advanced method dDMM that directly implements the dose definition according to (1), almost without numerical restrictions due to the high resolution of the sub-voxel grid. In all tests dDMM yields better results in comparison to any other model. Here, the error values are typically 3-5% better than the results generated by the energy models especially for small voxel sizes.
bETM's susceptibility to errors based on energy mass displacements (section 2.3.4) is especially visible for tests using the mean dose error η <D> . Beside dDMM, also EMCM presents suitable results. However, the clinical application gets complicated due to the increased number of computational steps, which are necessary to build the mass mapped structure. Furthermore, any energy model should be used carefully. The tests of this study were performed with a simplified dose algorithm that realizes coordinate based energy deposition points. A full Monte Carlo algorithm typically performs primary energy depositions along discrete path lengths. This discrepancy and a reduction of path lengths to single coordinates, which are needed for bETM and EMCM, could cause additional errors. These influences are not considered yet and have to be investigated further.
It is difficult to compare the findings of this manuscript with the results of previous studies. Since the DMH error metric was used for the first time, previous work may not detected energy mass displacements even though the effect occurred. Instead, other studies more often focused on integral dose metrics composed by summarized energy values. From that point of view regarding their specific definitions, among other metrics (Yan et al 2012) discovered an average accumulated mean dose mapping error of 5%. Li et al (2014) measured mean dose errors up to 11% for specific target volumes between their implementations of DIM and bETM. However, the mentioned studies used error metrics based on integral dose or integral energy values regarding the entire ROI. Hence, no spatial energy-mass relations were considered. This prevented conclusions about specific energy mass displacements in the past. Future calculations should be handled carefully regarding quality assessment and error determination.
We want to point out that the conclusions drawn in this paper are only proven in the clinical scenario of lung deformation. We can assume similar or less pronounced effects for other moving organs, e.g. breathing induced liver distortions. However, the impact of different mass heterogeneities and structures of different compressibility should be analyzed in further studies, which is beyond the scope of this paper.

Conclusions
This study presented tests to investigate the capabilities of several dose transformation models for dose accumulation based on different breathing phases for lung cancer patients. DIM is unsuitable for large voxel resolutions. bETM is prone to errors especially for energy mass displacements on high voxel resolutions with heterogeneous mass distributions. dDMM as well as EMCM are best suited for practical applications. Their algorithmic structure is consistent with the dose definition and empirical tests (η DMH and η < D > ) confirmed their practicality. EMCM is associated with a high technical effort, while dDMM is more straightforward and efficient.