|
| 1 | +.. _dosecalc: |
| 2 | + |
| 3 | +----------------------------------- |
| 4 | +Dose (influence matrix) calculation |
| 5 | +----------------------------------- |
| 6 | + |
| 7 | +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. |
| 8 | + |
| 9 | +Steering information |
| 10 | +-------------------- |
| 11 | + |
| 12 | +The steering information holds all geometric information about the irradiation within the :ref:`stf struct <stf>`. It is generated by calling :func:`matRad_generateStf`. |
| 13 | + |
| 14 | +The :ref:`stf struct <stf>` uses both the LPS coordinate system and a beam's eye view coordinate system in [mm]. Coordinates in beam's eye view coordinate system are labeled accordingly; coordinates in LPS coordinate system do not have an extra label. |
| 15 | + |
| 16 | +Ray and bixel concept |
| 17 | +~~~~~~~~~~~~~~~~~~~~~ |
| 18 | + |
| 19 | +The irradiation geometry is organized to a ray and bixel concept, which is schematically shown below. |
| 20 | + |
| 21 | +*Schematic visualization of the ray and bixel concept* |
| 22 | + |
| 23 | +.. image:: /images/rayBixelConcept.png |
| 24 | + |
| 25 | +From a virtual radiation source (yellow) the target volume (red) within the patient (green) is covered by equidistant rays (solid black). Note that only a two-dimensional cut through a three dimensional cone of rays is shown for clarity. In the isocenter plane (not shown) the distance of the individual rays corresponds to the bixel width (compare :ref:`pln struct <pln>`). For photons, the term bixel refers to a discrete rectangular fluence element (the limits of the individual bixels are shown in dashed black). Together, all bixels cover the entire target volume. |
| 26 | + |
| 27 | +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. |
| 28 | + |
| 29 | +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) <http://archiv.ub.uni-heidelberg.de/volltextserver/1802/>`_. |
| 30 | + |
| 31 | +Dose influence data |
| 32 | +------------------- |
| 33 | + |
| 34 | +Based on the steering information, matRad can compute dose for individual photon bixels or particle ion pencil beams. The required geometric and radiological distances facilitate the same set of matRad functions. |
| 35 | + |
| 36 | +Dose Calculation algorithms |
| 37 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 38 | + |
| 39 | +Starting with Version 3.1.0, matRad switched to an object-oriented Dose Engine concept. The dose engines are stored :mod:`matRad.doseCalc.+DoseEngines`, and all dose calculation algorithms are implemented as subclasses of the abstract class :class:`DoseEngines.matRad_DoseEngineBase`. |
| 40 | +These engines can be then used directly for configuration and dose calculation, or configured in ```pln.propDoseCalc`` to then be instantiated in the top-level functions :func:`matRad_calcDoseInfluence` or :func:`matRad_calcDoseForward` (see the page on the :ref:`pln-struct <pln>` for more details). |
| 41 | + |
| 42 | +.. code-block:: matlab |
| 43 | +
|
| 44 | + % Example of how to use the dose engine directly |
| 45 | + doseEngine = DoseEngines.matRad_PhotonPencilBeamSVDEngine(pln); |
| 46 | + doseEngine.enableDijSampling = false; %Configure a property of the dose engine |
| 47 | + dij = doseEngine.calcDoseInfluence(ct,cst,stf); |
| 48 | +
|
| 49 | + % Example for a configuration in pln.propDoseCalc |
| 50 | + pln.propDoseCalc.engine = 'SVDPB'; %The engine is equivalent to the DoseEngines "shortName" |
| 51 | + pln.propDoseCalc.enableDijSampling = false; %Configure the property in pln.propDoseCalc |
| 52 | + dij = matRad_calcDoseInfluence(ct, cst, stf, pln); |
| 53 | +
|
| 54 | +Dose Engines implement the abstract property :attr:`DoseEngines.matRad_DoseEngineBase.possibleRadiationModes` to tell their availability for the certain radiation modalities. |
| 55 | + |
| 56 | +Given a ``pln`` with set ``pln.radiationMode`` (and optionally ``pln.machine``), the static function :meth:`DoseEngines.matRad_DoseEngineBase.getAvailableDoseEngines` can be used to retrieve a list of available dose engines for the given radiation mode and machine. |
| 57 | + |
| 58 | + |
| 59 | +External Beam Therapy |
| 60 | +""""""""""""""""""""" |
| 61 | + |
| 62 | +matRad supports dose calculation for photons and charged particles (protons, helium, carbon ions) in external beam therapy. |
| 63 | + |
| 64 | +matRad's standard algorithms are based on the pencil-beam model, which is a fast and efficient way to compute dose distributions in a patient CT. The pencil beam model assumes that the dose at a certain voxel can be computed as the product of a depth-dependent part and a lateral part, which are computed separately. The model uses a pencil-beam kernel of an infinitely narrow beam, which is convolved with the primary fluence to compute the dose distribution. The pencil beam model is a fast and efficient way to compute dose distributions in a patient CT, but it is not as accurate as Monte Carlo simulations. |
| 65 | + |
| 66 | +The class :class:`DoseEngines.matRad_PencilBeamEngineAbstract` implements the abstract base class for all pencil beam dose engines, providing functionality for iterating through the "rays and bixels, geometry computations for distances and transformations like rotations, as well as filling of the sparse dose influence matrix. |
| 67 | + |
| 68 | +Monte Carlo dose calculation algorithms subclass from :class:`DoseEngines.matRad_MonteCarloEngineAbstract`, which provides less functionality than the base class for pencil beam algorithms, since Monte Carlo algorithms are usually integrated as third-party libraries and thus need individual implementations anyway. |
| 69 | + |
| 70 | +Ray Tracing |
| 71 | +~~~~~~~~~~~ |
| 72 | + |
| 73 | +Especially for pencil-beam algorithms, but also for other things like spot placement, matRad implements a voxel raytracing algorithm in :func:`matRad_siddonRayTracer` according to `Siddon (1985) Medical Physics <http://www.ncbi.nlm.nih.gov/pubmed/4000088>`_. Dose calculation algorithms can use :func:`matRad_rayTracing` to perform a volumetric tracing to obtain an image of the water equivalent depth from a source,according to `Siggel (2012) Physica Medica <http://www.sciencedirect.com/science/article/pii/S1120179711001359>`_. |
| 74 | + |
| 75 | +Geometric distances from the lateral distances to the central ray are computed by standard matrix vector algebra in :meth:`DoseEngines.matRad_PencilBeamEngineAbstract.calcGeoDists`. |
| 76 | + |
| 77 | +Photons |
| 78 | +~~~~~~~ |
| 79 | + |
| 80 | +Dose Calculation Algorithm |
| 81 | +########################## |
| 82 | + |
| 83 | +matRad's standard algorithm for photon dose calculation is a singular value decomposed pencil beam algorithm implemented according to `Bortfeld et al. (1993) Medical Physics <http://scitation.aip.org/content/aapm/journal/medphys/20/2/10.1118/1.597070>`_. It is implemented in :class:`DoseEngines.matRad_PhotonPencilBeamSVDEngine`. A more detailed description of the inner workings of the algorithm can also be found in the diploma thesis of Martin Siggel "Entwicklung einer schnellen Dosisberechnung auf Basis eines Pencil-Kernel Algorithmus für die Strahlentherapie mit Photonen" (Universität Heidelberg, 2008) which is available for download `here <https://raw.githubusercontent.com/wiki/e0404/matRad/resources/DiplomaMartinSiggel.pdf>`_. |
| 84 | + |
| 85 | +The dose delivered to a certain voxel *i* from bixel *j* is stored as dose influence matrix in the :ref:`dij struct <dij>` using MATLAB's built-in double precision sparse matrix format. |
| 86 | + |
| 87 | +An experimental Monte Carlo dose engine based on a the reimplementation `ompMC <https://github.com/edoerner/ompMC>`_ of the `EGSnrc <https://www.irs.inms.nrc.ca/EGSnrc/>`_ Monte Carlo code by Edgardo Doerner is available in :class:`DoseEngines.matRad_PhotonOmpMCEngine` class. It is not yet fully functional, especially the source characteristics are not fully compatible with the machine files / base data. |
| 88 | + |
| 89 | +Base data / Machine |
| 90 | +################### |
| 91 | + |
| 92 | +The necessary measured base data, namely the kernel functions as described by `Bortfeld et al. (1993) Medical Physics <http://scitation.aip.org/content/aapm/journal/medphys/20/2/10.1118/1.597070>`_ are supplied for a 6MV LINAC and stored in `photons_Generic.mat <https://github.com/e0404/matRad/blob/master/matRad/basedata/photons_Generic.mat>`_ as MATLAB piecewise polynomial. |
| 93 | + |
| 94 | +matRad's photon dose engine is calibrated such that a bixel intensity of all ones, i.e., ``w = ones(stf.totalNumOfBixels,1)``, yields a dose of roughly 1Gy in 5cm depth for a 5cm by 5cm field at SSD = 900mm. If you want to reproduce this, you can download a set of stf and pln structures that work with the TG119 phantom `here <https://raw.githubusercontent.com/wiki/e0404/matRad/resources/TG119absDosNorm.mat>`_. |
| 95 | + |
| 96 | +Approximations |
| 97 | +############## |
| 98 | + |
| 99 | +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 matrices (kernel1Mx, kernel2Mx, kernel3Mx). An extension to inhomogeneous primary fluences is easily possible. |
| 100 | + |
| 101 | +To model a virtual photon source with finite size, the primary fluence is convolved with a Gaussian filter. The size of this filter is determined by the measured penumbra at isocenter, which is related to the source size as depicted in the following images from `Martin Siggel <https://www.dlr.de/sc/desktopdefault.aspx/tabid-1192/1635_read-27193/sortby-lastname/>`_ (`Siggel M. Entwicklung einer schnellen Dosisberechnung auf Basis eines Pencil-Kernel Algorithmus für die Strahlentherapie mit Photonen. Diploma Thesis 2008, Universität Heidelberg <https://raw.githubusercontent.com/wiki/e0404/matRad/resources/DiplomaMartinSiggel.pdf>`_): |
| 102 | + |
| 103 | +.. image:: /images/source_penumbra.png |
| 104 | +.. image:: /images/penumbra_gaussian_convolution.png |
| 105 | + |
| 106 | +matRad initially used a hardcoded value of 5mm for the measured penumbra width. From version 2.10.1 on, the value can also be stored in the photon base data file so that it can also be translated to the virtual source size for, for example, Monte Carlo simulations. |
| 107 | + |
| 108 | +Computational bottlenecks |
| 109 | +######################### |
| 110 | + |
| 111 | +For the photon dose calculation, there are three main possibilities to speed up computations and reduce memory consumption. |
| 112 | + |
| 113 | +1. 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. |
| 114 | +2. You can increase the variable ``pln.bixelWidth`` in the `matRad script <https://github.com/e0404/matRad/blob/master/matRad.m>`_ which effectively reduces the number of bixels which have to be computed approximately quadratically. |
| 115 | +3. You can reduce the radius around the central ray, around which dose is computed by adjusting the variable ``lateralCutOff`` in ``pln.propDoseCalc``. Currently this variable is already set to 50mm. |
| 116 | + |
| 117 | +Particles |
| 118 | +~~~~~~~~~ |
| 119 | + |
| 120 | +Dose calculation algorithm |
| 121 | +########################## |
| 122 | + |
| 123 | +matRad's default conventional pencil beam model for particle dose calculations is implemented in :class:`DoseEngines.matRad_ParticleHongPencilBeamEngine` and similar to the work of `Hong et al. (1996) <http://iopscience.iop.org/0031-9155/41/8/005/>`_. The dose at a particular 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 broadening, matRad uses a depth-dependent sigma of a Gaussian profile, which is also tabulated versus depth for all available beam energies. |
| 124 | + |
| 125 | +The dose delivered to a certain voxel *i* from bixel *j* is stored as dose influence matrix in the `dij struct <dij>`_ using MATLAB's built-in double precision sparse matrix format. |
| 126 | + |
| 127 | +α- and β-matrix pre-computations |
| 128 | +################################ |
| 129 | + |
| 130 | +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. |
| 131 | + |
| 132 | +matRad only models variations of α and β with depth in `matRad_calcLQParameter <https://github.com/e0404/matRad/blob/master/matRad_calcLQParameter.m>`_. Potential dependencies in lateral direction are not explicitly modeled. |
| 133 | + |
| 134 | +Base data |
| 135 | +######### |
| 136 | + |
| 137 | +The base data files `protons_Generic <https://github.com/e0404/matRad/blob/master/matRad/basedata/protons_Generic.mat>`_ and `carbon_Generic <https://github.com/e0404/matRad/blob/master/matRad/basedata/carbon_Generic.mat>`_ 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 (`Bortfeld (1997) <http://www.ncbi.nlm.nih.gov/pubmed/9434986>`_) and Highland's approximation for multiple Coulomb scattering (`Gottschalk (1993) <http://www.sciencedirect.com/science/article/pii/0168583X9395944Z>`_). 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 <https://github.com/e0404/matRad/blob/master/matRad_calcParticleDoseBixel.m>`_. Given that the lateral components of the dose calculation have the unit [1/mm^2], the weight of the pencil beams in matRad directly corresponds to the number of particles in [1e6] while the dose is given in [Gy]. |
| 138 | + |
| 139 | +More detailed information on the structure of a particle base data file is given :ref:`here <basedata_particles>`! |
| 140 | + |
| 141 | +Approximations |
| 142 | +############## |
| 143 | + |
| 144 | +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. |
| 145 | + |
| 146 | +Computational bottlenecks |
| 147 | +######################### |
| 148 | + |
| 149 | +For the particle dose calculation, there are three main possibilities to speed up computations and reduce memory consumption. |
| 150 | + |
| 151 | +1. 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. |
| 152 | +2. You can increase the variable ``pln.propStf.bixelWidth`` in the :scpt:`matRad script <matRad.m>`_ which effectively reduces the number of pencil beams which have to be computed approximately quadratically. |
| 153 | +3. You can reduce the radius around the central ray where dose is computed by adjusting the :attr:`DoseEngines.matRad_PencilBeamEngineAbstract.dosimetricLateralCutOff` property. |
0 commit comments