Dose influence matrix calculation
matRad's dose influence matrix calculation algorithms are split into two parts: First we determine the irradiation geometry by generating the steering information for the desired beam setup. In a second step, we generate dosimetric information by pre-computing dose influence matrices for inverse planning.
<a name="stfInfo" /a> #Steering information
The steering information holds all geometric information about the irradiation within the stf struct. It is generated by calling matRad_generateStf.
The stf struct uses both the LPS coordinate system and a beam's eye view coordinate system in [mm]. Corrdinates in beam's eye view coordinate system are labeled accordingly; coordinates in LPS coordinate system do not have an extra label.
<a name="rayBixelConcept" /a> ##Ray and bixel concept
The irradiation geometry is organized to a ray and bixel concept which is schematically shown below.
Schematic visualization of the ray and bixel concept
From a virtual radiation source (yellow) the target volume (red) within the patient (green) is covered by equidistant rays (black). Note that only a two dimensional cut through a three dimensional cone of rays is shown for clarity. In the iso center plane (not shown) the distance of the individual rays correpsonds to the bixel width (compare pln struct. For photons, the term bixel refers to a discrete rectenagular fluence element (the limits of the individual bixels are shown in dashed black). Together, all bixels cover the entire target volume.
For 3D IMPT for particles, we have an additional degree of freedom, namely the particle energy to be considered. This is accounted for during the stf-struct generation by determining the depth of the target volume on individual rays and placing spots (black dots) accordingly.
More information about the ray and bixel concept (though with slight variations in nomenclature) can be found in sections 2.3 and 2.5 Nill (2001).
<a name="doseInfData" /a> #Dose influence data
Based on the steering information, matRad can compute dose for individual photon bixels or proton/carbon ion pencil beams. The required geometric and radiological distances facilitate the same set of matRad functions.
<a name="rayCasting" /a> ##Ray casting
We perform an exact calculation of the radiological depth on the central axis for every bixel or pencil beam with a ray tracing algorithm in matRad_siddonRayTracer according to Siddon (1985) Medical Physics.
Geometric distances from the virtual radiation source and lateral distances to the central ray are computed by standard matrix vector algebra in matRad_calcRadGeoDists.
<a name="photons" /a> ##Photons
<a name="photonDoseCalc" /a> ###Dose calculation algorithm
matRad uses a singular value decomposed pencil beam algorithm implemented according to Bortfeld et al. (1993) Medical Physics. The functionalities are organized in the matRad functions matRad_calcPhotonDose and matRad_calcPhotonDoseBixel.
The dose delivered to a certain voxel i from bixel j is stored as dose influence matrix in the dij struct using matlab's built-in double precision sparse matrix format.
<a name="photonBaseData" /a> ###Base data
The necessary measured base data, namely the kernel functions as described by Bortfeld et al. (1993) Medical Physics are supplied for a 6MV linac and stored in photonPencilBeamKernels_6MV.mat as matlab piecewise polynomial.
<a name="photonApprox" /a> ###Approximations
For the photon dose calculation we assume a uniform primary fluence. While this should not have any impact on the resulting dose distributions, as we allow for intensity-modulation, this reduces the computation time because we only have to perform one convolution and we can compute the individual contributions from different bixels with the same kernel martcies (kernel1Mx, kernel2Mx, kernel3Mx). An extension to inhomogeneous primary fluences is easily possible.
matRad only performs a ray tracing on the central ray of every bixel. For off axis voxels, the projected radiological depth on the central axis is used. While a ray tracing for every individual voxel would not increase computation time considerably (as the ray tracing could be re-used for all bixels of one beam direction) we chose this approach for consistency with particle dose calculation. Note that this approximation is naturally deteriorating with increasing bixel width; if you want to modify matRad's photon dose calculation algorithm for larger fields you would have to adjust the ray tracing strategy accordingly.
<a name="photonCompBottlenecks" /a> ###Computational bottlenecks
For the photon dose calculation there are three main possibilities to speed up computations and reduce memory consumption
- It is possible to reduce the spatial resolution of the dose calculation in the patient CT by downsampling the CT data upon import. If you do this as a post-processing step, be aware you need to adjust the binary segmentations in the cst cell array accordingly.
- You can increase the variable pln.bixelWidth in the matRad script which effectively reduces the number of bixels which have to be computed approximately quadratically.
- You can reduce the radius around the central ray, around which dose in computed by adjusting the variable lateralCutOff in matRad_calcPhotonDose. Currently this variable is already set to 20mm which is quite tight so we do not recommend this.
<a name="particles" /a> ##Particles
<a name="particleDoseCalc" /a> ###Dose calculation algorithm
matRad facilitates a simple pencil beam model for particle dose calculations similar to the work of Hong et al. (1996) Pysics in Medicine and Biology. The dose at a particualr voxel is given as the product of a depth dependent part and a lateral part. For the depth dependent part, matRad uses tabulated depth dose curves for individual particle energies. For lateral beam broading, matRad uses a depth-dependent sigma of a Gaussian profile, which is also tabulated versus depth for all available beam energies.
The dose delivered to a certain voxel i from bixel j is stored as dose influence matrix in the dij struct using matlab's built-in double precision sparse matrix format.
<a name="particleAlphaBeta" /a> ###α- and β-matrix pre-computations
For carbon ions, matRad also enables the computation of α- and β-matrices that can be used to compute three-dimensional dose weighted α- and β-distributions which can in turn be used to compute three-dimensional RBE distributions during inverse planning.
matRad only models variations of α and β with depth in matRad_calcLQParameter. Potential dependencies in lateral direction are not explicitly modeled.
<a name="particleBaseData" /a> ###Base data
The base data files protonBaseData and carbonBaseData required for particle dose calculation include depth dose curves and tabulated lateral beam widths (as Gaussian sigmas) for a library of different energies. The proton base data has been computed based on an analytical approximation for the Bragg curve and Highland's approximation for multiple Coulomb scattering. The carbon ion base data has been Monte Carlo simulated for an idealized beam line without monitoring devices. Besides this physical information, the carbon ion base data also includes α and β tables that have been computed based on LEM IV. Within the *.mat files, the depth is stored in [mm], α tables are stored in [1/Gy], β tables are stored in [1/Gy^2], and the integrated depth dose distribution is stored in [MeV cm^2 /(g * primary)]. Upon dose calculation, the integrated depth dose is converted to [Gy mm^2 /(1e6 primaries)] with a linear scaling in the function matRad_calcParticleDoseBixel.m.
<a name="particleApprox" /a> ###Approximations
Besides the standard approximations made by pencil beam algorithms (i.e. factorization of lateral and depth dependent part, ray tracing only on central ray), we do not make any approximations for particle dose calculation.
<a name="particleCompBottlenecks" /a> ###Computational bottlenecks
For the particle dose calculation there are three main possibilities to speed up computations and reduce memory consumption
- It is possible to reduce the spatial resolution of the dose calculation in the patient CT by downsampling the CT data upon import. If you do this as a post-processing step, be aware you need to adjust the binary segmentations in the cst cell array accordingly.
- You can increase the variable pln.bixelWidth in the matRad script which effectively reduces the number of pencil beams which have to be computed approximately quadratically.
- You can reduce the radius around the central ray where dose is computed by adjusting the variable lateralCutOff in matRad_calcParticleDose.
