From a97b6bead022b6398724fc0ea509923078dbebef Mon Sep 17 00:00:00 2001
From: Niklas Wahl <niklaswahl2011@gmail.com>
Date: Fri, 17 Sep 2021 17:15:03 +0200
Subject: [PATCH] the mean dose objective has now the option to choose
 optimization of a linear or quadratic difference to a reference value

---
 .../+DoseObjectives/matRad_MeanDose.m         | 67 ++++++++++++++++---
 1 file changed, 56 insertions(+), 11 deletions(-)

diff --git a/optimization/+DoseObjectives/matRad_MeanDose.m b/optimization/+DoseObjectives/matRad_MeanDose.m
index 43a1207d9..90163c417 100644
--- a/optimization/+DoseObjectives/matRad_MeanDose.m
+++ b/optimization/+DoseObjectives/matRad_MeanDose.m
@@ -20,19 +20,17 @@
     
     properties (Constant)
         name = 'Mean Dose';
-        parameterNames = {'d^{ref}'};
-        parameterTypes = {'dose'};
-        %parameterNames = {};
-        %parameterIsDose = [];
+        parameterNames = {'d^{ref}','f_{diff}'}; %When optimizing to a reference, one might consider using a quadratic relationship with a non-linear optimizer
+        parameterTypes = {'dose',{'Linear','Quadratic'}};
     end
     
     properties
-        parameters = {0};        
+        parameters = {0,1};        
         penalty = 1;
     end
     
     methods 
-        function obj = matRad_MeanDose(penalty,dMeanRef)
+        function obj = matRad_MeanDose(penalty,dMeanRef,fDiff)
            
             % if we have a struct in first argument
             if nargin == 1 && isstruct(penalty)
@@ -47,7 +45,22 @@
             obj@DoseObjectives.matRad_DoseObjective(inputStruct);
             
             if ~initFromStruct
-                if nargin == 2 && isscalar(dMeanRef)
+                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 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
 
@@ -60,14 +73,46 @@
         
         %% Calculates the Objective Function value
         function fDose = computeDoseObjectiveFunction(obj,dose)
-            %fDose = obj.penalty * abs(mean(dose(:)) - obj.parameters{1});
-            fDose = obj.penalty * (mean(dose(:)) - obj.parameters{1})^2;
+            switch obj.parameters{2}
+                case 1
+                    fDose = obj.penalty * obj.objectiveLinearDiff(dose);
+                case 2
+                    fDose = obj.penalty * obj.objectiveQuadraticDiff(dose);
+                otherwise
+                    matRad_cfg = MatRad_Config.instance();
+                    matRad_cfg.dispError('Invalid setting for %s in Mean Dose Objective!',obj.parameterNames{2});  
+            end
         end
         
         %% Calculates the Objective Function gradient
         function fDoseGrad   = computeDoseObjectiveGradient(obj,dose)
-            %fDoseGrad = (obj.penalty/numel(dose))*sign(dose(:)-obj.parameters{1});
-            fDoseGrad = obj.penalty*2*(mean(dose(:))-obj.parameters{1}) * ones(size(dose(:)))/numel(dose);
+            switch obj.parameters{2}
+                case 1
+                    fDoseGrad = obj.penalty * obj.gradientLinearDiff(dose);
+                case 2
+                    fDoseGrad = obj.penalty * obj.gradientQuadraticDiff(dose);
+                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 fDose = objectiveQuadraticDiff(obj,dose)
+            fDose = (mean(dose(:)) - obj.parameters{1})^2;
+        end
+
+        function fDoseGrad = gradientQuadraticDiff(obj,dose)
+            fDoseGrad = 2*(mean(dose(:))-obj.parameters{1}) * ones(size(dose(:)))/numel(dose);
+        end
+
+        function fDose = objectiveLinearDiff(obj,dose)
+            fDose = abs(mean(dose(:)) - obj.parameters{1});
+        end
+
+        function fDoseGrad = gradientLinearDiff(obj,dose)
+            fDoseGrad = (1/numel(dose))*sign(dose(:)-obj.parameters{1});
         end
     end