From f2a9687648dab4ebb4f4f693c1f8c9dec2e600d2 Mon Sep 17 00:00:00 2001 From: remocristoforetti Date: Fri, 30 Jun 2023 11:47:22 +0200 Subject: [PATCH 1/2] Changes to fix effect prescription in optimizationProblem --- matRad_fluenceOptimizationJO.m | 1 + .../matRad_OptimizationProblemMixMod.m | 1 + .../matRad_constraintFunctions.m | 6 ++++++ .../matRad_constraintJacobian.m | 6 ++++++ .../matRad_getConstraintBounds.m | 10 ++++++++-- .../matRad_objectiveFunction.m | 5 +++++ .../matRad_objectiveGradient.m | 8 +++++++- 7 files changed, 34 insertions(+), 3 deletions(-) diff --git a/matRad_fluenceOptimizationJO.m b/matRad_fluenceOptimizationJO.m index a21146d3f..0770bceb6 100644 --- 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 index a9221fb66..5f7fa2ee5 100644 --- 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 index c355e1da7..811b6afa0 100644 --- 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 index a411cc5a2..aa3bcbe42 100644 --- 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 index 000516e02..21dc77fdc 100644 --- a/optimization/@matRad_OptimizationProblemMixMod/matRad_getConstraintBounds.m +++ b/optimization/@matRad_OptimizationProblemMixMod/matRad_getConstraintBounds.m @@ -63,14 +63,20 @@ [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}))]; @@ -78,7 +84,7 @@ %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 index 8c8020c39..358bf9d30 100644 --- 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 index 35e9b2a16..ad4b63365 100644 --- 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 @@ -346,7 +352,7 @@ weightGradient = weightGradient + gProb{1}; end % code snippet to check the gradient - gradientChecker = 1; + gradientChecker = 0; if gradientChecker == 1 f = matRad_objectiveFunction(optiProb,w,dij,cst); epsilon = 1e-6; From 3613d60a9e65a124904c1ee51c8a21c3ef9c865c Mon Sep 17 00:00:00 2001 From: remocristoforetti Date: Fri, 30 Jun 2023 11:47:44 +0200 Subject: [PATCH 2/2] Add visualization for checking --- matRad_MixedModality.m | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/matRad_MixedModality.m b/matRad_MixedModality.m index 9e6963608..4709be303 100644 --- 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'); \ No newline at end of file