e0404 · JenHardt · Dec 3, 2021 · Dec 3, 2021 · Dec 3, 2021 · Dec 3, 2021
diff --git a/matRad/bioModels/matRad_BiologicalModel.m b/matRad/bioModels/matRad_BiologicalModel.m
@@ -12,7 +12,8 @@
     % input
     %   sRadiationMode:     radiation modality 'photons' 'protons' 'carbon'
     %   sQuntityOpt:        string to denote the quantity used for
-    %                       optimization 'physicalDose', 'RBExD', 'effect'
+    %                       optimization 'physicalDose', 'RBExD',
+    %                       'effect','BED'
     %   sModel:             string to denote which biological model is used
     %                       '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 
@@ -61,7 +62,7 @@
     properties(Constant = true)
         AvailableModels                 = {'none','constRBE','MCN','WED','LEM','HEL'};   % cell array determines available models - if cell is deleted then the corersponding model can not be generated
         AvailableradiationModealities   = {'photons','protons','helium','carbon'};
-        AvailableQuantitiesForOpt       = {'physicalDose','effect','RBExD'};
+        AvailableQuantitiesForOpt       = {'physicalDose','effect','RBExD','BED'};
 
         AvailableAlphaXBetaX = {[0.036 0.024],    'prostate';
             [0.089 0.287],    'rectum and normal tissue';
@@ -161,6 +162,15 @@
                                 else
                                     setDefaultValues = true;   
                                 end
+
+                            case {'BED'}
+                                if strcmp( this.model,'none')
+                                    boolCHECK        = true;
+                                    this.bioOpt      = true;
+                                    this.quantityVis = 'RBExD';
+                                else
+                                    setDefaultValues = true;   
+                                end
 
                             otherwise
                                 matRad_cfg.dispWarning('matRad: Invalid biological optimization quantity: %s; using physical dose instead.',this.quantityOpt);
@@ -211,6 +221,16 @@
                                     matRad_cfg.dispWarning(['matRad: Invalid biological model: ' this.model  '; using MCN Model instead.']);
                                     this.model  = 'MCN';
                                 end
+
+                            case {'BED'}
+                                if sum(strcmp( this.model, {'MCN','WED'})) == 1
+                                    boolCHECK           = true;
+                                    this.bioOpt         = true;
+                                    this.quantityVis    = 'RBExD';
+                                else
+                                    matRad_cfg.dispWarning(['matRad: Invalid biological model: ' this.model  '; using MCN Model instead.']);
+                                    this.model  = 'MCN';
+                                end
 
                             otherwise
                                 matRad_cfg.dispWarning(['matRad: Invalid biological optimization quantity: ' this.quantityOpt  '; using "none" instead.']);
@@ -230,7 +250,7 @@
                                     this.model  = 'none';
                                 end
 
-                            case {'effect','RBExD'}
+                            case {'effect','RBExD','BED'}
                                 if sum(strcmp(this.model,{'LEM','HEL'})) > 0
                                     boolCHECK           = true;
                                     this.bioOpt         = true;
@@ -240,6 +260,7 @@
                                     this.model = 'none';
                                 end
 
+
                             otherwise
                                 matRad_cfg.dispWarning(['matRad: Invalid biological optimization quantity: ' this.quantityOpt  '; using "none" instead.']);
                                 this.quantityOpt = 'physicalDose';
@@ -258,7 +279,7 @@
                                     this.model  = 'none';
                                 end
 
-                            case {'effect','RBExD'}
+                            case {'effect','RBExD','BED'}
                                 if strcmp(this.model,'LEM')
                                     boolCHECK           = true;
                                     this.bioOpt         = true;
@@ -282,7 +303,7 @@
 
             % check quantity for optimization
             if this.bioOpt 
-                if sum(strcmp(this.quantityOpt,{'physicalDose','RBExD','effect'})) == 0
+                if sum(strcmp(this.quantityOpt,this.AvailableQuantitiesForOpt)) == 0
                     matRad_cfg.dispError(['matRad: Invalid quantity for optimization: ' this.quantityOpt  ]);
                 end
             else
@@ -294,7 +315,7 @@
 
             % check quantity for visualization
             if this.bioOpt 
-                if sum(strcmp(this.quantityVis,{'physicalDose','RBExD','effect'})) == 0
+                if sum(strcmp(this.quantityVis,this.AvailableQuantitiesForOpt)) == 0
                     matRad_cfg.dispError(['matRad: Invalid quantity for visualization: ' this.quantityVis  ]);
                 end
             else

diff --git a/matRad/bioModels/matRad_bioModel.m b/matRad/bioModels/matRad_bioModel.m
@@ -13,7 +13,8 @@
 % input
 %   sRadiationMode:     radiation modality 'photons' 'protons' 'carbon'
 %   sQuantityOpt:       string to denote the quantity used for
-%                       optimization 'physicalDose', 'RBExD', 'effect'
+%                       optimization 'physicalDose', 'RBExD', 'effect',
+%                       'BED'
 %   sModel:             string to denote which biological model is used
 %                       '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

diff --git a/matRad/gui/widgets/matRad_PlanWidget.m b/matRad/gui/widgets/matRad_PlanWidget.m
@@ -30,6 +30,7 @@
 
     properties (Constant)
         Modalities = {'photons','protons','carbon', 'helium'};
+        availableProjections = {  'physicalDose'; 'RBExD'; 'effect'; 'BED'; }
     end
 
     methods
@@ -347,7 +348,7 @@
             h33 = uicontrol(...
                 'Parent',h12,...
                 'Units','normalized',...
-                'String',{  'physicalDose'; 'RBExD'; 'effect';  },...
+                'String',this.availableProjections,...
                 'TooltipString',txt,...
                 'Style','popupmenu',...
                 'Value',1,...

diff --git a/matRad/matRad_fluenceOptimization.m b/matRad/matRad_fluenceOptimization.m
@@ -189,8 +189,26 @@
         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
 
+    elseif strcmp(pln.bioParam.quantityOpt, 'BED')
+
+        if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose')
+            abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX;
+            meanBED = mean((dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^2)./cst{ixTarget,5}.alphaX);
+            BEDTarget = doseTarget.*(1 + doseTarget./abr);
+        elseif isfield(dij, 'RBE')
+            abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX;
+            meanBED = mean(dij.RBE.*dij.physicalDose{1}(V,:)*wOnes.*(1+dij.RBE.*dij.physicalDose{1}(V,:)*wOnes./abr));
+            BEDTarget = dij.RBE.*doseTarget.*(1 + dij.RBE.*doseTarget./abr);
+        else
+            abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX;
+            meanBED = mean(dij.physicalDose{1}(V,:)*wOnes.*(1+dij.physicalDose{1}(V,:)*wOnes./abr));
+            BEDTarget = doseTarget.*(1 + doseTarget./abr);
+        end
+
+        bixelWeight =  BEDTarget/meanBED;
+        wInit       = wOnes * bixelWeight;
+    end
     matRad_cfg.dispInfo('chosen weights adapted to biological dose calculation!\n');
 else
     doseTmp = dij.physicalDose{1}*wOnes;
@@ -257,6 +275,9 @@
         else
             backProjection = matRad_VariableRBEProjection;
         end
+    case 'BED'
+        backProjection = matRad_BEDProjection; 
+
     case 'physicalDose'
         backProjection = matRad_DoseProjection;
     otherwise

diff --git a/matRad/optimization/+DoseConstraints/matRad_MinMaxDVH.m b/matRad/optimization/+DoseConstraints/matRad_MinMaxDVH.m
@@ -29,6 +29,7 @@
         voxelScalingRatio = 1;
         referenceScalingVal = 0.01;
         parameters = {30,0,100};
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseConstraints/matRad_MinMaxDose.m b/matRad/optimization/+DoseConstraints/matRad_MinMaxDose.m
@@ -27,6 +27,7 @@
     properties
         parameters = {0,30,1};
         epsilon = 1e-3; %slack parameter for the logistic approximation
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseConstraints/matRad_MinMaxEUD.m b/matRad/optimization/+DoseConstraints/matRad_MinMaxEUD.m
@@ -27,6 +27,7 @@
 
     properties
         parameters = {5,0,30};
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseConstraints/matRad_MinMaxMeanDose.m b/matRad/optimization/+DoseConstraints/matRad_MinMaxMeanDose.m
@@ -27,6 +27,7 @@
 
     properties
         parameters = {0,30};
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseObjectives/matRad_EUD.m b/matRad/optimization/+DoseObjectives/matRad_EUD.m
@@ -27,6 +27,7 @@
     properties
         parameters = {0, 3.5};
         penalty = 1;
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m b/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m
@@ -27,6 +27,7 @@
     properties
         parameters = {30,95};
         penalty = 1;
+        numOfFractions = NaN;
     end
 
     methods 

diff --git a/matRad/optimization/+DoseObjectives/matRad_MeanDose.m b/matRad/optimization/+DoseObjectives/matRad_MeanDose.m
@@ -27,6 +27,7 @@
     properties
         parameters = {0,1};        
         penalty = 1;
+        numOfFractions = NaN;
     end
 
     methods 

diff --git a/matRad/optimization/+DoseObjectives/matRad_MinDVH.m b/matRad/optimization/+DoseObjectives/matRad_MinDVH.m
@@ -27,6 +27,7 @@
     properties
         parameters = {60,95};
         penalty = 1;
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredDeviation.m b/matRad/optimization/+DoseObjectives/matRad_SquaredDeviation.m
@@ -29,6 +29,7 @@
     properties
         parameters = {60};
         penalty = 1;
+        numOfFractions = NaN;
     end
 
 

diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredOverdosing.m b/matRad/optimization/+DoseObjectives/matRad_SquaredOverdosing.m
@@ -27,6 +27,7 @@
     properties
         parameters = {30};
         penalty = 1;
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m b/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m
@@ -27,6 +27,7 @@
     properties
         parameters = {60};
         penalty = 1;
+        numOfFractions = NaN;
     end
 
     methods

diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m
@@ -30,6 +30,7 @@
 
 BPtype = class(optiProb.BP);
 isEffectBP = strcmp(BPtype,'matRad_EffectProjection');
+isBEDBP = strcmp(BPtype,'matRad_BEDProjection');
 
 % Initialize bounds
 cl = [];
@@ -66,13 +67,20 @@
             if isa(optiFunc,'DoseConstraints.matRad_DoseConstraint')
 
                 if isEffectBP
-
-                    doses  = optiFunc.getDoseParameters();
+
+                    doses = optiFunc.getDoseParameters();
+
                     effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2;
-
                     optiFunc = optiFunc.setDoseParameters(effect);
                 end
-
+
+                if isBEDBP
+                    doses = optiFunc.getDoseParameters();
+                    abr = cst{i,5}.alphaX ./ cst{i,5}.betaX;
+                    BED = optiProb.BP.numOfFractions.*doses.*(1+doses./abr);
+                    optiFunc = optiFunc.setDoseParameters(BED);
+                end
+
                  cl = [cl;optiFunc.lowerBounds(numel(cst{i,4}{1}))];
                  cu = [cu;optiFunc.upperBounds(numel(cst{i,4}{1}))];
 

diff --git a/matRad/optimization/projections/matRad_BEDProjection.m b/matRad/optimization/projections/matRad_BEDProjection.m
@@ -0,0 +1,27 @@
+classdef matRad_BEDProjection < matRad_EffectProjection
+% matRad_BEDProjection class for BED-based optimization
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
+    methods
+        function obj = matRad_BEDProjection()
+        end
+    end
+
+    methods
+        function BED = computeSingleScenario(obj,dij,scen,w)
+                effect = computeSingleScenario@matRad_EffectProjection(obj,dij,scen,w);
+                BED = zeros(dij.doseGrid.numOfVoxels,1);
+                ix = ~(dij.ax{scen}==0);
+                BED(ix) = effect(ix)./ dij.ax{scen}(ix);
+                % photon equivalent BED = n * effect / alphax
+            end 
+
+        function wGrad = projectSingleScenarioGradient(obj,dij,doseGrad,scen,w)
+                ix = ~(dij.ax{scen}==0);
+                doseGradtmp{scen} = zeros(size(doseGrad{scen}));
+                doseGradtmp{scen}(ix) = doseGrad{scen}(ix)./dij.ax{scen}(ix);
+                wGradEffect = projectSingleScenarioGradient@matRad_EffectProjection(obj,dij,doseGradtmp,scen,w); 
+                wGrad = wGradEffect;
+            end
+        end
+end
+
diff --git a/matRad/optimization/projections/matRad_EffectProjection.m b/matRad/optimization/projections/matRad_EffectProjection.m
@@ -21,25 +21,42 @@
 
     methods
         function effect = computeSingleScenario(~,dij,scen,w)
-            if isempty(dij.mAlphaDose{scen}) || isempty(dij.mSqrtBetaDose{scen})
-                effect = [];
-                matRad_cfg = MatRad_Config.instance();
-                matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n');
-            else
+            if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose')
+                if isempty(dij.mAlphaDose{scen}) || isempty(dij.mSqrtBetaDose{scen})
+                    effect = [];
+                    matRad_cfg = MatRad_Config.instance();
+                    matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n');
+                else
                 effect = dij.mAlphaDose{scen}*w + (dij.mSqrtBetaDose{scen}*w).^2;
+                end
+            else 
+                if isfield(dij, 'RBE')
+                    dose = dij.RBE.* dij.physicalDose{scen}*w;
+                else
+                    dose = dij.physicalDose{scen}*w;
+                end
+                effect = dij.ax.*dose + dij.bx.*dose.^2; 
             end 
         end
 
         function wGrad = projectSingleScenarioGradient(~,dij,doseGrad,scen,w)
-            if isempty(dij.mAlphaDose{scen}) || isempty(dij.mSqrtBetaDose{scen})
-                wGrad = [];
-                matRad_cfg = MatRad_Config.instance();
-                matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n');
-            else
+            if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose')
+                if isempty(dij.mAlphaDose{scen}) || isempty(dij.mSqrtBetaDose{scen})
+                    wGrad = [];
+                    matRad_cfg = MatRad_Config.instance();
+                    matRad_cfg.dispWarning('Empty scenario in optimization detected! This should not happen...\n');
+                else
                 vBias = (doseGrad{scen}' * dij.mAlphaDose{scen})';
                 quadTerm = dij.mSqrtBetaDose{scen} * w;
                 mPsi = (2*(doseGrad{scen}.*quadTerm)' * dij.mSqrtBetaDose{scen})';
                 wGrad = vBias + mPsi;
+                end
+            elseif isfield(dij, 'RBE')
+                dose = dij.RBE.* dij.physicalDose{scen}*w;
+                wGrad = ((doseGrad{scen}.*dij.RBE.*(dij.ax + 2*dij.bx.*dose))'*dij.physicalDose{scen})';
+            else 
+                dose = dij.physicalDose{scen}*w;
+                wGrad = ((doseGrad{scen}.*(dij.ax + 2*dij.bx.*dose))'*dij.physicalDose{scen})';
             end
         end