diff --git a/AUTHORS.txt b/AUTHORS.txt index 603501ba8..35ef7dd14 100644 --- a/AUTHORS.txt +++ b/AUTHORS.txt @@ -30,6 +30,7 @@ * Giuseppe Pezzano * Daniel Ramirez * Carsten Scholz +* Lisa Seckler * Camilo Sevilla * Alexander Stadler * Uwe Titt diff --git a/examples/matRad_example12_simpleParticleMonteCarlo.m b/examples/matRad_example12_simpleParticleMonteCarlo.m index 1f9e207f4..be8072876 100644 --- a/examples/matRad_example12_simpleParticleMonteCarlo.m +++ b/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 diff --git a/examples/matRad_example16_LETd.m b/examples/matRad_example16_LETd.m new file mode 100644 index 000000000..482f167de --- /dev/null +++ b/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) diff --git a/examples/matRad_example17_LETxDose.m b/examples/matRad_example17_LETxDose.m new file mode 100644 index 000000000..1787d5997 --- /dev/null +++ b/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) diff --git a/examples/matRad_example18_DirtyDose.m b/examples/matRad_example18_DirtyDose.m new file mode 100644 index 000000000..45fc0621a --- /dev/null +++ b/examples/matRad_example18_DirtyDose.m @@ -0,0 +1,132 @@ +%% Welcome to a dirty dose 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.dirtyDose; +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('dirtyDose without DD objective') +zoom(4) + +%% Adding dirty dose objective +% adding a dirty dose objective to the Core +cst{1,6}{2} = struct(DirtyDoseObjectives.matRad_SquaredOverdosingDirtyDose(300,0)); + +%% Optimization +resultGUI = matRad_fluenceOptimization(dij,cst,pln); + +%% Well done your first dirty dose 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.dirtyDose; +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('dirtyDose with DD objective in Core') +zoom(4) + +%% Now with the Target +%% Adding dirty dose +% for adding a dirty dose objective to OuterTarget +cst{2,6}{2} = struct(DirtyDoseObjectives.matRad_SquaredUnderdosingDirtyDose(300,20)); +resultGUI = matRad_fluenceOptimization(dij,cst,pln); + +%% Now with the Core and Target +%% Adding dirty dose +% for adding a dirty dose objective choose two structures: Core and OuterTarget +cst{1,6}{2} = struct(DirtyDoseObjectives.matRad_SquaredOverdosingDirtyDose(300,0)); +cst{2,6}{2} = struct(DirtyDoseObjectives.matRad_SquaredUnderdosingDirtyDose(300,20)); +%% Optimization +resultGUI = matRad_fluenceOptimization(dij,cst,pln); diff --git a/examples/matRad_example7_carbon.m b/examples/matRad_example7_carbon.m index 0b3756881..256a1660a 100644 --- a/examples/matRad_example7_carbon.m +++ b/examples/matRad_example7_carbon.m @@ -118,7 +118,7 @@ % Let's plot the transversal iso-center LET slice slice = round(pln.propStf.isoCenter(3)./ct.resolution.z); figure; -imagesc(resultGUI.LET(:,:,slice)),colorbar, colormap(jet); +imagesc(resultGUI.LETd(:,:,slice)),colorbar, colormap(jet); %% Inverse Optimization for IMPT based on biological effect % To perform a dose optimization for carbon ions we can also use the diff --git a/matRad/doseCalc/matRad_calcDirtyDose.m b/matRad/doseCalc/matRad_calcDirtyDose.m new file mode 100644 index 000000000..fb9160575 --- /dev/null +++ b/matRad/doseCalc/matRad_calcDirtyDose.m @@ -0,0 +1,22 @@ +function dij = matRad_calcDirtyDose(LET_thres,dij) +% Calculates Dirty and Clean Dose by using LET threshold +% +% call +% dij = matRad_calcDirtyDose(LET_thres,dij) +% +% input +% LET_thres: LET threshold, above: dirty dose, below: clean dose +% dij: matRad dij struct +% +% output +% dij: matRad dij struct with dirty and clean dose +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +dij.dirtyDoseThreshold = LET_thres; +[dij.LETmaskDirty,dij.LETmaskClean,dij.mLET] = matRad_calcLETmask(dij); + +dij.dirtyDose = cellfun(@times,dij.LETmaskDirty, dij.physicalDose,'UniformOutput',false); +dij.cleanDose = cellfun(@times,dij.LETmaskClean, dij.physicalDose,'UniformOutput',false); + +end \ No newline at end of file diff --git a/matRad/doseCalc/matRad_calcLETmask.m b/matRad/doseCalc/matRad_calcLETmask.m new file mode 100644 index 000000000..5cdf84d24 --- /dev/null +++ b/matRad/doseCalc/matRad_calcLETmask.m @@ -0,0 +1,52 @@ +function [LETmaskDirty,LETmaskClean,mLET] = matRad_calcLETmask(dij) +% Calculates logical matrix where LET is above (LETmaskDirty) or below +% (LETmaskClean) a certain threshold which is set in matRad_calcDirtyDose +% +% call +% [LETmaskDirty,LETmaskClean,mLET] = matRad_calcLETmask(dij) +% +% input +% dij: matRad dij struct after the dose calculation using +% matRad_calcParticleDose +% +% output +% LETmaskDirty: logical matrix (0 and 1) for dirty dose +% -> it needs to go into matRad_calcDirtyDose to create +% dirty dose as a variable +% LETmaskClean: logical matrix (0 and 1) for clean dose -> into matRad_calcDirtyDose +% mLET: mLETDose matrix divided by the physicalDose matrix +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +try %Leverage similar Nonzero Structure + for scenIx = find(~cellfun(@isempty,dij.physicalDose)) + [i,j] = find(dij.mLETDose{scenIx}); + mLET{scenIx} = sparse(i,j,nonzeros(dij.mLETDose{scenIx}) ./ nonzeros(dij.physicalDose{scenIx}),dij.doseGrid.numOfVoxels,dij.totalNumOfBixels); + + LETmaskDirty{scenIx} = mLET{scenIx} > dij.dirtyDoseThreshold; + subIxClean = nonzeros(mLET{scenIx}) < dij.dirtyDoseThreshold; + + LETmaskClean{scenIx} = sparse(i(subIxClean),j(subIxClean),true,dij.doseGrid.numOfVoxels,dij.totalNumOfBixels); + + %{ + [i,j,v] = find(dij.physicalDose{1}); + idx = sub2ind(size(dij.physicalDose{1}),i,j); + + mLET{1} = full(dij.mLETDose{1}(idx) ./ v); + + subIxDirty = mLET{1} > dij.dirtyDoseThreshold; + subIxClean = ~subIxDirty; + + mLET{1} = sparse(i,j, mLET{1},dij.doseGrid.numOfVoxels,dij.totalNumOfBixels); + + LETmaskDirty{1} = sparse(i(subIxDirty),j(subIxDirty),true,dij.doseGrid.numOfVoxels,dij.totalNumOfBixels); + LETmaskClean{1} = sparse(i(subIxClean),j(subIxClean),true,dij.doseGrid.numOfVoxels,dij.totalNumOfBixels); + %} + end +catch ME + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Nonzero structure of LETxDose and physicalDose matrix not similar or matrices do not exist!'); +end + + +end \ No newline at end of file diff --git a/matRad/gui/widgets/matRad_OptimizationWidget.m b/matRad/gui/widgets/matRad_OptimizationWidget.m index af4bf5340..678d2499b 100644 --- a/matRad/gui/widgets/matRad_OptimizationWidget.m +++ b/matRad/gui/widgets/matRad_OptimizationWidget.m @@ -389,6 +389,45 @@ 'TooltipString','Objective Penalty', ... 'UserData',[i,j,0],... 'Callback',@(hObject,eventdata)editObjParam_Callback(this,hObject,eventdata)); + elseif isa(obj,'DirtyDoseObjectives.matRad_DirtyDoseObjective') + h = uicontrol(cstPanel,'Style','edit', ... + 'String',num2str(obj.penalty), ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) penaltyW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Tooltip','Objective Penalty', ... + 'UserData',[i,j,0],... + 'Callback',@(hObject,eventdata)editObjParam_Callback(this,hObject,eventdata)); + elseif isa(obj,'LETxDoseObjectives.matRad_LETxDoseObjective') + h = uicontrol(cstPanel,'Style','edit', ... + 'String',num2str(obj.penalty), ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) penaltyW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Tooltip','Objective Penalty', ... + 'UserData',[i,j,0],... + 'Callback',@(hObject,eventdata)editObjParam_Callback(this,hObject,eventdata)); + elseif isa(obj,'LETdObjectives.matRad_LETdObjective') + h = uicontrol(cstPanel,'Style','edit', ... + 'String',num2str(obj.penalty), ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) penaltyW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Tooltip','Objective Penalty', ... + 'UserData',[i,j,0],... + 'Callback',@(hObject,eventdata)editObjParam_Callback(this,hObject,eventdata)); else h = uicontrol(cstPanel,'Style','edit','String','----','Units','normalized','Position',[xPos ypos(cnt) penaltyW objHeight], 'FontSize',matRad_cfg.gui.fontSize,'FontName',matRad_cfg.gui.fontName,'FontWeight',matRad_cfg.gui.fontWeight,'BackgroundColor',matRad_cfg.gui.elementColor,'ForegroundColor',matRad_cfg.gui.textColor,'Enable','off'); end @@ -582,6 +621,11 @@ function changeObjFunction_Callback(this,hObject, ~) %if (isfield(currentObj,'penalty') || isobject (currentObj ) && isprop(currentObj,'penalty')) && isprop(newObj,'penalty') if (isfield(currentObj,'penalty') || isa(currentObj,'DoseObjectives.matRad_DoseObjective')) && isa(newObj,'DoseObjectives.matRad_DoseObjective') + elseif (isfield(currentObj,'penalty') || isa(currentObj,'DirtyDoseObjectives.matRad_DirtyDoseObjective')) && isa(newObj,'DirtyDoseObjectives.matRad_DirtyDoseObjective') + newObj.penalty = currentObj.penalty; + elseif (isfield(currentObj,'penalty') || isa(currentObj,'LETxDoseObjectives.matRad_LETxDoseObjective')) && isa(newObj,'LETxDoseObjectives.matRad_LETxDoseObjective') + newObj.penalty = currentObj.penalty; + elseif (isfield(currentObj,'penalty') || isa(currentObj,'LETdObjectives.matRad_LETdObjective')) && isa(newObj,'LETdObjectives.matRad_LETdObjective') newObj.penalty = currentObj.penalty; end diff --git a/matRad/matRad_fluenceOptimization.m b/matRad/matRad_fluenceOptimization.m index e5f9951d9..7b42262e0 100644 --- a/matRad/matRad_fluenceOptimization.m +++ b/matRad/matRad_fluenceOptimization.m @@ -189,6 +189,7 @@ maxCurrRBE = max(-cst{ixTarget,5}.alphaX + sqrt(cst{ixTarget,5}.alphaX^2 + ... 4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*doseTmp(V))); wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(doseTmp(V))))* wOnes; + end matRad_cfg.dispInfo('chosen weights adapted to biological dose calculation!\n'); @@ -260,7 +261,7 @@ case 'physicalDose' backProjection = matRad_DoseProjection; otherwise - warning(['Did not recognize bioloigcal setting ''' pln.probOpt.bioOptimization '''!\nUsing physical dose optimization!']); + matRad_cfg.dispWarning('Did not recognize biological setting ''%s''!\nUsing physical dose optimization!',pln.probOpt.bioOptimization); backProjection = matRad_DoseProjection; end @@ -276,6 +277,22 @@ optiProb.useLogSumExpForRobOpt = pln.propOpt.useLogSumExpForRobOpt; end +rowCst = size(cst); +for i = 1 : rowCst(1,1) + num = size(cst{i,6}); + for j = 1 : num(1,2) + if isa(cst{i,6}{j},'DirtyDoseObjectives.matRad_DirtyDoseObjective') + optiProb.dirtyDoseBP = matRad_DirtyDoseProjection; + end + if isa(cst{i,6}{j},'LETxDoseObjectives.matRad_LETxDoseObjective') + optiProb.LETxDoseBP = matRad_LETxDoseProjection; + end + if isa(cst{i,6}{j},'LETdObjectives.matRad_LETdObjective') + optiProb.LETdBP = matRad_LETdProjection; + end + end +end + %Get Bounds if ~isfield(pln.propOpt,'boundMU') pln.propOpt.boundMU = false; @@ -326,6 +343,7 @@ resultGUI = matRad_calcCubes(wOpt,dij); resultGUI.wUnsequenced = wOpt; + resultGUI.usedOptimizer = optimizer; resultGUI.info = info; diff --git a/matRad/optimization/+DirtyDoseObjectives/matRad_DirtyDoseObjective.m b/matRad/optimization/+DirtyDoseObjectives/matRad_DirtyDoseObjective.m new file mode 100644 index 000000000..545d84cb0 --- /dev/null +++ b/matRad/optimization/+DirtyDoseObjectives/matRad_DirtyDoseObjective.m @@ -0,0 +1,68 @@ +classdef (Abstract) matRad_DirtyDoseObjective < matRad_DoseOptimizationFunction +% matRad_DirtyDoseObjective: Interface for optimization objectives +% This abstract base class provides the structure of optimization +% objectives like mean dirty dose, squared deviation dirty dose, EUD dirty dose, dose-volume etc. +% Implementations can be found in the DirtyDoseObjectives package +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Abstract, Access = public) + penalty %Optimization penalty + + end + + methods (Static) + function rob = availableRobustness() + rob = {'none','STOCH','PROB','VWWC','VWWC_INV','COWC','OWC'}; %By default, no robustness is available + end + end + + %These should be abstract methods, however Octave can't parse them. As soon + %as Octave is able to do this, they should be made abstract again + methods %(Abstract) + + %returns the objective function value for the given dose vector. Needs to be implemented in sub-classes. + function fDirtyDose = computeDirtyDoseObjectiveFunction(obj,dirtyDose) + error('Function needs to be implemented!'); + end + + + %returns the dose-gradient for the given dose vector. Needs to be implemented in sub-classes. + function fDirtyDoseGrad = computeDirtyDoseObjectiveGradient(obj,dirtyDose) + error('Function needs to be implemented!'); + end + + end + + methods (Access = public) + + % constructor of matRad_DoseObjective + function obj = matRad_DirtyDoseObjective(varargin) + %default initialization from struct (parameters & penalty) + obj@matRad_DoseOptimizationFunction(varargin{:}); + end + + %Overloads the struct function to add Objective related information + %to output struct + function s = struct(obj) + s = struct@matRad_DoseOptimizationFunction(obj); + s.penalty = obj.penalty; + + end + + end +end + diff --git a/matRad/optimization/+DirtyDoseObjectives/matRad_DirtyDoseVariance.m b/matRad/optimization/+DirtyDoseObjectives/matRad_DirtyDoseVariance.m new file mode 100644 index 000000000..250af1af7 --- /dev/null +++ b/matRad/optimization/+DirtyDoseObjectives/matRad_DirtyDoseVariance.m @@ -0,0 +1,74 @@ +classdef matRad_DirtyDoseVariance < DirtyDoseObjectives.matRad_DirtyDoseObjective +% matRad_DirtyDoseVariance implements a variance objective for +% homogenous dirty dose distribution +% See matRad_DoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Variance'; + parameterNames = {}; + parameterTypes = {}; + end + + properties + parameters = {}; + penalty = 1; + end + + methods + function obj = matRad_DirtyDoseVariance(penalty) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@DirtyDoseObjectives.matRad_DirtyDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(dRef) + obj.parameters{1} = dRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fDirtyDose = computeDirtyDoseObjectiveFunction(obj,dirtyDose) + fDirtyDose = obj.penalty * var(dirtyDose); + + end + + %% Calculates the Objective Function gradient + function fDirtyDoseGrad = computeDirtyDoseObjectiveGradient(obj,dirtyDose) + fDirtyDoseGrad = obj.penalty * 2/(numel(dirtyDose) - 1) * (dirtyDose - mean(dirtyDose)); + end + end + +end \ No newline at end of file diff --git a/matRad/optimization/+DirtyDoseObjectives/matRad_MeanDirtyDose.m b/matRad/optimization/+DirtyDoseObjectives/matRad_MeanDirtyDose.m new file mode 100644 index 000000000..8b8075d7f --- /dev/null +++ b/matRad/optimization/+DirtyDoseObjectives/matRad_MeanDirtyDose.m @@ -0,0 +1,130 @@ +classdef matRad_MeanDirtyDose < DirtyDoseObjectives.matRad_DirtyDoseObjective +% matRad_MeanDirtyDose Implements a penalized MeanDirtyDose objective +% See matRad_DirtyDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Mean Dirty Dose'; + parameterNames = {'d^{ref}','f_{diff}'}; %When optimizing to a reference, one might consider using a quadratic relationship with a non-linear optimizer + parameterTypes = {'dirtyDose',{'Linear','Quadratic'}}; + end + + properties + parameters = {0,1}; + penalty = 1; + end + + methods + function obj = matRad_MeanDirtyDose(penalty,dMeanRef,fDiff) + + % if we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@DirtyDoseObjectives.matRad_DirtyDoseObjective(inputStruct); + + if ~initFromStruct + if nargin < 3 || ~ischar(fDiff) + fDiff = 'Linear'; + end + + fDiffIx = find(strcmp(fDiff,obj.parameterTypes{2})); + + if isempty(fDiffIx) || numel(fDiffIx) > 1 + fDiffIx = 1; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Mean dirty dose difference function can only be %s! Using %s difference.', strjoin(obj.parameterTypes{2},' or '), obj.parameterTypes{2}{fDiffIx}); + end + + obj.parameters{2} = fDiffIx; + + + if nargin >= 2 && isscalar(dMeanRef) + obj.parameters{1} = dMeanRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + + %% Downwards compatability / set default values + %TODO: maybe move into set method for parameters + if numel(obj.parameters) < 1 + obj.parameters{1} = 0; + end + + if numel(obj.parameters) < 2 + obj.parameters{2} = 1; + end + + end + + %% Calculates the Objective Function value + function fDirtyDose = computeDirtyDoseObjectiveFunction(obj,dirtyDose) + switch obj.parameters{2} + case 1 + fDirtyDose = obj.objectiveLinearDiff(dirtyDose); + case 2 + fDirtyDose = obj.objectiveQuadraticDiff(dirtyDose); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid setting for %s in Mean Dirty Dose Objective!',obj.parameterNames{2}); + end + end + + %% Calculates the Objective Function gradient + function fDirtyDoseGrad = computeDirtyDoseObjectiveGradient(obj,dirtyDose) + switch obj.parameters{2} + case 1 + fDirtyDoseGrad = obj.gradientLinearDiff(dirtyDose); + case 2 + fDirtyDoseGrad = obj.gradientQuadraticDiff(dirtyDose); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid setting for %s in Mean Dirty Dose Objective!',obj.parameterNames{2}); + end + end + end + + methods (Access = protected) + function fDirtyDose = objectiveQuadraticDiff(obj,dirtyDose) + fDirtyDose = (mean(dirtyDose(:)) - obj.parameters{1})^2; + end + + function fDirtyDoseGrad = gradientQuadraticDiff(obj,dirtyDose) + fDirtyDoseGrad = 2*(mean(dirtyDose(:))-obj.parameters{1}) * ones(size(dirtyDose(:)))/numel(dirtyDose); + end + + function fDirtyDose = objectiveLinearDiff(obj,dirtyDose) + fDirtyDose = abs(mean(dirtyDose(:)) - obj.parameters{1}); + end + + function fDirtyDoseGrad = gradientLinearDiff(obj,dirtyDose) + fDirtyDoseGrad = (1/numel(dirtyDose))*sign(dirtyDose(:)-obj.parameters{1}); + end + end + +end + diff --git a/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredDeviationDirtyDose.m b/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredDeviationDirtyDose.m new file mode 100644 index 000000000..582be1eca --- /dev/null +++ b/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredDeviationDirtyDose.m @@ -0,0 +1,80 @@ +classdef matRad_SquaredDeviationDirtyDose < DirtyDoseObjectives.matRad_DirtyDoseObjective +% matRad_SquaredDeviationDirtyDose Implements a penalized least squares dirtyDose objective +% See matRad_DirtyDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Deviation Dirty Dose'; + parameterNames = {'d^{ref}'}; + parameterTypes = {'dirtyDose'}; + end + + properties + parameters = {60}; + penalty = 1; + end + + + methods + function obj = matRad_SquaredDeviationDirtyDose(penalty,dRef) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@DirtyDoseObjectives.matRad_DirtyDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(dRef) + obj.parameters{1} = dRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fDirtyDose = computeDirtyDoseObjectiveFunction(obj,dirtyDose) + % deviation : dirtyDose minus prefered dose + deviation = dirtyDose - obj.parameters{1}; + + % calculate objective function + fDirtyDose = 1/numel(dirtyDose) * (deviation'*deviation); + end + + %% Calculates the Objective Function gradient + function fDirtyDoseGrad = computeDirtyDoseObjectiveGradient(obj,dirtyDose) + % deviation : dirtyDose minus prefered dose + deviation = dirtyDose - obj.parameters{1}; + + % calculate delta + fDirtyDoseGrad = 2 * 1/numel(dirtyDose) * deviation; + end + end +end \ No newline at end of file diff --git a/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredOverdosingDirtyDose.m b/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredOverdosingDirtyDose.m new file mode 100644 index 000000000..520e3b18a --- /dev/null +++ b/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredOverdosingDirtyDose.m @@ -0,0 +1,83 @@ +classdef matRad_SquaredOverdosingDirtyDose < DirtyDoseObjectives.matRad_DirtyDoseObjective +% matRad_SquaredOverdosingDirtyDose implements a penalized dirty dose +% See matRad_DirtyDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Overdosing Dirty Dose'; + parameterNames = {'d^{max}'}; + parameterTypes = {'dirtyDose'}; + end + + properties + parameters = {30}; + penalty = 1; + end + + methods + function obj = matRad_SquaredOverdosingDirtyDose(penalty,dMax) + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@DirtyDoseObjectives.matRad_DirtyDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin >= 2 && isscalar(dMax) + obj.parameters{1} = dMax; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fDirtyDose = computeDirtyDoseObjectiveFunction(obj,dirtyDose) + % overdose : dirtyDose minus prefered dose + overdose = dirtyDose - obj.parameters{1}; + + % apply positive operator + overdose(overdose<0) = 0; + + % calculate objective function + fDirtyDose = 1/numel(dirtyDose) * (overdose'*overdose); + end + + %% Calculates the Objective Function gradient + function fDirtyDoseGrad = computeDirtyDoseObjectiveGradient(obj,dirtyDose) + % overdose : dirtyDose minus prefered dose + overdose = dirtyDose - obj.parameters{1}; + + % apply positive operator + overdose(overdose<0) = 0; + + % calculate delta + fDirtyDoseGrad = 2 * 1/numel(dirtyDose) * overdose; + end + end + +end diff --git a/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredUnderdosingDirtyDose.m b/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredUnderdosingDirtyDose.m new file mode 100644 index 000000000..220d4566a --- /dev/null +++ b/matRad/optimization/+DirtyDoseObjectives/matRad_SquaredUnderdosingDirtyDose.m @@ -0,0 +1,83 @@ +classdef matRad_SquaredUnderdosingDirtyDose < DirtyDoseObjectives.matRad_DirtyDoseObjective +% matRad_SquaredUnderdosingDirtyDose implements a penalized dirty dose +% See matRad_DirtyDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Underdosing Dirty Dose'; + parameterNames = {'d^{min}'}; + parameterTypes = {'dirtyDose'}; + end + + properties + parameters = {60}; + penalty = 1; + end + + methods + function obj = matRad_SquaredUnderdosingDirtyDose(penalty,dMin) + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@DirtyDoseObjectives.matRad_DirtyDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin >= 2 && isscalar(dMin) + obj.parameters{1} = dMin; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fDirtyDose = computeDirtyDoseObjectiveFunction(obj,dirtyDose) + % underdose : dirtyDose minus prefered dose + underdose = dirtyDose - obj.parameters{1}; + + % apply positive operator + underdose(underdose>0) = 0; + + % calculate objective function + fDirtyDose = 1/numel(dirtyDose) * (underdose'*underdose); + end + + %% Calculates the Objective Function gradient + function fDirtyDoseGrad = computeDirtyDoseObjectiveGradient(obj,dirtyDose) + % underdose : dirtyDose minus prefered dose + underdose = dirtyDose - obj.parameters{1}; + + % apply positive operator + underdose(underdose>0) = 0; + + % calculate delta + fDirtyDoseGrad = 2/numel(dirtyDose) * underdose; + end + end + +end diff --git a/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m b/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m index 5322c97b6..f711fcade 100644 --- a/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m +++ b/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m @@ -21,6 +21,7 @@ properties (Abstract, Access = public) penalty %Optimization penalty + end methods (Static) @@ -43,6 +44,13 @@ function fDoseGrad = computeDoseObjectiveGradient(obj,dose) error('Function needs to be implemented!'); end + + % function fClusterDose = computeClusterDoseObjectiveFunction(obj,dose) + % error('Function needs to be implemented!') + % end + % function fClusterDoseGrad = computeClusterDoseObjectiveGradient(obj,dose) + % error('Function needs to be implemented!') + %end end @@ -59,6 +67,7 @@ function s = struct(obj) s = struct@matRad_DoseOptimizationFunction(obj); s.penalty = obj.penalty; + end end diff --git a/matRad/optimization/+DoseObjectives/matRad_DoseVariance.m b/matRad/optimization/+DoseObjectives/matRad_DoseVariance.m new file mode 100644 index 000000000..30e0e4e54 --- /dev/null +++ b/matRad/optimization/+DoseObjectives/matRad_DoseVariance.m @@ -0,0 +1,73 @@ +classdef matRad_DoseVariance < DoseObjectives.matRad_DoseObjective +% matRad_DoseVariance Implements a variance objective for homogenous dose +% See matRad_DoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Variance'; + parameterNames = {}; + parameterTypes = {}; + end + + properties + parameters = {}; + penalty = 1; + end + + methods + function obj = matRad_DoseVariance(penalty) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@DoseObjectives.matRad_DoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(dRef) + obj.parameters{1} = dRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fDose = computeDoseObjectiveFunction(obj,dose) + fDose = obj.penalty * var(dose); + + end + + %% Calculates the Objective Function gradient + function fDoseGrad = computeDoseObjectiveGradient(obj,dose) + fDoseGrad = obj.penalty * 2/(numel(dose) - 1) * (dose - mean(dose)); + end + end + +end \ No newline at end of file diff --git a/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m b/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m index bf505662c..1a92a24c2 100644 --- a/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m +++ b/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m @@ -20,7 +20,7 @@ properties (Constant) name = 'Max DVH'; - parameterNames = {'d', 'V^{max}'}; + parameterNames = {'dose', 'V^{max}'}; parameterTypes = {'dose','numeric'}; end @@ -73,7 +73,7 @@ deviation(dose < obj.parameters{1} | dose > d_ref2) = 0; - % claculate objective function + % calculate objective function fDose = (1/numel(dose))*(deviation'*deviation); end diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m b/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m index 006295b9a..50adcf6ae 100644 --- a/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m +++ b/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m @@ -69,7 +69,7 @@ %% Calculates the Objective Function gradient function fDoseGrad = computeDoseObjectiveGradient(obj,dose) - % overdose : dose minus prefered dose + % underdose : dose minus prefered dose underdose = dose - obj.parameters{1}; % apply positive operator diff --git a/matRad/optimization/+LETdObjectives/matRad_LETdObjective.m b/matRad/optimization/+LETdObjectives/matRad_LETdObjective.m new file mode 100644 index 000000000..dcff212d8 --- /dev/null +++ b/matRad/optimization/+LETdObjectives/matRad_LETdObjective.m @@ -0,0 +1,67 @@ +classdef (Abstract) matRad_LETdObjective < matRad_DoseOptimizationFunction +% matRad_LETdObjective: Interface for optimization objectives +% This abstract base class provides the structure of optimization +% objectives like mean LET, squared deviation, EULET, LET-volume etc. +% Implementations can be found in the mLETDoseObjectives package +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Abstract, Access = public) + penalty %Optimization penalty + + end + + methods (Static) + function rob = availableRobustness() + rob = {'none','STOCH','PROB','VWWC','VWWC_INV','COWC','OWC'}; %By default, no robustness is available + end + end + + %These should be abstract methods, however Octave can't parse them. As soon + %as Octave is able to do this, they should be made abstract again + methods %(Abstract) + + %returns the objective function value for the given LET vector. Needs to be implemented in sub-classes. + function fLETd = computeLETdObjectiveFunction(obj,LETd) + error('Function needs to be implemented!'); + end + + + %returns the LET-gradient for the given LET vector. Needs to be implemented in sub-classes. + function fLETdGrad = computeLETdObjectiveGradient(obj,LETd) + error('Function needs to be implemented!'); + end + end + + methods (Access = public) + + % constructor of matRad_LETObjective + function obj = matRad_LETdObjective(varargin) + %default initialization from struct (parameters & penalty) + obj@matRad_DoseOptimizationFunction(varargin{:}); + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('LETd SquaredOverdosing does not work!') + end + + %Overloads the struct function to add Objective related information + %to output struct + function s = struct(obj) + s = struct@matRad_DoseOptimizationFunction(obj); + s.penalty = obj.penalty; + + end + end +end \ No newline at end of file diff --git a/matRad/optimization/+LETdObjectives/matRad_LETdVariance.m b/matRad/optimization/+LETdObjectives/matRad_LETdVariance.m new file mode 100644 index 000000000..09350f3a9 --- /dev/null +++ b/matRad/optimization/+LETdObjectives/matRad_LETdVariance.m @@ -0,0 +1,75 @@ +classdef matRad_LETdVariance < LETdObjectives.matRad_LETdObjective +% matRad_LETdVariance Implements a variance objective for homogenous +% LETd +% See matRad_LETdObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Variance'; + parameterNames = {}; + parameterTypes = {}; + end + + properties + parameters = {}; + penalty = 1; + end + + methods + function obj = matRad_LETdVariance(penalty) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETdObjectives.matRad_LETdObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETdRef) + obj.parameters{1} = LETdRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETd = computeLETdObjectiveFunction(obj,LETd) + fLETd = obj.penalty * var(LETd); + + end + + %% Calculates the Objective Function gradient + function fLETdGrad = computeLETdObjectiveGradient(obj,LETd) + + fLETdGrad = obj.penalty * 2/(numel(LETd) - 1) * (LETd - mean(LETd)); + end + end + +end \ No newline at end of file diff --git a/matRad/optimization/+LETdObjectives/matRad_MeanLETd.m b/matRad/optimization/+LETdObjectives/matRad_MeanLETd.m new file mode 100644 index 000000000..8050942a1 --- /dev/null +++ b/matRad/optimization/+LETdObjectives/matRad_MeanLETd.m @@ -0,0 +1,130 @@ +classdef matRad_MeanLETd < LETdObjectives.matRad_LETdObjective +% matRad_MeanLETd Implements a penalized MeanLETd objective +% See matRad_LETdObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Mean LETd'; + parameterNames = {'LETd^{ref}','f_{diff}'}; %When optimizing to a reference, one might consider using a quadratic relationship with a non-linear optimizer + parameterTypes = {'LETd',{'Linear','Quadratic'}}; + end + + properties + parameters = {0,1}; + penalty = 1; + end + + methods + function obj = matRad_MeanLETd(penalty,LETdMeanRef,fDiff) + + % if we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETdObjectives.matRad_LETdObjective(inputStruct); + + if ~initFromStruct + if nargin < 3 || ~ischar(fDiff) + fDiff = 'Linear'; + end + + fDiffIx = find(strcmp(fDiff,obj.parameterTypes{2})); + + if isempty(fDiffIx) || numel(fDiffIx) > 1 + fDiffIx = 1; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Mean LET difference function can only be %s! Using %s difference.', strjoin(obj.parameterTypes{2},' or '), obj.parameterTypes{2}{fDiffIx}); + end + + obj.parameters{2} = fDiffIx; + + + if nargin >= 2 && isscalar(LETdMeanRef) + obj.parameters{1} = LETdMeanRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + + %% Downwards compatability / set default values + %TODO: maybe move into set method for parameters + if numel(obj.parameters) < 1 + obj.parameters{1} = 0; + end + + if numel(obj.parameters) < 2 + obj.parameters{2} = 1; + end + + end + + %% Calculates the Objective Function value + function fLETd = computeLETdObjectiveFunction(obj,LETd) + switch obj.parameters{2} + case 1 + fLETd = obj.objectiveLinearDiff(LETd); + case 2 + fLETd = obj.objectiveQuadraticDiff(LETd); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid setting for %s in Mean LET Objective!',obj.parameterNames{2}); + end + end + + %% Calculates the Objective Function gradient + function fLETdGrad = computeLETdObjectiveGradient(obj,LETd) + switch obj.parameters{2} + case 1 + fLETdGrad = obj.gradientLinearDiff(LETd); + case 2 + fLETdGrad = obj.gradientQuadraticDiff(LETd); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid setting for %s in Mean Dose Objective!',obj.parameterNames{2}); + end + end + end + + methods (Access = protected) + function fLETd = objectiveQuadraticDiff(obj,LETd) + fLETd = (mean(LETd(:)) - obj.parameters{1})^2; + end + + function fLETdGrad = gradientQuadraticDiff(obj,LETd) + fLETdGrad = 2*(mean(LETd(:))-obj.parameters{1}) * ones(size(LETd(:)))/numel(LETd); + end + + function fLETd = objectiveLinearDiff(obj,LETd) + fLETd = abs(mean(LETd(:)) - obj.parameters{1}); + end + + function fLETdGrad = gradientLinearDiff(obj,LETd) + fLETdGrad = (1/numel(LETd))*sign(LETd(:)-obj.parameters{1}); + end + end + +end + diff --git a/matRad/optimization/+LETdObjectives/matRad_SquaredDeviationLETd.m b/matRad/optimization/+LETdObjectives/matRad_SquaredDeviationLETd.m new file mode 100644 index 000000000..ff49fb6ef --- /dev/null +++ b/matRad/optimization/+LETdObjectives/matRad_SquaredDeviationLETd.m @@ -0,0 +1,80 @@ +classdef matRad_SquaredDeviationLETd < LETdObjectives.matRad_LETdObjective +% matRad_SquaredDeviationLETd Implements a penalized least squares objective +% See matRad_LETdObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Deviation LETd'; + parameterNames = {'LETd^{ref}'}; + parameterTypes = {'LETd'}; + end + + properties + parameters = {60}; + penalty = 1; + end + + + methods + function obj = matRad_SquaredDeviationLETd(penalty,LETdRef) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETdObjectives.matRad_LETdObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETdRef) + obj.parameters{1} = LETdRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETd = computeLETdObjectiveFunction(obj,LETd) + % deviation : LETd minus prefered LETd + deviation = LETd - obj.parameters{1}; + % calculate objective function + fLETd = 1/numel(LETd) * (deviation'*deviation); + end + + %% Calculates the Objective Function gradient + function fLETdGrad = computeLETdObjectiveGradient(obj,LETd) + % deviation : LETd minus prefered LETd + deviation = LETd - obj.parameters{1}; + + % calculate delta + fLETdGrad = 2/numel(LETd) * deviation; + end + end + +end diff --git a/matRad/optimization/+LETdObjectives/matRad_SquaredOverdosingLETd.m b/matRad/optimization/+LETdObjectives/matRad_SquaredOverdosingLETd.m new file mode 100644 index 000000000..2975b0a6c --- /dev/null +++ b/matRad/optimization/+LETdObjectives/matRad_SquaredOverdosingLETd.m @@ -0,0 +1,83 @@ +classdef matRad_SquaredOverdosingLETd < LETdObjectives.matRad_LETdObjective +% matRad_SquaredOverdosingLETd Implements a penalized squared overdosing LETd objective +% See matRad_LETdObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Overdosing LETd'; + parameterNames = {'LETd^{max}'}; + parameterTypes = {'LETd'}; + end + + properties + parameters = {30}; + penalty = 1; + end + + methods + function obj = matRad_SquaredOverdosingLETd(penalty,LETdMax) + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETdObjectives.matRad_LETdObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETdMax) + obj.parameters{1} = LETdMax; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETd = computeLETdObjectiveFunction(obj,LETd) + % overLETd : LETd minus prefered dose + overLETd = LETd - obj.parameters{1}; + + % apply positive operator + overLETd(overLETd<0) = 0; + + % calculate objective function + fLETd = 1/numel(LETd) * (overLETd'*overLETd); + end + + %% Calculates the Objective Function gradient + function fLETdGrad = computeLETdObjectiveGradient(obj,LETd) + % overLETd : LETd minus prefered dose + overLETd = LETd - obj.parameters{1}; + + % apply positive operator + overLETd(overLETd<0) = 0; + + % calculate delta + fLETdGrad = 2 * 1/numel(LETd) * overLETd; + end + end + +end diff --git a/matRad/optimization/+LETdObjectives/matRad_SquaredUnderdosingLETd.m b/matRad/optimization/+LETdObjectives/matRad_SquaredUnderdosingLETd.m new file mode 100644 index 000000000..dabf3f312 --- /dev/null +++ b/matRad/optimization/+LETdObjectives/matRad_SquaredUnderdosingLETd.m @@ -0,0 +1,83 @@ +classdef matRad_SquaredUnderdosingLETd < LETdObjectives.matRad_LETdObjective +% matRad_SquaredUnderdosingLETd Implements a penalized squared underdosing LETd objective +% See matRad_LETdObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Underdosing LETd'; + parameterNames = {'LETd^{min}'}; + parameterTypes = {'LETd'}; + end + + properties + parameters = {60}; + penalty = 1; + end + + methods + function obj = matRad_SquaredUnderdosingLETd(penalty,LETdMin) + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETdObjectives.matRad_LETdObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETdMin) + obj.parameters{1} = LETdMin; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETd = computeLETdObjectiveFunction(obj,LETd) + % underLETd : LETd minus prefered LETd + underLETd = LETd - obj.parameters{1}; + + % apply positive operator + underLETd(underLETd>0) = 0; + + % calculate objective function + fLETd = 1/numel(LETd) * (underLETd'*underLETd); + end + + %% Calculates the Objective Function gradient + function fLETdGrad = computeLETdObjectiveGradient(obj,LETd) + % underLETd : LETd minus prefered LETd + underLETd = LETd - obj.parameters{1}; + + % apply positive operator + underLETd(underLETd>0) = 0; + + % calculate delta + fLETdGrad = 2/numel(LETd) * underLETd; + end + end + +end diff --git a/matRad/optimization/+LETxDoseObjectives/matRad_LETxDoseObjective.m b/matRad/optimization/+LETxDoseObjectives/matRad_LETxDoseObjective.m new file mode 100644 index 000000000..1edc7e434 --- /dev/null +++ b/matRad/optimization/+LETxDoseObjectives/matRad_LETxDoseObjective.m @@ -0,0 +1,66 @@ +classdef (Abstract) matRad_LETxDoseObjective < matRad_DoseOptimizationFunction +% matRad_LETxDoseObjective: Interface for optimization objectives +% This abstract base class provides the structure of optimization +% objectives like mean LETxDose, squared deviation, EULETxDose, LETxDose-volume etc. +% Implementations can be found in the LETxDoseObjectives package +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Abstract, Access = public) + penalty %Optimization penalty + + end + + methods (Static) + function rob = availableRobustness() + rob = {'none','STOCH','PROB','VWWC','VWWC_INV','COWC','OWC'}; %By default, no robustness is available + end + end + + %These should be abstract methods, however Octave can't parse them. As soon + %as Octave is able to do this, they should be made abstract again + methods %(Abstract) + + %returns the objective function value for the given LET vector. Needs to be implemented in sub-classes. + function fLETxDose = computeLETxDoseObjectiveFunction(obj,LETxDose) + error('Function needs to be implemented!'); + end + + + %returns the LET-gradient for the given LET vector. Needs to be implemented in sub-classes. + function fLETxDoseGrad = computeLETxDoseObjectiveGradient(obj,LETxDose) + error('Function needs to be implemented!'); + end + + end + + methods (Access = public) + + % constructor of matRad_LETObjective + function obj = matRad_LETxDoseObjective(varargin) + %default initialization from struct (parameters & penalty) + obj@matRad_DoseOptimizationFunction(varargin{:}); + end + + %Overloads the struct function to add Objective related information + %to output struct + function s = struct(obj) + s = struct@matRad_DoseOptimizationFunction(obj); + s.penalty = obj.penalty; + + end + end +end \ No newline at end of file diff --git a/matRad/optimization/+LETxDoseObjectives/matRad_LETxDoseVariance.m b/matRad/optimization/+LETxDoseObjectives/matRad_LETxDoseVariance.m new file mode 100644 index 000000000..09385b404 --- /dev/null +++ b/matRad/optimization/+LETxDoseObjectives/matRad_LETxDoseVariance.m @@ -0,0 +1,74 @@ +classdef matRad_LETxDoseVariance < LETxDoseObjectives.matRad_LETxDoseObjective +% matRad_LETxDoseVariance Implements a variance objective for homogenous LETxDose +% See matRad_LETxDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Variance'; + parameterNames = {}; + parameterTypes = {}; + end + + properties + parameters = {}; + penalty = 1; + end + + methods + function obj = matRad_LETxDoseVariance(penalty) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETxDoseObjectives.matRad_LETxDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETxDoseRef) + obj.parameters{1} = LETxDoseRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETxDose = computeLETxDoseObjectiveFunction(obj,LETxDose) + fLETxDose = obj.penalty * var(LETxDose); + + end + + %% Calculates the Objective Function gradient + function fLETxDoseGrad = computeLETxDoseObjectiveGradient(obj,LETxDose) + + fLETxDoseGrad = obj.penalty * 2/(numel(LETxDose) - 1) * (LETxDose - mean(LETxDose)); + end + end + +end \ No newline at end of file diff --git a/matRad/optimization/+LETxDoseObjectives/matRad_MeanLETxDose.m b/matRad/optimization/+LETxDoseObjectives/matRad_MeanLETxDose.m new file mode 100644 index 000000000..5f81552d9 --- /dev/null +++ b/matRad/optimization/+LETxDoseObjectives/matRad_MeanLETxDose.m @@ -0,0 +1,130 @@ +classdef matRad_MeanLETxDose < LETxDoseObjectives.matRad_LETxDoseObjective +% matRad_MeanLETxDose implements a penalized MeanLETxDose objective +% See matRad_LETxDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Mean LETxDose'; + parameterNames = {'LETxDose^{ref}','f_{diff}'}; %When optimizing to a reference, one might consider using a quadratic relationship with a non-linear optimizer + parameterTypes = {'LETxDose',{'Linear','Quadratic'}}; + end + + properties + parameters = {0,1}; + penalty = 1; + end + + methods + function obj = matRad_MeanLETxDose(penalty,LETxDoseRef,fDiff) + + % if we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETxDoseObjectives.matRad_LETxDoseObjective(inputStruct); + + if ~initFromStruct + if nargin < 3 || ~ischar(fDiff) + fDiff = 'Linear'; + end + + fDiffIx = find(strcmp(fDiff,obj.parameterTypes{2})); + + if isempty(fDiffIx) || numel(fDiffIx) > 1 + fDiffIx = 1; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Mean LETxDose difference function can only be %s! Using %s difference.', strjoin(obj.parameterTypes{2},' or '), obj.parameterTypes{2}{fDiffIx}); + end + + obj.parameters{2} = fDiffIx; + + + if nargin >= 2 && isscalar(LETxDoseRef) + obj.parameters{1} = LETxDoseRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + + %% Downwards compatability / set default values + + if numel(obj.parameters) < 1 + obj.parameters{1} = 0; + end + + if numel(obj.parameters) < 2 + obj.parameters{2} = 1; + end + + end + + %% Calculates the Objective Function value + function fLETxDose = computeLETxDoseObjectiveFunction(obj,LETxDose) + switch obj.parameters{2} + case 1 + fLETxDose = obj.objectiveLinearDiff(LETxDose); + case 2 + fLETxDose = obj.objectiveQuadraticDiff(LETxDose); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid setting for %s in Mean LETxDose Objective!',obj.parameterNames{2}); + end + end + + %% Calculates the Objective Function gradient + function fLETxDoseGrad = computeLETxDoseObjectiveGradient(obj,LETxDose) + switch obj.parameters{2} + case 1 + fLETxDoseGrad = obj.gradientLinearDiff(LETxDose); + case 2 + fLETxDoseGrad = obj.gradientQuadraticDiff(LETxDose); + otherwise + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid setting for %s in Mean LETxDose Objective!',obj.parameterNames{2}); + end + end + end + + methods (Access = protected) + function fLETxDose = objectiveQuadraticDiff(obj,LETxDose) + fLETxDose = (mean(LETxDose(:)) - obj.parameters{1})^2; + end + + function fLETxDoseGrad = gradientQuadraticDiff(obj,LETxDose) + fLETxDoseGrad = 2*(mean(LETxDose(:))-obj.parameters{1}) * ones(size(LETxDose(:)))/numel(LETxDose); + end + + function fLETxDose = objectiveLinearDiff(obj,LETxDose) + fLETxDose = abs(mean(LETxDose(:)) - obj.parameters{1}); + end + + function fLETxDoseGrad = gradientLinearDiff(obj,LETxDose) + fLETxDoseGrad = (1/numel(mLETxDose))*sign(LETxDose(:)-obj.parameters{1}); + end + end + +end + diff --git a/matRad/optimization/+LETxDoseObjectives/matRad_SquaredDeviationLETxDose.m b/matRad/optimization/+LETxDoseObjectives/matRad_SquaredDeviationLETxDose.m new file mode 100644 index 000000000..ea1b43323 --- /dev/null +++ b/matRad/optimization/+LETxDoseObjectives/matRad_SquaredDeviationLETxDose.m @@ -0,0 +1,79 @@ +classdef matRad_SquaredDeviationLETxDose < LETxDoseObjectives.matRad_LETxDoseObjective +% matRad_SquaredDeviationLETxDose implements a penalized least squares objective +% See matRad_LETxDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Deviation LETxDose'; + parameterNames = {'LETxd^{ref}'}; + parameterTypes = {'LETxd'}; + end + + properties + parameters = {60}; + penalty = 1; + end + + methods + function obj = matRad_SquaredDeviationLETxDose(penalty,LETxDoseRef) + + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETxDoseObjectives.matRad_LETxDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETxDoseRef) + obj.parameters{1} = LETxDoseRef; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETxDose = computeLETxDoseObjectiveFunction(obj,LETxDose) + % deviation : mLETDose minus prefered mLETDose + deviation = LETxDose - obj.parameters{1}; + % calculate objective function + fLETxDose = 1/numel(LETxDose) * (deviation'*deviation); + end + + %% Calculates the Objective Function gradient + function fLETxDoseGrad = computeLETxDoseObjectiveGradient(obj,LETxDose) + % deviation : LETxDose minus prefered LETxDose + deviation = LETxDose - obj.parameters{1}; + + % calculate delta + fLETxDoseGrad = 2 * 1/numel(LETxDose) * deviation; + end + end + +end \ No newline at end of file diff --git a/matRad/optimization/+LETxDoseObjectives/matRad_SquaredOverdosingLETxDose.m b/matRad/optimization/+LETxDoseObjectives/matRad_SquaredOverdosingLETxDose.m new file mode 100644 index 000000000..e9ac62b86 --- /dev/null +++ b/matRad/optimization/+LETxDoseObjectives/matRad_SquaredOverdosingLETxDose.m @@ -0,0 +1,83 @@ +classdef matRad_SquaredOverdosingLETxDose < LETxDoseObjectives.matRad_LETxDoseObjective +% matRad_SquaredOverdosingLETxDose implements a penalized squared overdosing LETxDose objective +% See matRad_LETxDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Overdosing LETxDose'; + parameterNames = {'LETxd^{max}'}; + parameterTypes = {'LETxd'}; + end + + properties + parameters = {30}; + penalty = 1; + end + + methods + function obj = matRad_SquaredOverdosingLETxDose(penalty,LETxDoseMax) + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETxDoseObjectives.matRad_LETxDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin == 2 && isscalar(LETxDoseMax) + obj.parameters{1} = LETxDoseMax; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETxDose = computeLETxDoseObjectiveFunction(obj,LETxDose) + % overLETxDose : LETxDose minus prefered LETxDose + overLETxDose = LETxDose - obj.parameters{1}; + + % apply positive operator + overLETxDose(overLETxDose<0) = 0; + + % calculate objective function + fLETxDose = 1/numel(LETxDose) * (overLETxDose'*overLETxDose); + end + + %% Calculates the Objective Function gradient + function fLETxDoseGrad = computeLETxDoseObjectiveGradient(obj,LETxDose) + % overLETxDose : LETxDose minus prefered LETxDose + overLETxDose = LETxDose - obj.parameters{1}; + + % apply positive operator + overLETxDose(overLETxDose<0) = 0; + + % calculate delta + fLETxDoseGrad = 2 * 1/numel(LETxDose) * overLETxDose; + end + end + +end diff --git a/matRad/optimization/+LETxDoseObjectives/matRad_SquaredUnderdosingLETxDose.m b/matRad/optimization/+LETxDoseObjectives/matRad_SquaredUnderdosingLETxDose.m new file mode 100644 index 000000000..51fbc65b3 --- /dev/null +++ b/matRad/optimization/+LETxDoseObjectives/matRad_SquaredUnderdosingLETxDose.m @@ -0,0 +1,83 @@ +classdef matRad_SquaredUnderdosingLETxDose < LETxDoseObjectives.matRad_LETxDoseObjective +% matRad_SquaredUnderdosingLETxDose implements a penalized squared underdosing LETxDose objective +% See matRad_LETxDoseObjective for interface description +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (Constant) + name = 'Squared Underdosing LETxDose'; + parameterNames = {'LETxd^{min}'}; + parameterTypes = {'LETxd'}; + end + + properties + parameters = {60}; + penalty = 1; + end + + methods + function obj = matRad_SquaredUnderdosingLETxDose(penalty,LETxDoseMin) + %If we have a struct in first argument + if nargin == 1 && isstruct(penalty) + inputStruct = penalty; + initFromStruct = true; + else + initFromStruct = false; + inputStruct = []; + end + + %Call Superclass Constructor (for struct initialization) + obj@LETxDoseObjectives.matRad_LETxDoseObjective(inputStruct); + + %now handle initialization from other parameters + if ~initFromStruct + if nargin >= 2 && isscalar(LETxDoseMin) + obj.parameters{1} = LETxDoseMin; + end + + if nargin >= 1 && isscalar(penalty) + obj.penalty = penalty; + end + end + end + + %% Calculates the Objective Function value + function fLETxDose = computeLETxDoseObjectiveFunction(obj,LETxDose) + % underLETxDose : LETxDose minus prefered LETxDose + underLETxDose = LETxDose - obj.parameters{1}; + + % apply positive operator + underLETxDose(underLETxDose>0) = 0; + + % calculate objective function + fLETxDose = 1/numel(LETxDose) * (underLETxDose'*underLETxDose); + end + + %% Calculates the Objective Function gradient + function fLETxDoseGrad = computeLETxDoseObjectiveGradient(obj,LETxDose) + % underLETxDose : LETxDose minus prefered LETxDose + underLETxDose = LETxDose - obj.parameters{1}; + + % apply positive operator + underLETxDose(underLETxDose>0) = 0; + + % calculate delta + fLETxDoseGrad = 2/numel(LETxDose) * underLETxDose; + end + end + +end \ No newline at end of file diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m index 529e20f3b..6146c8192 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m @@ -30,12 +30,19 @@ minimumW = NaN; maximumW = NaN; + + dirtyDoseBP + LETxDoseBP + LETdBP end methods function obj = matRad_OptimizationProblem(backProjection) obj.BP = backProjection; - end + obj.dirtyDoseBP = []; + obj.LETxDoseBP = []; + obj.LETdBP = []; + end %Objective function declaration fVal = matRad_objectiveFunction(optiProb,w,dij,cst) diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m index 17d4abab3..ff67f4a7b 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m @@ -34,7 +34,6 @@ % get current dose / effect / RBExDose vector optiProb.BP.compute(dij,w); d = optiProb.BP.GetResult(); - % get probabilistic quantities (nearly no overhead if empty) [dExp,dOmega] = optiProb.BP.GetResultProb(); @@ -48,6 +47,32 @@ [fullScen{:}] = ind2sub(size(d),useScen); contourScen = fullScen{1}; +%Manage other quantities + +if ~isempty(optiProb.dirtyDoseBP) + optiProb.dirtyDoseBP.compute(dij,w); + dD = optiProb.dirtyDoseBP.GetResult(); + + % get probabilistic quantities (nearly no overhead if empty) + [dDExp,dDOmega] = optiProb.dirtyDoseBP.GetResultProb(); +end + +if ~isempty(optiProb.LETxDoseBP) + optiProb.LETxDoseBP.compute(dij,w); + LxD = optiProb.LETxDoseBP.GetResult(); + + % also get probabilistic quantities (nearly no overhead if empty) + [LxDExp,LxDOmega] = optiProb.LETxDoseBP.GetResultProb(); +end + +if ~isempty(optiProb.LETdBP) + optiProb.LETdBP.compute(dij,w); + Ld = optiProb.LETdBP.GetResult(); + + % also get probabilistic quantities (nearly no overhead if empty) + [LdExp,LdOmega] = optiProb.LETdBP.GetResultProb(); +end + % initialize f f = 0; @@ -56,24 +81,24 @@ % compute objective function for every VOI. for i = 1:size(cst,1) - + % Only take OAR or target VOI. if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) - + % loop over the number of constraints for the current VOI for j = 1:numel(cst{i,6}) - + objective = cst{i,6}{j}; - + % only perform gradient computations for objectiveectives if isa(objective,'DoseObjectives.matRad_DoseObjective') - + % rescale dose parameters to biological optimization quantity if required objective = optiProb.BP.setBiologicalDosePrescriptions(objective,cst{i,5}.alphaX,cst{i,5}.betaX); % retrieve the robustness type robustness = objective.robustness; - + switch robustness case 'none' % if conventional opt: just sum objectives of nominal dose for ixScen = useNominalCtScen @@ -170,7 +195,7 @@ case 'OWC' % objective-wise worst case considers the worst individual objective function value f_OWC = zeros(numel(useScen),1); - + for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); @@ -192,14 +217,391 @@ fMax = max(f_OWC); end f = f + fMax; - + otherwise matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); - end %robustness type - end % objective check - end %objective loop - end %empty check + end %robustness type + end % objective check + + if ~isempty(optiProb.dirtyDoseBP) && isa(objective,'DirtyDoseObjectives.matRad_DirtyDoseObjective') + + % retrieve the robustness type + robustness = objective.robustness; + + switch robustness + case 'none' % if conventional opt: just sum objectives of nominal dirty dose + for ixScen = useNominalCtScen + d_i = dD{ixScen}(cst{i,4}{useScen(ixScen)}); + f = f + objective.penalty * objective.computeDirtyDoseObjectiveFunction(d_i); + end + + case 'STOCH' % if prob opt: sum up expectation value of objectives + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = dD{ixScen}(cst{i,4}{ixContour}); + + f = f + scenProb(s) * objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + + end + + case 'PROB' % if prob opt: sum up expectation value of objectives + + d_i = dDExp{1}(cst{i,4}{1}); + + f = f + objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + + p = objective.penalty/numel(cst{i,4}{1}); + + % only one variance term per VOI + if j == 1 + f = f + p * w' * dDOmega{i,1}; + end + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [dD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_max; + elseif isequal(cst{i,3},'TARGET') + d_i = d_min; + end + + f = f + objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + + case 'VWWC_INV' %inverse voxel-wise conformitiy - consider the maximum and minimum dose in the target and optimize the dose conformity + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispWarning('4D inverted VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [dD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_min; + elseif isequal(cst{i,3},'TARGET') + d_i = d_max; + end + + f = f + objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = dD{ixScen}(cst{i,4}{ixContour}); + + f_COWC(s) = f_COWC(s) + objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + end + + case 'OWC' % objective-wise worst case considers the worst individual objective function value + + f_OWC = zeros(numel(useScen),1); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = dD{ixScen}(cst{i,4}{ixContour}); + f_OWC(s) = objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + end + + % compute the maximum objective function value + switch optiProb.useMaxApprox + case 'logsumexp' + fMax = optiProb.logSumExp(f_OWC); + case 'pnorm' + fMax = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + fMax = max(f_OWC); + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + fMax = max(f_OWC); + end + f = f + fMax; + end + end + if ~isempty(optiProb.LETxDoseBP) && isa(objective,'LETxDoseObjectives.matRad_LETxDoseObjective') + + % retrieve the robustness type + robustness = objective.robustness; + + switch robustness + case 'none' % if conventional opt: just sum objectives of nominal dirty dose + for ixScen = useNominalCtScen + d_i = LxD{ixScen}(cst{i,4}{useScen(ixScen)}); + f = f + objective.penalty * objective.computeLETxDoseObjectiveFunction(d_i); + end + + case 'STOCH' % if prob opt: sum up expectation value of objectives + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + + f = f + scenProb(s) * objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + + end + + case 'PROB' % if prob opt: sum up expectation value of objectives + + d_i = LxDExp{1}(cst{i,4}{1}); + + f = f + objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + + p = objective.penalty/numel(cst{i,4}{1}); + + % only one variance term per VOI + if j == 1 + f = f + p * w' * LxDOmega{i,1}; + end + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [LxD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_max; + elseif isequal(cst{i,3},'TARGET') + d_i = d_min; + end + + f = f + objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + + case 'VWWC_INV' %inverse voxel-wise conformitiy - consider the maximum and minimum dose in the target and optimize the dose conformity + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispWarning('4D inverted VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [LxD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_min; + elseif isequal(cst{i,3},'TARGET') + d_i = d_max; + end + + f = f + objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + + f_COWC(s) = f_COWC(s) + objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + end + + case 'OWC' % objective-wise worst case considers the worst individual objective function value + + f_OWC = zeros(numel(useScen),1); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + f_OWC(s) = objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + end + + % compute the maximum objective function value + switch optiProb.useMaxApprox + case 'logsumexp' + fMax = optiProb.logSumExp(f_OWC); + case 'pnorm' + fMax = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + fMax = max(f_OWC); + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + fMax = max(f_OWC); + end + f = f + fMax; + end + end + if ~isempty(optiProb.LETdBP) && isa(objective,'LETdObjectives.matRad_LETdObjective') + + % retrieve the robustness type + robustness = objective.robustness; + + switch robustness + case 'none' % if conventional opt: just sum objectives of nominal dirty dose + for ixScen = useNominalCtScen + d_i = Ld{ixScen}(cst{i,4}{useScen(ixScen)}); + f = f + objective.penalty * objective.computeLETdObjectiveFunction(d_i); + end + + case 'STOCH' % if prob opt: sum up expectation value of objectives + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + + f = f + scenProb(s) * objective.penalty*objective.computeLETdObjectiveFunction(d_i); + + end + + case 'PROB' % if prob opt: sum up expectation value of objectives + + d_i = LdExp{1}(cst{i,4}{1}); + + f = f + objective.penalty*objective.computeLETdObjectiveFunction(d_i); + + p = objective.penalty/numel(cst{i,4}{1}); + + % only one variance term per VOI + if j == 1 + f = f + p * w' * LdOmega{i,1}; + end + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [Ld{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_max; + elseif isequal(cst{i,3},'TARGET') + d_i = d_min; + end + + f = f + objective.penalty*objective.computeLETdObjectiveFunction(d_i); + + case 'VWWC_INV' %inverse voxel-wise conformitiy - consider the maximum and minimum dose in the target and optimize the dose conformity + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispWarning('4D inverted VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [Ld{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_min; + elseif isequal(cst{i,3},'TARGET') + d_i = d_max; + end + + f = f + objective.penalty*objective.computeLETdObjectiveFunction(d_i); + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + + f_COWC(s) = f_COWC(s) + objective.penalty*objective.computeLETdObjectiveFunction(d_i); + end + + case 'OWC' % objective-wise worst case considers the worst individual objective function value + + f_OWC = zeros(numel(useScen),1); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + f_OWC(s) = objective.penalty*objective.computeLETdObjectiveFunction(d_i); + end + + % compute the maximum objective function value + switch optiProb.useMaxApprox + case 'logsumexp' + fMax = optiProb.logSumExp(f_OWC); + case 'pnorm' + fMax = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + fMax = max(f_OWC); + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + fMax = max(f_OWC); + end + f = f + fMax; + end + end + + end %objective loop + end %empty check end %cst structure loop diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m index 1d92ff56c..757561d2d 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m @@ -56,6 +56,40 @@ doseGradient = cell(size(dij.physicalDose)); doseGradient(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; +%Manage other quantites +if ~isempty(optiProb.dirtyDoseBP) + optiProb.dirtyDoseBP.compute(dij,w); + dD = optiProb.dirtyDoseBP.GetResult(); + + % also get probabilistic quantities (nearly no overhead if empty) + [dDExp,dDOmega] = optiProb.dirtyDoseBP.GetResultProb(); + + dirtyDoseGradient = cell(size(dij.dirtyDose)); + dirtyDoseGradient(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; +end + +if ~isempty(optiProb.LETxDoseBP) + optiProb.LETxDoseBP.compute(dij,w); + LxD = optiProb.LETxDoseBP.GetResult(); + + % also get probabilistic quantities (nearly no overhead if empty) + [LxDExp,LxDOmega] = optiProb.LETxDoseBP.GetResultProb(); + + LETxDoseGradient = cell(size(dij.mLETDose)); + LETxDoseGradient(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; +end + +if ~isempty(optiProb.LETdBP) + optiProb.LETdBP.compute(dij,w); % the problem is here + Ld = optiProb.LETdBP.GetResult(); + + % also get probabilistic quantities (nearly no overhead if empty) + [LdExp,LdOmega] = optiProb.LETdBP.GetResultProb(); + + LETdGradient = cell(size(dij.mLETDose)); + LETdGradient(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; +end + %For probabilistic optimization vOmega = 0; @@ -64,25 +98,25 @@ % compute objective function for every VOI. for i = 1:size(cst,1) - + % Only take OAR or target VOI. if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) - + % loop over the number of constraints and objectives for the current VOI for j = 1:numel(cst{i,6}) - + %Get current optimization function objective = cst{i,6}{j}; - + % only perform gradient computations for objectives if isa(objective,'DoseObjectives.matRad_DoseObjective') - + % retrieve the robustness type robustness = objective.robustness; - + % rescale dose parameters to biological optimization quantity if required objective = optiProb.BP.setBiologicalDosePrescriptions(objective,cst{i,5}.alphaX,cst{i,5}.betaX); - + switch robustness case 'none' % if conventional opt: just sum objectiveectives of nominal dose for s = useNominalCtScen @@ -96,28 +130,28 @@ for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); - + d_i = d{ixScen}(cst{i,4}{ixContour}); - + doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + ... (objective.penalty*objective.computeDoseObjectiveGradient(d_i) * scenProb(s)); - + end - + case 'PROB' % use the expectation value and the integral variance influence matrix %First check the speficic cache for probabilistic if ~exist('doseGradientExp','var') doseGradientExp{1} = zeros(dij.doseGrid.numOfVoxels,1); end - + d_i = dExp{1}(cst{i,4}{1}); - + doseGradientExp{1}(cst{i,4}{1}) = doseGradientExp{1}(cst{i,4}{1}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); - + p = objective.penalty/numel(cst{i,4}{1}); - + vOmega = vOmega + p * dOmega{i,1}; - + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR contourIx = unique(contourScen); if ~isscalar(contourIx) @@ -125,41 +159,41 @@ % not yet implemented matRad_cfg.dispError('4D VWWC optimization is currently not supported'); end - + % prepare min/max dose vector for voxel-wise worst case if ~exist('d_tmp','var') d_tmp = [d{useScen}]; end - + d_Scen = d_tmp(cst{i,4}{contourIx},:); [d_max,max_ix] = max(d_Scen,[],2); [d_min,min_ix] = min(d_Scen,[],2); - + if isequal(cst{i,3},'OAR') d_i = d_max; elseif isequal(cst{i,3},'TARGET') d_i = d_min; end - + if any(isnan(d_i)) matRad_cfg.dispWarning('%d NaN values in gradient.',numel(isnan(d_i))); end - + deltaTmp = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - + for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); - + if isequal(cst{i,3},'OAR') currWcIx = double(max_ix == s); elseif isequal(cst{i,3},'TARGET') currWcIx = double(min_ix == s); end - + doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; end - + case 'VWWC_INV' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR contourIx = unique(contourScen); if ~isscalar(contourIx) @@ -167,79 +201,80 @@ % not yet implemented matRad_cfg.dispError('4D VWWC optimization is currently not supported'); end - + % prepare min/max dose vector for voxel-wise worst case if ~exist('d_tmp','var') d_tmp = [d{useScen}]; end - + d_Scen = d_tmp(cst{i,4}{1},:); [d_max,max_ix] = max(d_Scen,[],2); [d_min,min_ix] = min(d_Scen,[],2); - + if isequal(cst{i,3},'OAR') d_i = d_min; elseif isequal(cst{i,3},'TARGET') d_i = d_max; end - + if any(isnan(d_i)) matRad_cfg.dispWarning('%d NaN values in gradFuncWrapper.',numel(isnan(d_i))); end - + deltaTmp = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - + for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); - + if isequal(cst{i,3},'OAR') currWcIx = double(min_ix == s); elseif isequal(cst{i,3},'TARGET') currWcIx = double(max_ix == s); end - + doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; end - + case 'COWC' % composite worst case consideres ovarall the worst objective function value %First check the speficic cache for COWC if ~exist('delta_COWC','var') delta_COWC = cell(size(doseGradient)); delta_COWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; end - + for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); - + d_i = d{ixScen}(cst{i,4}{ixContour}); - + f_COWC(ixScen) = f_COWC(ixScen) + objective.penalty*objective.computeDoseObjectiveFunction(d_i); delta_COWC{ixScen}(cst{i,4}{ixContour}) = delta_COWC{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); end - + case 'OWC' % objective-wise worst case consideres the worst individual objective function value %First check the speficic cache for COWC f_OWC = zeros(size(doseGradient)); - + if ~exist('delta_OWC','var') delta_OWC = cell(size(doseGradient)); delta_OWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; end - + for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); - + d_i = d{ixScen}(cst{i,4}{ixContour}); - + f_OWC(ixScen) = objective.penalty*objective.computeDoseObjectiveFunction(d_i); - + delta_OWC{ixScen}(cst{i,4}{ixContour}) = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - + end - + + switch optiProb.useMaxApprox case 'logsumexp' [~,fGrad] = optiProb.logSumExp(f_OWC); @@ -255,7 +290,7 @@ fGrad = zeros(size(f_OWC)); fGrad(ix) = 1; end - + for s = 1:numel(useScen) ixScen = useScen(s); ixContour = contourScen(s); @@ -263,54 +298,695 @@ doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + fGrad(ixScen)*delta_OWC{ixScen}(cst{i,4}{ixContour}); end end - + otherwise matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); - - end %robustness type - end % objective check - end %objective loop - end %empty check -end %cst structure loop -if exist('delta_COWC','var') - switch optiProb.useMaxApprox - case 'logsumexp' - [~,fGrad] = optiProb.logSumExp(f_COWC); - case 'pnorm' - [~,fGrad] = optiProb.pNorm(f_COWC,numel(useScen)); - case 'none' - [~,ixCurrWC] = max(f_COWC(:)); - fGrad = zeros(size(f_COWC)); - fGrad(ixCurrWC) = 1; - case 'otherwise' - matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); - [~,ixCurrWC] = max(f_COWC(:)); - fGrad = zeros(size(f_COWC)); - fGrad(ixCurrWC) = 1; - end - - for s = 1:numel(useScen) - ixScen = useScen(s); - if fGrad(ixScen) ~= 0 - doseGradient{ixScen} = doseGradient{ixScen} + fGrad(ixScen)*delta_COWC{ixScen}; - end - end -end + end %robustness type + end + if ~isempty(optiProb.dirtyDoseBP) && isa(objective,'DirtyDoseObjectives.matRad_DirtyDoseObjective') + % only perform gradient computations for objectives -weightGradient = zeros(dij.totalNumOfBixels,1); + % retrieve the robustness type + robustness = objective.robustness; -optiProb.BP.computeGradient(dij,doseGradient,w); -g = optiProb.BP.GetGradient(); + switch robustness + case 'none' % if conventional opt: just sum objectiveectives of nominal dose + for s = useNominalCtScen + ixScen = useScen(s); + ixContour = contourScen(s); + d_i = dD{ixScen}(cst{i,4}{ixContour}); + %add to dose gradient + dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) = dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i); + end + case 'STOCH' % perform stochastic optimization with weighted / random scenarios + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); -for s = 1:numel(useScen) - weightGradient = weightGradient + g{useScen(s)}; -end + d_i = dD{ixScen}(cst{i,4}{ixContour}); -if vOmega ~= 0 - optiProb.BP.computeGradientProb(dij,doseGradientExp,vOmega,w); - gProb = optiProb.BP.GetGradientProb(); - - %Only implemented for first scenario now - weightGradient = weightGradient + gProb{1}; -end + dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) = dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) + ... + (objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i) * scenProb(s)); + + end + + case 'PROB' % use the expectation value and the integral variance influence matrix + %First check the speficic cache for probabilistic + if ~exist('dirtyDoseGradientExp','var') + dirtyDoseGradientExp{1} = zeros(dij.doseGrid.numOfVoxels,1); + end + + d_i = dDExp{1}(cst{i,4}{1}); + + dirtyDoseGradientExp{1}(cst{i,4}{1}) = dirtyDoseGradientExp{1}(cst{i,4}{1}) + objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i); + + p = objective.penalty/numel(cst{i,4}{1}); + + vOmega = vOmega + p * dDOmega{i,1}; + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [dD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_max; + elseif isequal(cst{i,3},'TARGET') + d_i = d_min; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradient.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{i,3},'OAR') + currWcIx = double(max_ix == s); + elseif isequal(cst{i,3},'TARGET') + currWcIx = double(min_ix == s); + end + + dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) = dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'VWWC_INV' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [dD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{1},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_min; + elseif isequal(cst{i,3},'TARGET') + d_i = d_max; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradFuncWrapper.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{i,3},'OAR') + currWcIx = double(min_ix == s); + elseif isequal(cst{i,3},'TARGET') + currWcIx = double(max_ix == s); + end + + dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) = dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + %First check the speficic cache for COWC + if ~exist('delta_COWC','var') + delta_COWC = cell(size(dirtyDoseGradient)); + delta_COWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = dD{ixScen}(cst{i,4}{ixContour}); + + f_COWC(ixScen) = f_COWC(ixScen) + objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + delta_COWC{ixScen}(cst{i,4}{ixContour}) = delta_COWC{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i); + end + + case 'OWC' % objective-wise worst case consideres the worst individual objective function value + %First check the speficic cache for COWC + f_OWC = zeros(size(dirtyDoseGradient)); + + if ~exist('delta_OWC','var') + delta_OWC = cell(size(dirtyDoseGradient)); + delta_OWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = dD{ixScen}(cst{i,4}{ixContour}); + + f_OWC(ixScen) = objective.penalty*objective.computeDirtyDoseObjectiveFunction(d_i); + + delta_OWC{ixScen}(cst{i,4}{ixContour}) = objective.penalty*objective.computeDirtyDoseObjectiveGradient(d_i); + + end + + switch optiProb.useMaxApprox + case 'logsumexp' + [~,fGrad] = optiProb.logSumExp(f_OWC); + case 'pnorm' + [~,fGrad] = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + if fGrad(ixScen ) ~= 0 + dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) = dirtyDoseGradient{ixScen}(cst{i,4}{ixContour}) + fGrad(ixScen)*delta_OWC{ixScen}(cst{i,4}{ixContour}); + end + end + + otherwise + matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); + + end %robustness type + end + if ~isempty(optiProb.LETxDoseBP) && isa(objective,'LETxDoseObjectives.matRad_LETxDoseObjective') + % only perform gradient computations for objectives + + % retrieve the robustness type + robustness = objective.robustness; + + switch robustness + case 'none' % if conventional opt: just sum objectiveectives of nominal dose + for s = useNominalCtScen + ixScen = useScen(s); + ixContour = contourScen(s); + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + %add to dose gradient + LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) = LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i); + end + case 'STOCH' % perform stochastic optimization with weighted / random scenarios + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + + LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) = LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) + ... + (objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i) * scenProb(s)); + + end + + case 'PROB' % use the expectation value and the integral variance influence matrix + %First check the speficic cache for probabilistic + if ~exist('LETxDoseGradientExp','var') + LETxDoseGradientExp{1} = zeros(dij.doseGrid.numOfVoxels,1); + end + + d_i = LxDExp{1}(cst{i,4}{1}); + + LETxDoseGradientExp{1}(cst{i,4}{1}) = LETxDoseGradientExp{1}(cst{i,4}{1}) + objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i); + + p = objective.penalty/numel(cst{i,4}{1}); + + vOmega = vOmega + p * LxDOmega{i,1}; + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [LxD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_max; + elseif isequal(cst{i,3},'TARGET') + d_i = d_min; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradient.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{i,3},'OAR') + currWcIx = double(max_ix == s); + elseif isequal(cst{i,3},'TARGET') + currWcIx = double(min_ix == s); + end + + LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) = LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'VWWC_INV' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [LxD{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{1},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_min; + elseif isequal(cst{i,3},'TARGET') + d_i = d_max; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradFuncWrapper.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{i,3},'OAR') + currWcIx = double(min_ix == s); + elseif isequal(cst{i,3},'TARGET') + currWcIx = double(max_ix == s); + end + + LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) = LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + %First check the speficic cache for COWC + if ~exist('delta_COWC','var') + delta_COWC = cell(size(LETxDoseGradient)); + delta_COWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + + f_COWC(ixScen) = f_COWC(ixScen) + objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + delta_COWC{ixScen}(cst{i,4}{ixContour}) = delta_COWC{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i); + end + + case 'OWC' % objective-wise worst case consideres the worst individual objective function value + %First check the speficic cache for COWC + f_OWC = zeros(size(LETxDoseGradient)); + + if ~exist('delta_OWC','var') + delta_OWC = cell(size(LETxDoseGradient)); + delta_OWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = LxD{ixScen}(cst{i,4}{ixContour}); + + f_OWC(ixScen) = objective.penalty*objective.computeLETxDoseObjectiveFunction(d_i); + + delta_OWC{ixScen}(cst{i,4}{ixContour}) = objective.penalty*objective.computeLETxDoseObjectiveGradient(d_i); + + end + + switch optiProb.useMaxApprox + case 'logsumexp' + [~,fGrad] = optiProb.logSumExp(f_OWC); + case 'pnorm' + [~,fGrad] = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + if fGrad(ixScen ) ~= 0 + LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) = LETxDoseGradient{ixScen}(cst{i,4}{ixContour}) + fGrad(ixScen)*delta_OWC{ixScen}(cst{i,4}{ixContour}); + end + end + + otherwise + matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); + + end %robustness type + end + if ~isempty(optiProb.LETdBP) && isa(objective,'LETdObjectives.matRad_LETdObjective') + % only perform gradient computations for objectives + + % retrieve the robustness type + robustness = objective.robustness; + + switch robustness + case 'none' % if conventional opt: just sum objectiveectives of nominal dose + for s = useNominalCtScen + ixScen = useScen(s); + ixContour = contourScen(s); + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + %add to dose gradient + LETdGradient{ixScen}(cst{i,4}{ixContour}) = LETdGradient{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeLETdObjectiveGradient(d_i); + end + case 'STOCH' % perform stochastic optimization with weighted / random scenarios + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + + LETdGradient{ixScen}(cst{i,4}{ixContour}) = LETdGradient{ixScen}(cst{i,4}{ixContour}) + ... + (objective.penalty*objective.computeLETdObjectiveGradient(d_i) * scenProb(s)); + + end + + case 'PROB' % use the expectation value and the integral variance influence matrix + %First check the speficic cache for probabilistic + if ~exist('LETdGradientExp','var') + LETdGradientExp{1} = zeros(dij.doseGrid.numOfVoxels,1); + end + + d_i = LdExp{1}(cst{i,4}{1}); + + LETdGradientExp{1}(cst{i,4}{1}) = LETdGradientExp{1}(cst{i,4}{1}) + objective.penalty*objective.computeLETdObjectiveGradient(d_i); + + p = objective.penalty/numel(cst{i,4}{1}); + + vOmega = vOmega + p * LdOmega{i,1}; + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [Ld{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{contourIx},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_max; + elseif isequal(cst{i,3},'TARGET') + d_i = d_min; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradient.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeLETdObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{i,3},'OAR') + currWcIx = double(max_ix == s); + elseif isequal(cst{i,3},'TARGET') + currWcIx = double(min_ix == s); + end + + LETdGradient{ixScen}(cst{i,4}{ixContour}) = LETdGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'VWWC_INV' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [Ld{useScen}]; + end + + d_Scen = d_tmp(cst{i,4}{1},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{i,3},'OAR') + d_i = d_min; + elseif isequal(cst{i,3},'TARGET') + d_i = d_max; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradFuncWrapper.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeLETdObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{i,3},'OAR') + currWcIx = double(min_ix == s); + elseif isequal(cst{i,3},'TARGET') + currWcIx = double(max_ix == s); + end + + LETdGradient{ixScen}(cst{i,4}{ixContour}) = LETdGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + %First check the speficic cache for COWC + if ~exist('delta_COWC','var') + delta_COWC = cell(size(LETdGradient)); + delta_COWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + + f_COWC(ixScen) = f_COWC(ixScen) + objective.penalty*objective.computeLETdObjectiveFunction(d_i); + delta_COWC{ixScen}(cst{i,4}{ixContour}) = delta_COWC{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeLETdObjectiveGradient(d_i); + end + + case 'OWC' % objective-wise worst case consideres the worst individual objective function value + %First check the speficic cache for COWC + f_OWC = zeros(size(LETdGradient)); + + if ~exist('delta_OWC','var') + delta_OWC = cell(size(LETdGradient)); + delta_OWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = Ld{ixScen}(cst{i,4}{ixContour}); + + f_OWC(ixScen) = objective.penalty*objective.computeLETdObjectiveFunction(d_i); + + delta_OWC{ixScen}(cst{i,4}{ixContour}) = objective.penalty*objective.computeLETdObjectiveGradient(d_i); + + end + + switch optiProb.useMaxApprox + case 'logsumexp' + [~,fGrad] = optiProb.logSumExp(f_OWC); + case 'pnorm' + [~,fGrad] = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + if fGrad(ixScen ) ~= 0 + LETdGradient{ixScen}(cst{i,4}{ixContour}) = LETdGradient{ixScen}(cst{i,4}{ixContour}) + fGrad(ixScen)*delta_OWC{ixScen}(cst{i,4}{ixContour}); + end + end + + otherwise + matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); + + end %robustness type + end % objective check + + end %objective loop + end %empty check +end %cst structure loop + +if exist('delta_COWC','var') + switch optiProb.useMaxApprox + case 'logsumexp' + [~,fGrad] = optiProb.logSumExp(f_COWC); + case 'pnorm' + [~,fGrad] = optiProb.pNorm(f_COWC,numel(useScen)); + case 'none' + [~,ixCurrWC] = max(f_COWC(:)); + fGrad = zeros(size(f_COWC)); + fGrad(ixCurrWC) = 1; + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + [~,ixCurrWC] = max(f_COWC(:)); + fGrad = zeros(size(f_COWC)); + fGrad(ixCurrWC) = 1; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + if fGrad(ixScen) ~= 0 + if ~isempty(optiProb.dirtyDoseBP) + dirtyDoseGradient{ixScen} = dirtyDoseGradient{ixScen} + fGrad{ixScen} * delta_COWC{ixScen}; + end + if ~isempty(optiProb.LETxDoseBP) + LETxDoseGradient{ixScen} = LETxDoseGradient{ixScen} + fGrad{ixScen} * delta_COWC{ixScen}; + end + if ~isempty(optiProb.LETdBP) + LETdGradient{ixScen} = LETdGradient{ixScen} + fGrad{ixScen} * delta_COWC{ixScen}; + end + if ~isempty(optiProb.BP) + doseGradient{ixScen} = doseGradient{ixScen} + fGrad(ixScen)*delta_COWC{ixScen}; + end + end + end +end + +weightGradient = zeros(dij.totalNumOfBixels,1); + +%Manage Dose Gradient +optiProb.BP.computeGradient(dij,doseGradient,w); +g = optiProb.BP.GetGradient(); + +for s = 1:numel(useScen) + weightGradient = weightGradient + g{useScen(s)}; +end + +if vOmega ~= 0 + optiProb.BP.computeGradientProb(dij,doseGradientExp,vOmega,w); + gProb = optiProb.BP.GetGradientProb(); + + %Only implemented for first scenario now + weightGradient = weightGradient + gProb{1}; +end + +%Manage Dirty Dose Gradient +if ~isempty(optiProb.dirtyDoseBP) + optiProb.dirtyDoseBP.computeGradient(dij,dirtyDoseGradient,w); + g = optiProb.dirtyDoseBP.GetGradient(); + + for s = 1:numel(useScen) + weightGradient = weightGradient + g{useScen(s)}; + end + + if vOmega ~= 0 + optiProb.dirtyDoseBP.computeGradientProb(dij,dirtyDoseGradientExp,vOmega,w); + gProb = optiProb.dirtyDoseBP.GetGradientProb(); + + %Only implemented for first scenario now + weightGradient = weightGradient + gProb{1}; + end +end +if ~isempty(optiProb.LETxDoseBP) + optiProb.LETxDoseBP.computeGradient(dij,LETxDoseGradient,w); + g = optiProb.LETxDoseBP.GetGradient(); + + for s = 1:numel(useScen) + weightGradient = weightGradient + g{useScen(s)}; + end + + if vOmega ~= 0 + optiProb.LETxDoseBP.computeGradientProb(dij,LETxDoseGradientExp,vOmega,w); + gProb = optiProb.LETxDoseBP.GetGradientProb(); + + %Only implemented for first scenario now + weightGradient = weightGradient + gProb{1}; + end +end +if ~isempty(optiProb.LETdBP) + optiProb.LETdBP.computeGradient(dij,LETdGradient,w); + g = optiProb.LETdBP.GetGradient(); + + for s = 1:numel(useScen) + weightGradient = weightGradient + g{useScen(s)}; + end + + if vOmega ~= 0 + optiProb.LETdBP.computeGradientProb(dij,LETdGradientExp,vOmega,w); + gProb = optiProb.LETdBP.GetGradientProb(); + + %Only implemented for first scenario now + weightGradient = weightGradient + gProb{1}; + end +end + +end \ No newline at end of file diff --git a/matRad/optimization/matRad_DoseOptimizationFunction.m b/matRad/optimization/matRad_DoseOptimizationFunction.m index b377c54a3..42bab81a9 100644 --- a/matRad/optimization/matRad_DoseOptimizationFunction.m +++ b/matRad/optimization/matRad_DoseOptimizationFunction.m @@ -56,17 +56,58 @@ % Helper methods methods (Access = public) - function doseParams = getDoseParameters(obj) + function objParams = getDoseParameters(obj) % get only the dose related parameters. - ix = cellfun(@(c) isequal('dose',c),obj.parameterTypes); - doseParams = [obj.parameters{ix}]; + if isempty(obj.parameterTypes) + objParams = []; + end + + if cellfun(@(c) isequal('dose',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('dose',c),obj.parameterTypes); + objParams = [obj.parameters{ix}]; + end + + if cellfun(@(c) isequal('dirtyDose',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('dirtyDose',c),obj.parameterTypes); + objParams = [obj.parameters{ix}]; + end + + if cellfun(@(c) isequal('LETxd',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('LETxd',c),obj.parameterTypes); + objParams = [obj.parameters{ix}]; + end + + if cellfun(@(c) isequal('LETd',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('LETd',c),obj.parameterTypes); + objParams = [obj.parameters{ix}]; + end + end - - function obj = setDoseParameters(obj,doseParams) + + function obj = setDoseParameters(obj,objParams) % set only the dose related parameters. - ix = cellfun(@(c) isequal('dose',c),obj.parameterTypes); - obj.parameters(ix) = num2cell(doseParams); - end + + if cellfun(@(c) isequal('dose',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('dose',c),obj.parameterTypes); + obj.parameters(ix) = num2cell(objParams); + end + + if cellfun(@(c) isequal('dirtyDose',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('dirtyDose',c),obj.parameterTypes); + obj.parameters(ix) = num2cell(objParams); + end + + if cellfun(@(c) isequal('LETxd',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('LETxd',c),obj.parameterTypes); + obj.parameters(ix) = num2cell(objParams); + end + + if cellfun(@(c) isequal('LETd',c),obj.parameterTypes) + ix = cellfun(@(c) isequal('LETd',c),obj.parameterTypes); + obj.parameters(ix) = num2cell(objParams); + end + + end end methods (Access = private) diff --git a/matRad/optimization/matRad_getObjectivesAndConstraints.m b/matRad/optimization/matRad_getObjectivesAndConstraints.m index 9cc9734eb..d43c7b317 100644 --- a/matRad/optimization/matRad_getObjectivesAndConstraints.m +++ b/matRad/optimization/matRad_getObjectivesAndConstraints.m @@ -29,8 +29,12 @@ if strcmp(env,'MATLAB') %use the package methodology for Matlab (stable) mpkgObjectives = meta.package.fromName('DoseObjectives'); + mpkgDirtyObjectives = meta.package.fromName('DirtyDoseObjectives'); + mpkgLETxDoseObjectives = meta.package.fromName('LETxDoseObjectives'); + mpkgLETdObjectives = meta.package.fromName('LETdObjectives'); mpkgConstraints = meta.package.fromName('DoseConstraints'); - classList = [mpkgObjectives.ClassList; mpkgConstraints.ClassList]; + + classList = [mpkgObjectives.ClassList;mpkgDirtyObjectives.ClassList;mpkgLETxDoseObjectives.ClassList;mpkgLETdObjectives.ClassList;mpkgConstraints.ClassList]; classList = classList(not([classList.Abstract])); @@ -46,12 +50,12 @@ end else currDir = fileparts(mfilename('fullpath')); - objectiveDir = [currDir filesep '+DoseObjectives']; + objectiveDir = [currDir filesep '+DoseObjectives';currDir filesep '+DirtyDoseObjectives'; currDir filesep '+LETxDoseObjectives'; currDir filesep '+LETdObjectives']; constraintDir = [currDir filesep '+DoseConstraints']; objFiles = dir(objectiveDir); for i=1:numel(objFiles) - objFiles(i).pkgName = 'DoseObjectives'; + objFiles(i).pkgName = 'DoseObjectives'; end constrFiles = dir(constraintDir); for i=1:numel(objFiles) diff --git a/matRad/optimization/projections/matRad_BackProjection.m b/matRad/optimization/projections/matRad_BackProjection.m index b6299cda8..d62f511c8 100644 --- a/matRad/optimization/projections/matRad_BackProjection.m +++ b/matRad/optimization/projections/matRad_BackProjection.m @@ -77,7 +77,7 @@ dOmegaV = obj.dOmegaV; end - function wGrad = GetGradient(obj) + function wGrad = GetGradient(obj) wGrad = obj.wGrad; end @@ -105,7 +105,7 @@ function wGrad = projectGradient(obj,dij,doseGrad,w) wGrad = cell(size(dij.physicalDose)); wGrad(obj.scenarios) = arrayfun(@(scen) projectSingleScenarioGradient(obj,dij,doseGrad,scen,w),obj.scenarios,'UniformOutput',false); - end + end function wGrad = projectGradientProb(obj,dij,dExpGrad,dOmegaVgrad,w) wGrad = cell(size(dij.physicalDose)); diff --git a/matRad/optimization/projections/matRad_DirtyDoseProjection.m b/matRad/optimization/projections/matRad_DirtyDoseProjection.m new file mode 100644 index 000000000..647bf1901 --- /dev/null +++ b/matRad/optimization/projections/matRad_DirtyDoseProjection.m @@ -0,0 +1,70 @@ +classdef matRad_DirtyDoseProjection < matRad_BackProjection +% matRad_DirtyDoseProjection class to compute dirty dose during optimization +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2019 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + methods + function obj = matRad_DirtyDoseProjection() + + end + end + + methods + function dD = computeSingleScenario(~,dij,scen,w) + if ~isempty(dij.dirtyDose{scen}) + dD = dij.dirtyDose{scen}*w; + else + dD = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + + function [dDExp,dDOmegaV] = computeSingleScenarioProb(~,dij,scen,w) + if ~isempty(dij.dirtyDoseExp{scen}) + dDExp = dij.dirtyDoseExp{scen}*w; + + for i = 1:size(dij.dirtyDoseOmega,2) + dDOmegaV{scen,i} = dij.dirtyDoseOmega{scen,i} * w; + end + else + dDExp = []; + dDOmegaV = []; + end + end + + function wGrad = projectSingleScenarioGradient(~,dij,dirtyDoseGrad,scen,~) + if ~isempty(dij.dirtyDose{scen}) + wGrad = (dirtyDoseGrad{scen}' * dij.dirtyDose{scen})'; + else + wGrad = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + + function wGrad = projectSingleScenarioGradientProb(~,dij,dDExpGrad,dDOmegaVgrad,scen,~) + if ~isempty(dij.dirtyDoseExp{scen}) + wGrad = (dDExpGrad{scen}' * dij.dirtyDoseExp{scen})'; + wGrad = wGrad + 2 * dDOmegaVgrad; + else + wGrad = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + end +end + + diff --git a/matRad/optimization/projections/matRad_LETdProjection.m b/matRad/optimization/projections/matRad_LETdProjection.m new file mode 100644 index 000000000..340ce1648 --- /dev/null +++ b/matRad/optimization/projections/matRad_LETdProjection.m @@ -0,0 +1,96 @@ +classdef matRad_LETdProjection < matRad_BackProjection +% matRad_LETdProjection class to compute physical dose during optimization +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2019 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + methods + function obj = matRad_LETdProjection() + + end + end + + methods + function LETd = computeSingleScenario(~,dij,scen,w) + + if ~isempty(dij.mLETDose{scen}) + + d = dij.physicalDose{scen}*w; + + % Computation of the LETd + + LETD = dij.mLETDose{scen} * w; % computes the nominator quickly (should be correct) + LETD(d > 0) = LETD(d > 0)./d(d > 0); %dose averging -> avoid div by zero + + LETd = LETD; + + else + LETd = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + + end + + function [LdExp,LdOmegaV] = computeSingleScenarioProb(~,dij,scen,w) + if ~isempty(dij.LETdExp{scen}) + LdExp = dij.LETdExp{scen}*w; + + for i = 1:size(dij.LETdOmega,2) + LdOmegaV{scen,i} = dij.LETdOmega{scen,i} * w; + end + else + LdExp = []; + LdOmegaV = []; + end + end + + function wGrad = projectSingleScenarioGradient(~,dij,doseGrad,scen,w) + if ~isempty(dij.mLETDose{scen}) + + % LETd = this.computeSingleScenario(dij,scen,w); + + d = dij.physicalDose{scen} * w; + mLD = dij.mLETDose{scen} * w; + + % doseGrad * u'v/v^2 + vox_tmp = zeros(dij.doseGrid.numOfVoxels,1); + vox_tmp(d>0) = (doseGrad{scen}(d>0)./d(d>0)); + firstDerivativeterm = (dij.mLETDose{scen}' * vox_tmp); + + vox_tmp = zeros(dij.doseGrid.numOfVoxels,1); + vox_tmp(d>0) = (mLD(d>0)./(d(d>0).^2)).* doseGrad{scen}(d>0); + secondDerivativeterm = (vox_tmp' * dij.physicalDose{scen})'; + + wGrad = (firstDerivativeterm - secondDerivativeterm); + + else + wGrad = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + + end + + function wGrad = projectSingleScenarioGradientProb(~,dij,LdExpGrad,LdOmegaVgrad,scen,~) + if ~isempty(dij.LETdExp{scen}) + wGrad = (LdExpGrad{scen}' * dij.LETdExp{scen})'; + wGrad = wGrad + 2 * LdOmegaVgrad; + else + wGrad = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + end +end \ No newline at end of file diff --git a/matRad/optimization/projections/matRad_LETxDoseProjection.m b/matRad/optimization/projections/matRad_LETxDoseProjection.m new file mode 100644 index 000000000..3f5066b00 --- /dev/null +++ b/matRad/optimization/projections/matRad_LETxDoseProjection.m @@ -0,0 +1,70 @@ +classdef matRad_LETxDoseProjection < matRad_BackProjection +% matRad_LETxProjection class to compute physical dose during optimization +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2019 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + methods + function obj = matRad_LETxDoseProjection() + + end + end + + methods + function LxD = computeSingleScenario(~,dij,scen,w) + if ~isempty(dij.mLETDose{scen}) + LxD = dij.mLETDose{scen}*w; + else + LxD = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + + function [LxDExp,LxDOmegaV] = computeSingleScenarioProb(~,dij,scen,w) + if ~isempty(dij.mLETDoseExp{scen}) + LxDExp = dij.mLETDoseExp{scen}*w; + + for i = 1:size(dij.mLETDoseOmega,2) + LxDOmegaV{scen,i} = dij.mLETDoseOmega{scen,i} * w; + end + else + LxDExp = []; + LxDOmegaV = []; + end + end + + function wGrad = projectSingleScenarioGradient(~,dij,mLETDoseGrad,scen,~) + if ~isempty(dij.mLETDose{scen}) + wGrad = (mLETDoseGrad{scen}' * dij.mLETDose{scen})'; + else + wGrad = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + + function wGrad = projectSingleScenarioGradientProb(~,dij,LxDExpGrad,LxDOmegaVgrad,scen,~) + if ~isempty(dij.mLETDoseExp{scen}) + wGrad = (LxDExpGrad{scen}' * dij.mLETDoseExp{scen})'; + wGrad = wGrad + 2 * LxDOmegaVgrad; + else + wGrad = []; + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n'); + end + end + end +end + + diff --git a/matRad/plotting/matRad_plotDirtynCleanDoseShare.m b/matRad/plotting/matRad_plotDirtynCleanDoseShare.m new file mode 100644 index 000000000..b2bb6d4c5 --- /dev/null +++ b/matRad/plotting/matRad_plotDirtynCleanDoseShare.m @@ -0,0 +1,72 @@ +function [stackedbarDose,physDoseInDepth,RBExDcurve,RBEcurve,h,d] = matRad_plotDirtynCleanDoseShare(definedEnd,index,ct,ctCube,dij,resultGUI,add,LET_thres,k,displayComparison) +% +% call +% [LETbeamletSpectrum, PhysDose_LET] = matRad_DoseComparison(definedEnd,index,ct,ctCube,dij,resultGUI,add,LET_thres,k) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_DoseComparison(definedEnd,index,ct,ctCube,dij,resultGUI,add,LET_thres,k,displayComparison) +% +% input +% definedEnd: last voxel to consider +% index: voxel coordinates from the cube index in matRadGUI in [y x z] +% ct: matRad ct struct +% ctCube: matRad ct.cube in ct struct +% dij: matRad dij struct +% resultGUI: matRad resultGUI struct +% add: add [y x z] to your index: 0 for no change, 1 for one step and so on +% LET_thres: LET threshold: above this = dirty dose, below this = clean dose +% k: k the dimension that should change +% displayComparisony: displays a comparison of dose and RBE +% +% output +% stackedbarDose: bar plot with stacked dose shares adding up to the total dose +% physDoseInDepth: total physical dose plotted in a line depending on the penetration depth +% RBExDcurve: curve to show the RBExD distribution of different voxels +% RBEcurve: curve to show the RBE distribution of different voxels +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Initialize arrays + h = []; + l = []; + t = []; + d = []; + RBE = []; + RBExD = []; + + for i = 1:definedEnd % loop to describe a line of voxels within a certain area defined by the index + + % calculating dirty and clean dose for each voxel; don't display plots of matRad_returnDirtyandCleanDose, use [] + [highLETphysDose,lowLETphysDose,totalphysDose] = matRad_returnDirtyandCleanDose(index,ct,ctCube,dij,resultGUI,LET_thres,[],[],[]); + + depth = index(k); + + % safe results + h = [h highLETphysDose]; + l = [l lowLETphysDose]; + t = [t totalphysDose]; + d = [d depth]; + + % calculating RBExD and RBE for each voxel and storing them into arrays + rowRBExD = resultGUI.RBExD(index(1),index(2),index(3)); + rowRBE = resultGUI.RBE(index(1),index(2),index(3)); + + RBExD = [RBExD rowRBExD]; + RBE = [RBE rowRBE]; + + index = index + add; % describes next voxel + end + + if ~exist('displayComparison','var') || isempty(displayComparison) % for showing image, displayComparison must be 1 + displayComparison = []; + elseif displayComparison == 1 + sumDose = [h' l']; % stores the dirty and clean dose into one array + figure; + hold on + stackedbarDose = bar(d,sumDose,"stacked",'EdgeColor',[0 1 0],'LineWidth',1); % plots stacked bars of dirty and clean dose to see the respective share of the total dose + physDoseInDepth = plot(d,t,"Color",[0 0 1],'LineWidth',1); % plots the total physical dose + RBExDcurve = plot(d,RBExD, "Color", [1 0 0], 'LineWidth', 1); % plots the RBExD + RBEcurve = plot(d,RBE, "Color", [1 0 1], 'LineWidth', 1); % plots the RBE + xlabel('depth in mm'); ylabel('physical dose in Gy'); + legend('dirty share of physical dose','clean share of physical dose','total physical dose', 'RBExD','RBE',Location='best'); + end + +end \ No newline at end of file diff --git a/matRad/plotting/matRad_plotLETbeamletSpectrumInVoxel.m b/matRad/plotting/matRad_plotLETbeamletSpectrumInVoxel.m new file mode 100644 index 000000000..40eb70fb6 --- /dev/null +++ b/matRad/plotting/matRad_plotLETbeamletSpectrumInVoxel.m @@ -0,0 +1,53 @@ +function [wPhysDose,LET] = matRad_plotLETbeamletSpectrumInVoxel(index, ct, ctCube, dij, resultGUI, bins, displayfigures) +% matRad histogram plot for LET beamlet Spectrum for a certain Voxel & +% bivariate histogram plot to see the dose distribution and associated LET +% +% call +% [LETbeamletSpectrum, PhysDose_LET] = matRad_plotLETbeamletSpectrumInVoxel(index,ct,ctCube,dij,resultGUI) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_plotLETbeamletSpectrumInVoxel(index,ct,ctCube,dij,resultGUI,displayfigures) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_plotLETbeamletSpectrumInVoxel(index,ct,ctCube,dij,resultGUI,bins) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_plotLETbeamletSpectrumInVoxel(index,ct,ctCube,dij,resultGUI,bins,displayfigures) +% +% input +% index: voxel coordinates from the cube index in matRadGUI in [y x z] +% ct: matRad ct struct +% ctCube: matRad ct.cube in ct struct +% dij: matRad dij struct +% resultGUI: matRad resultGUI struct +% bins: (optional) specifies the number of bins +% displayfigures: (optional) displays the figures +% +% output +% wPhysDose: row of physicalDose multipled with w vector +% LET: LET in a certain Voxel from every bixel +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + ind = sub2ind(size(ct.cube{ctCube}),index(1),index(2),index(3)); % choose one voxel + rowOfmLETDose = full(dij.mLETDose{ctCube}(ind,:)); % returns row in the mLETDose matrix for the voxel + rowDijPhysicalDose = full(dij.physicalDose{ctCube}(ind,:)); % returns row in the physicalDose matrix for the voxel + + wLETDose = rowOfmLETDose.*resultGUI.w'; % calculates weighted dose averaged LET + wPhysDose = rowDijPhysicalDose.*resultGUI.w'; % calculates weighted physical dose + + LET = wLETDose./wPhysDose; % calculates pure LET + + if ~exist('displayfigures','var') || isempty(displayfigures) % only if displayfigures is 1 it will return the LET beamlet Spectrum as histogram + displayfigures = []; + elseif displayfigures == 1 + figure + LETbeamletSpectrum = histogram(LET); + xlabel('LET'); ylabel('number of bixel') + + if ~exist('bins', 'var') || isempty(bins) % here the number of bins is defined + bins = []; + else + Nbins = morebins(LETbeamletSpectrum); + LETbeamletSpectrum.NumBins = bins; + end + + figure + PhysDose_LET = histogram2(LET,wPhysDose); % bivariate histogram to see the dose distribution and associated LET + xlabel('LET'); ylabel('physical Dose'); zlabel('number of bixel') + end +end diff --git a/matRad/plotting/matRad_returnDirtyandCleanDose.m b/matRad/plotting/matRad_returnDirtyandCleanDose.m new file mode 100644 index 000000000..981ff1ad3 --- /dev/null +++ b/matRad/plotting/matRad_returnDirtyandCleanDose.m @@ -0,0 +1,83 @@ +function [highLETphysDose,lowLETphysDose,totalphysDose] = matRad_returnDirtyandCleanDose(index,ct,ctCube,dij,resultGUI,LET_thres,displayfigures,maxDirtyDose,bins) +% matRad function to return histogram plots for the dirty and clean dose share of physical dose +% +% call +% [LETbeamletSpectrum, PhysDose_LET] = matRad_returnDirtyandCleanDose(index,ct,ctCube,dij,resultGUI,LET_thres) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_returnDirtyandCleanDose(index,ct,ctCube,dij,resultGUI,LET_thres,displayfigures) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_returnDirtyandCleanDose(index,ct,ctCube,dij,resultGUI,LET_thres,displayfigures,maxDirtyDose) +% [LETbeamletSpectrum, PhysDose_LET] = matRad_returnDirtyandCleanDose(index,ct,ctCube,dij,resultGUI,LET_thres,displayfigures,maxDirtyDose,bins) +% +% input +% index: voxel coordinates from the cube index in matRadGUI in [y x z] +% ct: matRad ct struct +% ctCube: matRad ct.cube in ct struct +% dij: matRad dij struct +% resultGUI: matRad resultGUI struct +% LET_thres: LET threshold: above this = dirty dose, below this = clean dose +% displayfigures: (optional) displays the figures +% maxDirtyDose: (optional) maximum allowed dirty dose, dirty dose in a voxel exceeding this value gets penalized +% bins: (optional) specifies the number of bins +% +% output +% shareDose: bar plot for comparing dirty, clean and the total dose; each summed up to one value +% dirtyDose: histogram plot for the dirty dose share of physical dose +% cleanDose: histogram plot for the clean dose share of physical dose +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + [wPhysDose,LET] = matRad_plotLETbeamletSpectrumInVoxel(index, ct, ctCube, dij, resultGUI,[],[]); + + highLET = LET > LET_thres; % defining high LET above a certain threshold + + % share of high LET contribution in physical dose + highLETphysDose = sum(highLET.*wPhysDose); % sum of all the physical dose with LET value higher than the threshold + lowLETphysDose = sum(~highLET.*wPhysDose); % sum of all the physical dose with LET value lower than the threshold + totalphysDose = sum(wPhysDose); % sum of the total dose deposit in this voxel + + if ~exist('displayfigures','var') || isempty(displayfigures) % only if displayfigures is 1, it will return a bar plot + displayfigures = []; + elseif displayfigures == 1 + figure + dose = [highLETphysDose lowLETphysDose totalphysDose]; + C2 = {'dirty share of physical dose';'clean share of physical dose';'total physical dose in this voxel'}; + shareDose = bar(dose); % summary of the composition of physical dose + xticklabels(C2) + end + + dirtydose = wPhysDose(highLET); % dirty dose defined as the physical dose deposited by the bixels with high LET + cleandose = wPhysDose(~highLET); % clean dose defined as the physical dose deposited by the bixels with low LET + + if ~exist('maxDirtyDose', 'var') || isempty(maxDirtyDose) % if there is a maximum dose up to which the dirty dose should be considered, remove all the dose values above + maxDirtyDose = []; + else + dirtydose(dirtydose>maxDirtyDose) = []; + end + + if ~exist('displayfigures','var') || isempty(displayfigures) + displayfigures = []; + elseif displayfigures == 1 + figure + dirtyDose = histogram(dirtydose); % plot the dirty dose as histogram + + if ~exist('bins', 'var') || isempty(bins) + bins = []; + else + Nbins = morebins(dirtyDose); % define the number of bins + dirtyDose.NumBins = bins; + end + xlabel('dirty dose'); ylabel('number of bixel') + + figure + cleanDose = histogram(cleandose); % plot the clean dose as histogram + + if ~exist('bins', 'var') || isempty(bins) + bins = []; + else + Nbins = morebins(cleanDose); % define the number of bins + cleanDose.NumBins = bins; + end + xlabel('clean dose'); ylabel('number of bixel') + end + +end \ No newline at end of file diff --git a/matRad/util/matRad_calcCubes.m b/matRad/util/matRad_calcCubes.m index d63b5e12a..7f0a4526b 100644 --- a/matRad/util/matRad_calcCubes.m +++ b/matRad/util/matRad_calcCubes.m @@ -7,9 +7,9 @@ % resultGUI = matRad_calcCubes(w,dij,scenNum) % % input -% w: bixel weight vector -% dij: dose influence matrix -% scenNum: optional: number of scenario to calculated (default 1) +% w: bixel weight vector +% dij: dose influence matrix +% scenNum: optional: number of scenario to calculated (default 1) % % output % resultGUI: matRad result struct @@ -75,18 +75,33 @@ end end - %% LET -% consider LET if isfield(dij,'mLETDose') + %Apply a threshold to avoid division by doses close to zero + if ~isfield(dij,'doseLETthreshold') + doseLETthreshold = 0.01 * resultGUI.physicalDose; + else + doseLETthreshold = dij.doseLETthreshold; + end + + % compute LETd and LETxDose for all beams individually and together for i = 1:length(beamInfo) - LETDoseCube = reshape(full(dij.mLETDose{scenNum} * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); - resultGUI.(['LET', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); - ix = resultGUI.(['physicalDose', beamInfo(i).suffix]) > 0; - resultGUI.(['LET', beamInfo(i).suffix])(ix) = LETDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix); + resultGUI.(['LETxDose', beamInfo(i).suffix]) = reshape(full(dij.mLETDose{scenNum}*(resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); + resultGUI.(['LETd', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + ix = resultGUI.(['physicalDose', beamInfo(i).suffix]) > doseLETthreshold; + resultGUI.(['LETd', beamInfo(i).suffix])(ix) = resultGUI.(['LETxDose', beamInfo(i).suffix])(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix); end end +%% dirty & clean dose +if isfield(dij,'dirtyDoseThreshold') + for i = 1:length(beamInfo) + % consider dirty dose + resultGUI.(['dirtyDose', beamInfo(i).suffix]) = reshape(full(dij.dirtyDose{1}*(resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); + % consider clean dose + resultGUI.(['cleanDose', beamInfo(i).suffix]) = reshape(full(dij.cleanDose{1}*(resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); + end +end %% RBE weighted dose % consider RBE for protons and skip varRBE calculation diff --git a/test/autoExampleTest/test_examples.m b/test/autoExampleTest/test_examples.m index d40dbf42e..1b03adb5c 100644 --- a/test/autoExampleTest/test_examples.m +++ b/test/autoExampleTest/test_examples.m @@ -21,6 +21,9 @@ 'examples/matRad_example10_4DphotonRobust.m',... 'examples/matRad_example11_helium.m',... 'examples/matRad_example12_simpleParticleMonteCarlo.m',... + 'examples/matRad_example16_LETd.m',... + 'examples/matRad_example17_LETxDose.m',... + 'examples/matRad_example18_DirtyDose.m',... 'matRad.m',... };