Implementation of LET modifying objectives

AUTHORS.txt

@@ -30,6 +30,7 @@
 * Giuseppe Pezzano
 * Daniel Ramirez
 * Carsten Scholz
+* Lisa Seckler
 * Camilo Sevilla
 * Alexander Stadler
 * Uwe Titt

examples/matRad_example12_simpleParticleMonteCarlo.m

@@ -108,7 +108,7 @@
 
 %% Compare LET
 if isfield(resultGUI,'LET') && isfield(resultGUI_MC,'LET')
-    matRad_compareDose(resultGUI.LET, resultGUI.(['LET_' pln.propDoseCalc.engine]), ct, cst, [1, 1, 0] , 'off', pln, [2, 2], 1, 'global');    
+    matRad_compareDose(resultGUI.LETd, resultGUI.(['LET_' pln.propDoseCalc.engine]), ct, cst, [1, 1, 0] , 'off', pln, [2, 2], 1, 'global');    
 end
 
 %% GUI

examples/matRad_example16_LETd.m

@@ -0,0 +1,170 @@
+%% Welcome to a LETd example script
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright 2017 the matRad development team. 
+% 
+% This file is part of the matRad project. It is subject to the license 
+% terms in the LICENSE file found in the top-level directory of this 
+% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 
+% of the matRad project, including this file, may be copied, modified, 
+% propagated, or distributed except according to the terms contained in the 
+% LICENSE file.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% open matRad files and delete everything in the workspace
+matRad_rc
+
+%% load an easy phantom like the TG119.mat
+load("TG119.mat")
+
+%% Plan and Geometry
+pln.radiationMode   = 'protons';           % either photons / protons / helium / carbon
+pln.machine         = 'Generic';
+pln.numOfFractions  = 30;
+
+
+% beam geometry settings
+pln.propStf.bixelWidth      = 5; % [mm] / also corresponds to lateral spot spacing for particles
+pln.propStf.gantryAngles    = [45 0 -45];
+pln.propStf.couchAngles     = zeros(numel(pln.propStf.gantryAngles),1);  
+pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
+pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
+% optimization settings
+pln.propDoseCalc.calcLET = 1;   % very important, don't forget that one!
+
+pln.propOpt.runDAO          = false;      % 1/true: run DAO, 0/false: don't / will be ignored for particles
+pln.propOpt.runSequencing   = false;      % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below
+pln.propOpt.spatioTemp      = 0;
+
+% dose calculation settings
+pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
+% pln(1).propDoseCalc.doseGrid.resolution = ct.resolution;
+quantityOpt  = 'RBExD';     % options: physicalDose, effect, RBExD
+%=======================================> Model check error in bioModel
+modelName    = 'constRBE';             % none: for photons, protons, carbon            % constRBE: constant RBE for photons and protons 
+                                   % MCN: McNamara-variable RBE model for protons  % WED: Wedenberg-variable RBE model for protons 
+                                   % LEM: Local Effect Model for carbon ions
+
+
+scenGenType  = 'nomScen';          % scenario creation type 'nomScen'  'wcScen' 'impScen' 'rndScen'                                          
+
+% retrieve bio model parameters
+pln.bioParam = matRad_bioModel(pln.radiationMode,quantityOpt, modelName);
+
+% retrieve scenarios for dose calculation and optimziation
+pln.multScen = matRad_multScen(ct,scenGenType);
+stf = matRad_generateStf(ct,cst,pln);
+
+%% Dose calculation
+% Dij Calculation --> only for dose
+dij = matRad_calcParticleDose(ct,stf,pln,cst);
+
+% Dirty Dose Calculation --> compute dirty dose influence matrix
+%  arg: LET threshold (keV/um)
+dij = matRad_calcDirtyDose(2,dij);
+
+%% Optimization
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+
+%% Plotting
+% Let's see how it looks
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETd;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETd without LETd objective')
+zoom(4)
+
+%% Adding LETd objective
+% adding a LETd objective to the Core
+cst{1,6}{2} = struct(LETdObjectives.matRad_SquaredOverdosingLETd(100,0));
+
+%% Optimization
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+
+%% Well done your first LETd calculation is ready!
+% Let's see how it looks
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETd;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETd with LETd objective in Core')
+zoom(4)
+
+%% Now with the Target
+%% Adding LETd
+% for adding a LETd objective to OuterTarget
+cst{2,6}{2} = struct(LETdObjectives.matRad_SquaredUnderdosingLETd(100,20));
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+%% Plotting
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETd;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETd with LETd objective in Target')
+zoom(4)
+
+%% Now with the Core and Target
+%% Adding LETd
+% for adding a LETd objective choose two structures: Core and OuterTarget
+cst{1,6}{2} = struct(LETdObjectives.matRad_SquaredOverdosingLETd(100,0));
+cst{2,6}{2} = struct(LETdObjectives.matRad_SquaredUnderdosingLETd(100,20));
+%% Optimization
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+%% Plotting
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETd;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETd with LETd objective in Core and Target')
+zoom(4)

examples/matRad_example17_LETxDose.m

@@ -0,0 +1,170 @@
+%% Welcome to a LETxDose example script
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright 2017 the matRad development team. 
+% 
+% This file is part of the matRad project. It is subject to the license 
+% terms in the LICENSE file found in the top-level directory of this 
+% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 
+% of the matRad project, including this file, may be copied, modified, 
+% propagated, or distributed except according to the terms contained in the 
+% LICENSE file.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% open matRad files and delete everything in the workspace
+matRad_rc
+
+%% load an easy phantom like the TG119.mat
+load("TG119.mat")
+
+%% Plan and Geometry
+pln.radiationMode   = 'protons';           % either photons / protons / helium / carbon
+pln.machine         = 'Generic';
+pln.numOfFractions  = 30;
+
+
+% beam geometry settings
+pln.propStf.bixelWidth      = 5; % [mm] / also corresponds to lateral spot spacing for particles
+pln.propStf.gantryAngles    = [45 0 -45];
+pln.propStf.couchAngles     = zeros(numel(pln.propStf.gantryAngles),1);  
+pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
+pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
+% optimization settings
+pln.propDoseCalc.calcLET = 1;   % very important, don't forget that one!
+
+pln.propOpt.runDAO          = false;      % 1/true: run DAO, 0/false: don't / will be ignored for particles
+pln.propOpt.runSequencing   = false;      % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below
+pln.propOpt.spatioTemp      = 0;
+
+% dose calculation settings
+pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
+% pln(1).propDoseCalc.doseGrid.resolution = ct.resolution;
+quantityOpt  = 'RBExD';     % options: physicalDose, effect, RBExD
+%=======================================> Model check error in bioModel
+modelName    = 'constRBE';             % none: for photons, protons, carbon            % constRBE: constant RBE for photons and protons 
+                                   % MCN: McNamara-variable RBE model for protons  % WED: Wedenberg-variable RBE model for protons 
+                                   % LEM: Local Effect Model for carbon ions
+
+
+scenGenType  = 'nomScen';          % scenario creation type 'nomScen'  'wcScen' 'impScen' 'rndScen'                                          
+
+% retrieve bio model parameters
+pln.bioParam = matRad_bioModel(pln.radiationMode,quantityOpt, modelName);
+
+% retrieve scenarios for dose calculation and optimziation
+pln.multScen = matRad_multScen(ct,scenGenType);
+stf = matRad_generateStf(ct,cst,pln);
+
+%% Dose calculation
+% Dij Calculation --> only for dose
+dij = matRad_calcParticleDose(ct,stf,pln,cst);
+
+% Dirty Dose Calculation --> compute dirty dose influence matrix
+%  arg: LET threshold (keV/um)
+dij = matRad_calcDirtyDose(2,dij);
+
+%% Optimization
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+
+%% Plotting
+% Let's see how it looks
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETxDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETxDose without LETxDose objective')
+zoom(4)
+
+%% Adding LETxDose objective
+% adding a LETxDose objective to the Core
+cst{1,6}{2} = struct(LETxDoseObjectives.matRad_SquaredOverdosingLETxDose(6,0));
+
+%% Optimization
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+
+%% Well done your first LETxDose calculation is ready!
+% Let's see how it looks
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETxDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETxDose with LETxDose objective in Core')
+zoom(4)
+
+%% Now with the Target
+%% Adding LETxDose
+% for adding a LETxDose objective to OuterTarget
+cst{2,6}{2} = struct(LETxDoseObjectives.matRad_SquaredUnderdosingLETxDose(6,20));
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+%% Plotting
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETxDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETxDose with LETxDose objective in Target')
+zoom(4)
+
+%% Now with the Core and Target
+%% Adding LETxDose
+% for adding a LETxDose objective choose two structures: Core and OuterTarget
+cst{1,6}{2} = struct(LETxDoseObjectives.matRad_SquaredOverdosingLETxDose(6,0));
+cst{2,6}{2} = struct(LETxDoseObjectives.matRad_SquaredUnderdosingLETxDose(6,20));
+%% Optimization
+resultGUI = matRad_fluenceOptimization(dij,cst,pln);
+%% Plotting
+cube = resultGUI.physicalDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+figure
+subplot(2,1,1)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('physicalDose')
+zoom(4)
+
+cube = resultGUI.LETxDose;
+plane = 3;
+slice = 80;
+doseWindow = [min(cube(:)) max(cube(:))];
+subplot(2,1,2)
+matRad_plotSliceWrapper(gca,ct,cst,1,cube,plane,slice,[],[],colorcube,[],doseWindow,[],[]);
+title('LETxDose with LETxDose objective in Core and Target')
+zoom(4)