Merge pull request #639 from remocristoforetti/dev_MixModality_effect…

…PrescriptionFix Dev mix modality effect prescription fix
e0404 · Aug 10, 2023 · 0a204f2 · 0a204f2
2 parents da087e3 + 3613d60
commit 0a204f2
Show file tree

Hide file tree

Showing 8 changed files with 58 additions and 6 deletions.
diff --git a/matRad_MixedModality.m b/matRad_MixedModality.m
@@ -29,9 +29,9 @@
 pln(1).propDoseCalc.doseGrid.resolution.y = 8; % [mm]
 pln(1).propDoseCalc.doseGrid.resolution.z = 8; % [mm]
 % pln(1).propDoseCalc.doseGrid.resolution = ct.resolution;
-quantityOpt  = 'physicalDose';     % options: physicalDose, effect, RBExD
+quantityOpt  = 'effect';     % options: physicalDose, effect, RBExD
 %=======================================> Model check error in bioModel
-modelName    = 'none';             % none: for photons, protons, carbon            % constRBE: constant RBE for photons and protons 
+modelName    = 'MCN';             % 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
 
@@ -68,7 +68,7 @@
 pln(2).propDoseCalc.doseGrid.resolution.z = 8; % [mm]
 % pln(2).propDoseCalc.doseGrid.resolution = ct.resolution;
 
-quantityOpt  = 'physicalDose';     % options: physicalDose, effect, RBExD
+quantityOpt  = 'effect';     % options: physicalDose, effect, RBExD
 modelName    = 'none';             % 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
@@ -95,4 +95,25 @@
 % Dij Calculation
 dij = matRad_calcCombiDose(ct,stf,plnJO,cst,false);
 % Fluence optimization 
-resultGUI = matRad_fluenceOptimizationJO(dij,cst,plnJO)
+resultGUI = matRad_fluenceOptimizationJO(dij,cst,plnJO);
+
+%% Visualization
+slice = 59;
+
+photon_plan = resultGUI{2};
+proton_plan = resultGUI{1};
+totalPlan = pln(1).numOfFractions.*proton_plan.(quantityOpt) + pln(2).numOfFractions.*photon_plan.(quantityOpt);
+
+f = figure;
+subplot(1,3,1);
+    imagesc(proton_plan.(quantityOpt)(:,:,slice));
+    matRad_plotVoiContourSlice(gca(f), cst,ct.cube, 1, 1,3,slice);
+    title('Proton Plan');
+subplot(1,3,2);
+    imagesc(photon_plan.(quantityOpt)(:,:,slice));
+    matRad_plotVoiContourSlice(gca(f), cst,ct.cube, 1, 1,3,slice);
+    title('Photon Plan');
+subplot(1,3,3);
+    imagesc(totalPlan(:,:,slice));
+    matRad_plotVoiContourSlice(gca(f), cst,ct.cube, 1, 1,3,slice);
+    title('Total Plan');
diff --git a/matRad_fluenceOptimizationJO.m b/matRad_fluenceOptimizationJO.m
@@ -316,6 +316,7 @@
 
 if strcmp(pln.radiationMode, 'MixMod')
    optiProb = matRad_OptimizationProblemMixMod(backProjection);
+   optiProb.nFractions = pln.numOfFractions;
 else
 
    optiProb = matRad_OptimizationProblem(backProjection);

diff --git a/optimization/@matRad_OptimizationProblemMixMod/matRad_OptimizationProblemMixMod.m b/optimization/@matRad_OptimizationProblemMixMod/matRad_OptimizationProblemMixMod.m
@@ -30,6 +30,7 @@
 
         minimumW = NaN;
         maximumW = NaN;
+        nFractions;
     end
 
     methods

diff --git a/optimization/@matRad_OptimizationProblemMixMod/matRad_constraintFunctions.m b/optimization/@matRad_OptimizationProblemMixMod/matRad_constraintFunctions.m
@@ -64,8 +64,14 @@
          if isa(constraint,'DoseConstraints.matRad_DoseConstraint')
 
             % rescale dose parameters to biological optimization quantity if required
+            doseParameter = constraint.getDoseParameters();
+            constraint = constraint.setDoseParameters(doseParameter./dij.totalNumOfFractions);
+
             constraint = optiProb.BP.setBiologicalDosePrescriptions(constraint,cst{i,5}.alphaX,cst{i,5}.betaX);
 
+            doseParameter = constraint.getDoseParameters();
+            constraint = constraint.setDoseParameters(doseParameter.*dij.totalNumOfFractions);
+
             % retrieve the robustness type
             robustness = constraint.robustness;
 

diff --git a/optimization/@matRad_OptimizationProblemMixMod/matRad_constraintJacobian.m b/optimization/@matRad_OptimizationProblemMixMod/matRad_constraintJacobian.m
@@ -76,8 +76,14 @@
             robustness = constraint.robustness;
 
             % rescale dose parameters to biological optimization quantity if required
+            doseParameter = constraint.getDoseParameters();
+            constraint = constraint.setDoseParameters(doseParameter./dij.totalNumOfFractions);
+
             constraint = optiProb.BP.setBiologicalDosePrescriptions(constraint,cst{i,5}.alphaX,cst{i,5}.betaX);
 
+            doseParameter = constraint.getDoseParameters();
+            constraint = constraint.setDoseParameters(doseParameter.*dij.totalNumOfFractions); 
+
             switch robustness
 
                case 'none' % if conventional opt: just sum objectiveectives of nominal dose

diff --git a/optimization/@matRad_OptimizationProblemMixMod/matRad_getConstraintBounds.m b/optimization/@matRad_OptimizationProblemMixMod/matRad_getConstraintBounds.m
@@ -63,22 +63,28 @@
                     [clTmp,cuTmp] = matRad_getConstBounds(cst{i,6}(j),param);
 %}
             %if ~isempty(strfind(cst{i,6}{j}.type,'constraint'))
+
+
             if isa(optiFunc,'DoseConstraints.matRad_DoseConstraint')
 
                 if isEffectBP
-
+                    doseParameter = optiFunc.getDoseParameters();
+                    optiFunc = optiFunc.setDoseParameters(doseParameter./optiProb.nFractions);
+
                     doses  = optiFunc.getDoseParameters();
                     effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2;
 
                     optiFunc = optiFunc.setDoseParameters(effect);
+                    doseParameters = optiFunc.getDoseParameters();
+                    optiFunc = optiFunc.setDoseParameters(doseParameters*optiProb.nFractions);
                 end
 
                  cl = [cl;optiFunc.lowerBounds(numel(cst{i,4}{1}))];
                  cu = [cu;optiFunc.upperBounds(numel(cst{i,4}{1}))];
 
                 %end
             end
-
+            
         end % over all objectives of structure
 
     end % if structure not empty and target or oar

diff --git a/optimization/@matRad_OptimizationProblemMixMod/matRad_objectiveFunction.m b/optimization/@matRad_OptimizationProblemMixMod/matRad_objectiveFunction.m
@@ -90,8 +90,13 @@
             if isa(objective,'DoseObjectives.matRad_DoseObjective')
 
                 % rescale dose parameters to biological optimization quantity if required
+                doseParameter = objective.getDoseParameters();
+                objective = objective.setDoseParameters(doseParameter./dij.totalNumOfFractions);
+
                 objective = optiProb.BP.setBiologicalDosePrescriptions(objective,cst{i,5}.alphaX,cst{i,5}.betaX);
 
+                doseParameter = objective.getDoseParameters();
+                objective = objective.setDoseParameters(doseParameter.*dij.totalNumOfFractions);                
                 % retrieve the robustness type
                 robustness = objective.robustness;
 

diff --git a/optimization/@matRad_OptimizationProblemMixMod/matRad_objectiveGradient.m b/optimization/@matRad_OptimizationProblemMixMod/matRad_objectiveGradient.m
@@ -100,7 +100,13 @@
                     robustness = objective.robustness;
 
                     % rescale dose parameters to biological optimization quantity if required
+                    doseParameter = objective.getDoseParameters();
+                    objective = objective.setDoseParameters(doseParameter./dij.totalNumOfFractions);
+
                     objective = optiProb.BP.setBiologicalDosePrescriptions(objective,cst{i,5}.alphaX,cst{i,5}.betaX);
+
+                    doseParameter = objective.getDoseParameters();
+                    objective = objective.setDoseParameters(doseParameter.*dij.totalNumOfFractions);
 
                     switch robustness
                         case 'none' % if conventional opt: just sum objectiveectives of nominal dose
-Original file line number
+Diff line change
@@ Expand Up / @@ -30,6 +30,7 @@ @@
             minimumW = NaN;
             maximumW = NaN;
+            nFractions;
         end
         methods
@@ Expand Down @@