From d342341dde834d7f2d3e4697f93837f4024204a7 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Fri, 3 Dec 2021 13:46:47 +0100 Subject: [PATCH 01/14] added number of fractionations as parameter --- optimization/+DoseConstraints/matRad_MinMaxDVH.m | 1 + optimization/+DoseConstraints/matRad_MinMaxDose.m | 1 + optimization/+DoseConstraints/matRad_MinMaxEUD.m | 1 + optimization/+DoseConstraints/matRad_MinMaxMeanDose.m | 1 + optimization/+DoseObjectives/matRad_EUD.m | 1 + optimization/+DoseObjectives/matRad_MaxDVH.m | 1 + optimization/+DoseObjectives/matRad_MeanDose.m | 1 + optimization/+DoseObjectives/matRad_MinDVH.m | 1 + optimization/+DoseObjectives/matRad_SquaredDeviation.m | 1 + optimization/+DoseObjectives/matRad_SquaredOverdosing.m | 1 + optimization/+DoseObjectives/matRad_SquaredUnderdosing.m | 1 + optimization/projections/matRad_BackProjection.m | 1 + 12 files changed, 12 insertions(+) diff --git a/optimization/+DoseConstraints/matRad_MinMaxDVH.m b/optimization/+DoseConstraints/matRad_MinMaxDVH.m index 7215f4458..85769282f 100644 --- a/optimization/+DoseConstraints/matRad_MinMaxDVH.m +++ b/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/optimization/+DoseConstraints/matRad_MinMaxDose.m b/optimization/+DoseConstraints/matRad_MinMaxDose.m index d75196dd0..004ce2f3f 100644 --- a/optimization/+DoseConstraints/matRad_MinMaxDose.m +++ b/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/optimization/+DoseConstraints/matRad_MinMaxEUD.m b/optimization/+DoseConstraints/matRad_MinMaxEUD.m index 1f6346294..b92abb538 100644 --- a/optimization/+DoseConstraints/matRad_MinMaxEUD.m +++ b/optimization/+DoseConstraints/matRad_MinMaxEUD.m @@ -27,6 +27,7 @@ properties parameters = {5,0,30}; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseConstraints/matRad_MinMaxMeanDose.m b/optimization/+DoseConstraints/matRad_MinMaxMeanDose.m index b9f536782..85b61f367 100644 --- a/optimization/+DoseConstraints/matRad_MinMaxMeanDose.m +++ b/optimization/+DoseConstraints/matRad_MinMaxMeanDose.m @@ -27,6 +27,7 @@ properties parameters = {0,30}; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_EUD.m b/optimization/+DoseObjectives/matRad_EUD.m index 15e73e9cf..17b005d27 100644 --- a/optimization/+DoseObjectives/matRad_EUD.m +++ b/optimization/+DoseObjectives/matRad_EUD.m @@ -27,6 +27,7 @@ properties parameters = {0, 3.5}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_MaxDVH.m b/optimization/+DoseObjectives/matRad_MaxDVH.m index 2d5a2bc8f..eb626dd2a 100644 --- a/optimization/+DoseObjectives/matRad_MaxDVH.m +++ b/optimization/+DoseObjectives/matRad_MaxDVH.m @@ -27,6 +27,7 @@ properties parameters = {30,95}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_MeanDose.m b/optimization/+DoseObjectives/matRad_MeanDose.m index 43a1207d9..19ec08f3d 100644 --- a/optimization/+DoseObjectives/matRad_MeanDose.m +++ b/optimization/+DoseObjectives/matRad_MeanDose.m @@ -29,6 +29,7 @@ properties parameters = {0}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_MinDVH.m b/optimization/+DoseObjectives/matRad_MinDVH.m index 2f49dceaf..249c70a2e 100644 --- a/optimization/+DoseObjectives/matRad_MinDVH.m +++ b/optimization/+DoseObjectives/matRad_MinDVH.m @@ -27,6 +27,7 @@ properties parameters = {60,95}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_SquaredDeviation.m b/optimization/+DoseObjectives/matRad_SquaredDeviation.m index 43fc442f4..e79cb1976 100644 --- a/optimization/+DoseObjectives/matRad_SquaredDeviation.m +++ b/optimization/+DoseObjectives/matRad_SquaredDeviation.m @@ -29,6 +29,7 @@ properties parameters = {60}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_SquaredOverdosing.m b/optimization/+DoseObjectives/matRad_SquaredOverdosing.m index b2e25e40f..7370f0885 100644 --- a/optimization/+DoseObjectives/matRad_SquaredOverdosing.m +++ b/optimization/+DoseObjectives/matRad_SquaredOverdosing.m @@ -27,6 +27,7 @@ properties parameters = {30}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m b/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m index 5c0a84aa7..ceb33dce3 100644 --- a/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m +++ b/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m @@ -27,6 +27,7 @@ properties parameters = {60}; penalty = 1; + numOfFractions = NaN; end methods diff --git a/optimization/projections/matRad_BackProjection.m b/optimization/projections/matRad_BackProjection.m index 7ce56b5eb..0cd113754 100644 --- a/optimization/projections/matRad_BackProjection.m +++ b/optimization/projections/matRad_BackProjection.m @@ -24,6 +24,7 @@ properties dij %reference to matRad dij struct (to enable local changes) + numOfFractions = NaN; end From 9aef7123932ad388f5b360b1085c57524106ccf8 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Fri, 3 Dec 2021 13:47:01 +0100 Subject: [PATCH 02/14] calculation of BED cube --- matRad_calcCubes.m | 35 +++++++++++++++++++++++++++++++++-- 1 file changed, 33 insertions(+), 2 deletions(-) diff --git a/matRad_calcCubes.m b/matRad_calcCubes.m index 5eb0ec253..2d20cd61d 100644 --- a/matRad_calcCubes.m +++ b/matRad_calcCubes.m @@ -1,4 +1,4 @@ -function resultGUI = matRad_calcCubes(w,dij,scenNum) +function resultGUI = matRad_calcCubes(w,dij,pln,scenNum) % matRad computation of all cubes for the resultGUI struct % which is used as result container and for visualization in matRad's GUI % @@ -30,7 +30,7 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if nargin < 3 +if ~exist('sceneNum', 'var') scenNum = 1; end @@ -100,6 +100,37 @@ end end +% CACLCULATE BED (biological effective dose) BED +if exist('pln','var') + % When depth Dependent alpha beta values are calculated in dij calculation + alphax = reshape(dij.ax, dij.doseGrid.dimensions); + ix = ~(alphax == 0); + if isfield(dij,'mAlphaDose') && isfield(dij,'mSqrtBetaDose') + for i = 1:length(beamInfo) + % photon equivaluent BED = n * effect / alphax + resultGUI.(['BED', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + resultGUI.(['BED', beamInfo(i).suffix])(ix) = full(pln.numOfFractions.*resultGUI.(['effect', beamInfo(i).suffix])(ix) ./alphax(ix)); + resultGUI.(['BED', beamInfo(i).suffix]) = reshape(resultGUI.(['BED', beamInfo(i).suffix]), dij.doseGrid.dimensions); + end + + else + % Get Alpha and Beta Values form dij.ax and dij.bx + alphax = reshape(dij.ax, dij.doseGrid.dimensions); + betax = reshape(dij.bx, dij.doseGrid.dimensions); + for i = 1:length(beamInfo) + ix = ~isnan(dij.ax./dij.bx); + if isfield(resultGUI, 'RBExDose') + Dose = resultGUI.(['RBExDose', beamInfo(i).suffix]); + else + Dose = resultGUI.(['physicalDose', beamInfo(i).suffix]); + end + effect = alphax.* Dose + betax.*Dose.^2; + resultGUI.(['BED', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + resultGUI.(['BED', beamInfo(i).suffix])(ix) = full(pln.numOfFractions.*effect(ix)./alphax(ix)); + resultGUI.(['BED', beamInfo(i).suffix]) = reshape(resultGUI.(['BED', beamInfo(i).suffix]), dij.doseGrid.dimensions); + end + end +end % group similar fields together resultGUI = orderfields(resultGUI); From 8fc0b16433f22579048f05537f27d3863fd9ad00 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Fri, 3 Dec 2021 13:47:34 +0100 Subject: [PATCH 03/14] calculation of dij.ax and dij.bx for BED projection in calculation of Dose --- matRad_calcDoseFillDij.m | 8 ++++---- matRad_calcParticleDose.m | 27 +++++++++++++++------------ matRad_calcPhotonDose.m | 13 +++++++++++++ 3 files changed, 32 insertions(+), 16 deletions(-) diff --git a/matRad_calcDoseFillDij.m b/matRad_calcDoseFillDij.m index 79f05cb11..87ce86807 100644 --- a/matRad_calcDoseFillDij.m +++ b/matRad_calcDoseFillDij.m @@ -12,8 +12,8 @@ dij.mLETDose{1}(:,i) = dij.mLETDose{1}(:,i) + stf(i).ray(j).weight(k) * letDoseTmpContainer{1,1}; end - if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... - && strcmp(pln.radiationMode,'carbon') + if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... + ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') % score alpha and beta matrices dij.mAlphaDose{1}(:,i) = dij.mAlphaDose{1}(:,i) + stf(i).ray(j).weight(k) * alphaDoseTmpContainer{1,1}; @@ -33,8 +33,8 @@ dij.mLETDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [letDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}]; end - if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... - && strcmp(pln.radiationMode,'carbon') + if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... + ||isequal(pln.propOpt.bioOptimization,'BED') )&& strcmp(pln.radiationMode,'carbon') dij.mAlphaDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [alphaDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}]; dij.mSqrtBetaDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [betaDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}]; diff --git a/matRad_calcParticleDose.m b/matRad_calcParticleDose.m index 9efdb8aff..63f033ec5 100644 --- a/matRad_calcParticleDose.m +++ b/matRad_calcParticleDose.m @@ -46,7 +46,8 @@ % helper function for energy selection round2 = @(a,b)round(a*10^b)/10^b; -if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... +if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... + ||isequal(pln.propOpt.bioOptimization,'BED')) ... && strcmp(pln.radiationMode,'carbon') alphaDoseTmpContainer = cell(numOfBixelsContainer,dij.numOfScenarios); @@ -56,7 +57,8 @@ dij.mSqrtBetaDose{i} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1); end -elseif isequal(pln.propOpt.bioOptimization,'const_RBExD') && strcmp(pln.radiationMode,'protons') +elseif (isequal(pln.propOpt.bioOptimization,'const_RBExD') ||isequal(pln.propOpt.bioOptimization,'BED') )... + && strcmp(pln.radiationMode,'protons') dij.RBE = 1.1; matRad_cfg.dispInfo('matRad: Using a constant RBE of %g\n',dij.RBE); end @@ -76,12 +78,8 @@ end % generates tissue class matrix for biological optimization -if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... - && strcmp(pln.radiationMode,'carbon') - - if isfield(machine.data,'alphaX') && isfield(machine.data,'betaX') - - matRad_cfg.dispInfo('matRad: loading biological base data... '); + +matRad_cfg.dispInfo('matRad: loading biological base data... '); vTissueIndex = zeros(size(VdoseGrid,1),1); dij.ax = zeros(dij.doseGrid.numOfVoxels,1); dij.bx = zeros(dij.doseGrid.numOfVoxels,1); @@ -93,6 +91,11 @@ dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z); % retrieve photon LQM parameter for the current dose grid voxels [dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid); + +if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')... + ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') + + if isfield(machine.data,'alphaX') && isfield(machine.data,'betaX') for i = 1:size(cst,1) @@ -171,8 +174,8 @@ radDepths = radDepthVdoseGrid{1}(ix); % just use tissue classes of voxels found by ray tracer - if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... - && strcmp(pln.radiationMode,'carbon') + if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')... + ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') vTissueIndex_j = vTissueIndex(ix,:); end @@ -276,8 +279,8 @@ letDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,1} = sparse(VdoseGrid(ix(currIx)),1,bixelLET.*bixelDose,dij.doseGrid.numOfVoxels,1); end - if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... - && strcmp(pln.radiationMode,'carbon') + if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... + ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') % calculate alpha and beta values for bixel k on ray j of [bixelAlpha, bixelBeta] = matRad_calcLQParameter(... currRadDepths,... diff --git a/matRad_calcPhotonDose.m b/matRad_calcPhotonDose.m index f6468d269..e0c33b65c 100644 --- a/matRad_calcPhotonDose.m +++ b/matRad_calcPhotonDose.m @@ -308,6 +308,19 @@ end end +%% calculate Alpha and Beta Values for BED projection + +dij.ax = zeros(dij.doseGrid.numOfVoxels,1); +dij.bx = zeros(dij.doseGrid.numOfVoxels,1); + +cst = matRad_setOverlapPriorities(cst); + +% resizing cst to dose cube resolution +cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,... + dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z); +% retrieve photon LQM parameter for the current dose grid voxels +[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid); + %Close Waitbar if ishandle(figureWait) delete(figureWait); From c483793fb11eb7ea0028b0632aacd1e0a6c1aee9 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Fri, 3 Dec 2021 13:48:02 +0100 Subject: [PATCH 04/14] incorperated BED projection --- matRad_fluenceOptimization.m | 24 ++++++++++-- .../projections/matRad_BEDProjection.m | 25 +++++++++++++ .../projections/matRad_EffectProjection.m | 37 ++++++++++++++----- 3 files changed, 73 insertions(+), 13 deletions(-) create mode 100644 optimization/projections/matRad_BEDProjection.m diff --git a/matRad_fluenceOptimization.m b/matRad_fluenceOptimization.m index aa7d3cc3e..3b7910cf5 100644 --- a/matRad_fluenceOptimization.m +++ b/matRad_fluenceOptimization.m @@ -64,6 +64,7 @@ end obj = obj.setDoseParameters(obj.getDoseParameters()/pln.numOfFractions); + obj.numOfFractions = pln.numOfFractions; cst{i,6}{j} = obj; end @@ -164,7 +165,22 @@ 4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*(dij.physicalDose{1}(V,:)*wOnes))); wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(dij.physicalDose{1}(V,:)*wOnes)))* wOnes; end - +elseif strcmp(pln.propOpt.bioOptimization, 'BED') + if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose') + abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX; + meanBED = mean(pln.numOfFractions.*(dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^2)./cst{ixTarget,5}.alphaX); + BEDTarget = pln.numOfFractions.*doseTarget.*(1 + doseTarget./abr); + elseif isfield(dij, 'RBE') + abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX; + meanBED = mean(pln.numOfFractions.*dij.RBE.*dij.physicalDose{1}(V,:)*wOnes.*(1+dij.RBE.*dij.physicalDose{1}(V,:)*wOnes./abr)); + BEDTarget = pln.numOfFractions.*dij.RBE.*doseTarget.*(1 + dij.RBE.*doseTarget./abr); + else + abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX; + meanBED = mean(pln.numOfFractions.*dij.physicalDose{1}(V,:)*wOnes.*(1+dij.physicalDose{1}(V,:)*wOnes./abr)); + BEDTarget = pln.numOfFractions.*doseTarget.*(1 + doseTarget./abr); + end + bixelWeight = BEDTarget/meanBED; + wInit = wOnes * bixelWeight; else bixelWeight = (doseTarget)/(mean(dij.physicalDose{1}(V,:)*wOnes)); wInit = wOnes * bixelWeight; @@ -186,6 +202,8 @@ backProjection = matRad_ConstantRBEProjection; case 'LEMIV_RBExD' backProjection = matRad_VariableRBEProjection; + case 'BED' + backProjection = matRad_BEDProjection; case 'none' backProjection = matRad_DoseProjection; otherwise @@ -193,7 +211,7 @@ backProjection = matRad_DoseProjection; end - +backProjection.numOfFractions = pln.numOfFractions; %backProjection = matRad_DoseProjection(); optiProb = matRad_OptimizationProblem(backProjection); @@ -221,7 +239,7 @@ wOpt = optimizer.wResult; info = optimizer.resultInfo; -resultGUI = matRad_calcCubes(wOpt,dij); +resultGUI = matRad_calcCubes(wOpt,dij,pln); resultGUI.wUnsequenced = wOpt; resultGUI.usedOptimizer = optimizer; resultGUI.info = info; diff --git a/optimization/projections/matRad_BEDProjection.m b/optimization/projections/matRad_BEDProjection.m new file mode 100644 index 000000000..d6e530edd --- /dev/null +++ b/optimization/projections/matRad_BEDProjection.m @@ -0,0 +1,25 @@ +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==0); + BED(ix) = obj.numOfFractions.*effect(ix)./ dij.ax(ix); + % photon equivalent BED = n * effect / alphax + end + + function wGrad = projectSingleScenarioGradient(obj,dij,doseGrad,scen,w) + doseGrad{scen} = doseGrad{scen}./dij.ax; + wGradEffect = projectSingleScenarioGradient@matRad_EffectProjection(obj,dij,doseGrad,scen,w); + wGrad = obj.numOfFractions.*wGradEffect; + end + end +end + diff --git a/optimization/projections/matRad_EffectProjection.m b/optimization/projections/matRad_EffectProjection.m index 329280329..ed6aa31e0 100644 --- a/optimization/projections/matRad_EffectProjection.m +++ b/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 end From 30fb5ac6c4cbf5be4c716d5d6906e1d9f446fa3e Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Fri, 3 Dec 2021 13:48:28 +0100 Subject: [PATCH 05/14] chnage dose parameters to BED parameters --- .../matRad_constraintFunctions.m | 15 +++++++- .../matRad_constraintJacobian.m | 38 ++++++++++++++++--- .../matRad_getConstraintBounds.m | 9 ++++- .../matRad_objectiveFunction.m | 13 +++++-- .../matRad_objectiveGradient.m | 13 +++++-- 5 files changed, 73 insertions(+), 15 deletions(-) diff --git a/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m b/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m index f31bbd857..4299b0365 100644 --- a/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m +++ b/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m @@ -61,7 +61,8 @@ if (~isequal(obj.name, 'max dose constraint') && ~isequal(obj.name, 'min dose constraint') &&... ~isequal(obj.name, 'max mean dose constraint') && ~isequal(obj.name, 'min mean dose constraint') && ... ~isequal(obj.name, 'min EUD constraint') && ~isequal(obj.name, 'max EUD constraint')) && ... - (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')) + (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')... + && ~isa(optiProb.BP, 'matRad_BEDProjection')) doses = obj.getDoseParameters(); @@ -69,6 +70,18 @@ obj = obj.setDoseParameters(effect); end + + % if we have BED optimization, temporarily replace doses with BED + % Maybe we should put some switch into the classes for that + if (~isequal(obj.name, 'min EUD constraint') && ~isequal(obj.name, 'max EUD constraint'))&& ... + (isa(optiProb.BP,'matRad_BEDProjection')) + + doses = obj.getDoseParameters(); + abr = cst{i,5}.alphaX ./ cst{i,5}.betaX; + BED = optiProb.BP.numOfFractions.*doses.*(1+doses./abr); + obj = obj.setDoseParameters(BED); + end + % if conventional opt: just add constraints of nominal dose %if strcmp(cst{i,6}(j).robustness,'none') diff --git a/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m b/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m index d0ee2dd22..d03b77cb0 100644 --- a/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m +++ b/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m @@ -70,15 +70,22 @@ if (~isequal(obj.name, 'max dose constraint') && ~isequal(obj.name, 'min dose constraint') &&... ~isequal(obj.name, 'max mean dose constraint') && ~isequal(obj.name, 'min mean dose constraint') && ... ~isequal(obj.name, 'min EUD constraint') && ~isequal(obj.name, 'max EUD constraint')) && ... - (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')) + (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')... + && ~isa(optiProb.BP, 'matRad_BEDProjection')) doses = obj.getDoseParameters(); - effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2; - obj = obj.setDoseParameters(effect); end - + % if we have BED optimization, temporarily replace doses with BED + if (~isequal(obj.name, 'min EUD constraint') && ~isequal(obj.name, 'max EUD constraint')) && ... + (isa(optiProb.BP,'matRad_BEDProjection')) + + doses = obj.getDoseParameters(); + abr = cst{i,5}.alphaX ./ cst{i,5}.betaX; + BED = optiProb.BP.numOfFractions.*doses.*(1+doses./abr); + obj = obj.setDoseParameters(BED); + end % if conventional opt: just add constraints of nominal dose %if strcmp(cst{i,6}{j}.robustness,'none') @@ -94,7 +101,8 @@ %Iterate through columns of the sub-jacobian %TODO: Maybe this could all be function of the projection %Objects??? - if isa(optiProb.BP,'matRad_DoseProjection') && ~isempty(jacobSub) || isa(optiProb.BP,'matRad_ConstantRBEProjection') + if isa(optiProb.BP,'matRad_DoseProjection') && ~isempty(jacobSub) || isa(optiProb.BP,'matRad_ConstantRBEProjection')... + || (isa(optiProb.BP, 'matRad_BEDProjection') && ~isfield(dij, 'mAlphaDose') && ~isfield(dij, 'mSqrtBetaDose')) startIx = size(DoseProjection{1},2) + 1; %First append the Projection matrix with sparse zeros @@ -188,7 +196,7 @@ jacob = DoseProjection{1}' * dij.RBE * dij.physicalDose{1}; end -elseif isa(optiProb.BP,'matRad_EffectProjection') +elseif isa(optiProb.BP,'matRad_EffectProjection')&& ~isa(optiProb.BP,'matRad_BEDProjection') if ~isempty(mSqrtBetaDoseProjection{1}) && ~isempty(mAlphaDoseProjection{1}) mSqrtBetaDoseProjection{1} = mSqrtBetaDoseProjection{1}' * dij.mSqrtBetaDose{1} * w; @@ -199,5 +207,23 @@ mSqrtBetaDoseProjection{1}' * dij.mSqrtBetaDose{1}; end +elseif isa(optiProb.BP,'matRad_BEDProjection') + if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose') ... + &&~isempty(mSqrtBetaDoseProjection{1}) && ~isempty(mAlphaDoseProjection{1}) + mSqrtBetaDoseProjection{1} = mSqrtBetaDoseProjection{1}' * dij.mSqrtBetaDose{1} * w; + mSqrtBetaDoseProjection{1} = sparse(voxelID,constraintID,mSqrtBetaDoseProjection{1},... + size(mAlphaDoseProjection{1},1),size(mAlphaDoseProjection{1},2)); + + jacob = obj.numOfFractions.*((mAlphaDoseProjection{1}./dij.ax)' * dij.mAlphaDose{1} +... + (mSqrtBetaDoseProjection{1}./dij.ax)' * dij.mSqrtBetaDose{1}); + elseif isfield(dij, 'RBE') && ~isempty(DoseProjection{1}) + dose = dij.physicalDose{1}*w; + jacob = obj.numOfFractions.*((DoseProjection{1}.*dij.RBE.*(1 + 2*dij.RBE.*(dij.bx./dij.ax).*dose))'*dij.physicalDose{1}); + else + if ~isempty(DoseProjection{1}) + dose = dij.physicalDose{1}*w; + jacob = obj.numOfFractions.*((DoseProjection{1}.*(1+ 2*(dij.bx./dij.ax).*dose))'*dij.physicalDose{1}); + end + end end end diff --git a/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m b/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m index 81b21b7f4..5d54a67cd 100644 --- a/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m +++ b/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 = []; @@ -53,11 +54,15 @@ if isEffectBP 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}))]; diff --git a/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m b/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m index a8ef54158..05deffc67 100644 --- a/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m +++ b/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m @@ -55,15 +55,22 @@ % if we have effect optimization, temporarily replace doses with effect if (~isequal(objective.name, 'Mean Dose') && ~isequal(objective.name, 'EUD')) &&... - (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')) + (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')... + && ~isa(optiProb.BP, 'matRad_BEDProjection')) doses = objective.getDoseParameters(); - effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2; - objective = objective.setDoseParameters(effect); end + % if we have BED optimization, temporarily replace doses with BED + if isa(optiProb.BP, 'matRad_BEDProjection') + doses = objective.getDoseParameters(); + abr = cst{i,5}.alphaX ./ cst{i,5}.betaX; + BED = optiProb.BP.numOfFractions.*doses.*(1+doses./abr); + objective = objective.setDoseParameters(BED); + end + % if conventional opt: just sum objectiveectives of nominal dose %if strcmp(cst{i,6}{j}.robustness,'none') diff --git a/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m b/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m index da4d95a1b..ef7cec931 100644 --- a/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m +++ b/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m @@ -58,15 +58,22 @@ if isa(objective,'DoseObjectives.matRad_DoseObjective') % if we have effect optimization, temporarily replace doses with effect if (~isequal(objective.name, 'Mean Dose') && ~isequal(objective.name, 'EUD')) &&... - (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')) + (isa(optiProb.BP,'matRad_EffectProjection') && ~isa(optiProb.BP,'matRad_VariableRBEProjection')... + && ~isa(optiProb.BP, 'matRad_BEDProjection')) doses = objective.getDoseParameters(); - effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2; - objective = objective.setDoseParameters(effect); end + % if we have BED optimization, temporarily replace doses with BED + if isa(optiProb.BP, 'matRad_BEDProjection') + doses = objective.getDoseParameters(); + abr = cst{i,5}.alphaX ./ cst{i,5}.betaX; + BED = optiProb.BP.numOfFractions.*doses.*(1+doses./abr); + objective = objective.setDoseParameters(BED); + end + %dose in VOI d_i = d{1}(cst{i,4}{1}); From 64fb4709ee4737912585235513eef0bee10fb624 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Fri, 3 Dec 2021 13:48:50 +0100 Subject: [PATCH 06/14] showinf DVH of BED --- matRadGUI.m | 13 ++++++++++++- matRad_indicatorWrapper.m | 2 ++ matRad_showDVH.m | 2 ++ 3 files changed, 16 insertions(+), 1 deletion(-) diff --git a/matRadGUI.m b/matRadGUI.m index 9061d4516..0e7147a7f 100644 --- a/matRadGUI.m +++ b/matRadGUI.m @@ -182,6 +182,10 @@ end set(handles.popUpMachine,'String',handles.Machines); +%Add Bed as Biogical Optimization method +stringBioOpt = get(handles.popMenuBioOpt, 'String'); +stringBioOpt{5} = 'BED'; +set(handles.popMenuBioOpt, 'String', stringBioOpt); vChar = get(handles.editGantryAngle,'String'); if strcmp(vChar(1,1),'0') && length(vChar)==6 @@ -2178,6 +2182,8 @@ function getPlnFromGUI(handles) if (~strcmp(pln.radiationMode,'photons')) contentBioOpt = get(handles.popMenuBioOpt,'String'); pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:}; +elseif strcmp(pln.propOpt.bioOptimization, 'BED') + pln.propOpt.bioOptimization = 'BED'; else pln.propOpt.bioOptimization = 'none'; end @@ -2322,7 +2328,12 @@ function btnDVH_Callback(~, ~, handles) %check if one of the default fields is selected if sum(strcmp(SelectedCube,{'physicalDose','effect','RBE,','RBExDose','alpha','beta'})) > 0 resultGUI_SelectedCube.physicalDose = resultGUI.physicalDose; - resultGUI_SelectedCube.RBExDose = resultGUI.RBExDose; + if isfield(resultGUI,'RBExDose') + resultGUI_SelectedCube.RBExDose = resultGUI.RBExDose; + end + elseif strcmp(SelectedCube,'BED') + resultGUI_SelectedCube.physicalDose = resultGUI.physicalDose; + resultGUI_SelectedCube.BED = resultGUI.BED; else Idx = find(SelectedCube == '_'); SelectedSuffix = SelectedCube(Idx(1):end); diff --git a/matRad_indicatorWrapper.m b/matRad_indicatorWrapper.m index 0bd00f64e..be72fec73 100644 --- a/matRad_indicatorWrapper.m +++ b/matRad_indicatorWrapper.m @@ -38,6 +38,8 @@ if isfield(resultGUI,'RBExDose') doseCube = resultGUI.RBExDose; +elseif isfield(resultGUI,'BED') + doseCube = resultGUI.BED; else doseCube = resultGUI.physicalDose; end diff --git a/matRad_showDVH.m b/matRad_showDVH.m index 312949d16..c2f7115ff 100644 --- a/matRad_showDVH.m +++ b/matRad_showDVH.m @@ -96,6 +96,8 @@ function matRad_showDVH(dvh,cst,pln,lineStyleIndicator) if exist('pln','var') && ~isempty(pln) if strcmp(pln.propOpt.bioOptimization,'none') xlabel('Dose [Gy]','FontSize',fontSizeValue); + elseif strcmp(pln.propOpt.bioOptimization,'BED') + xlabel('Dose Gy]','FontSize',fontSizeValue); else xlabel('RBE x Dose [Gy(RBE)]','FontSize',fontSizeValue); end From fc8358d11ed4e463fe1ffcd7f00346b1a6a73bd4 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Thu, 17 Mar 2022 14:00:24 +0100 Subject: [PATCH 07/14] Update matRad_calcParticleDoseMC.m --- matRad_calcParticleDoseMC.m | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/matRad_calcParticleDoseMC.m b/matRad_calcParticleDoseMC.m index bd705be3f..10f2fc845 100644 --- a/matRad_calcParticleDoseMC.m +++ b/matRad_calcParticleDoseMC.m @@ -288,6 +288,17 @@ dij.physicalDose{1} = dij.physicalDose{1}(:,MCsquareOrder); end +% calculate Alpha and Beta Values for BED projection +dij.ax = zeros(dij.doseGrid.numOfVoxels,1); +dij.bx = zeros(dij.doseGrid.numOfVoxels,1); +cst = matRad_setOverlapPriorities(cst); +% resizing cst to dose cube resolution +cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,... + dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z); +% retrieve photon LQM parameter for the current dose grid voxels +[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid); + + matRad_cfg.dispInfo('matRad: done!\n'); try From f91cb0d20ef9f0535872f87e9e91e70c14241134 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Thu, 17 Mar 2022 14:00:52 +0100 Subject: [PATCH 08/14] Update matRad_calcPhotonDoseMC.m --- matRad_calcPhotonDoseMC.m | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/matRad_calcPhotonDoseMC.m b/matRad_calcPhotonDoseMC.m index 1fbc0595c..15acea838 100644 --- a/matRad_calcPhotonDoseMC.m +++ b/matRad_calcPhotonDoseMC.m @@ -328,6 +328,17 @@ end end +% calculate Alpha and Beta Values for BED projection +dij.ax = zeros(dij.doseGrid.numOfVoxels,1); +dij.bx = zeros(dij.doseGrid.numOfVoxels,1); +cst = matRad_setOverlapPriorities(cst); +% resizing cst to dose cube resolution +cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,... + dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z); +% retrieve photon LQM parameter for the current dose grid voxels +[dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid); + + matRad_cfg.dispInfo('matRad: MC photon dose calculation done!\n'); matRad_cfg.dispInfo(evalc('toc')); From 4f208395332baaa0c391139b986a19fccdb09ac0 Mon Sep 17 00:00:00 2001 From: JenHardt <84377305+JenHardt@users.noreply.github.com> Date: Mon, 2 May 2022 10:02:22 +0200 Subject: [PATCH 09/14] Update matRadGUI.m --- matRadGUI.m | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/matRadGUI.m b/matRadGUI.m index df2cc7c98..e0bc1b530 100644 --- a/matRadGUI.m +++ b/matRadGUI.m @@ -720,10 +720,10 @@ function popupRadMode_Callback(hObject, eventdata, handles) switch RadIdentifier case 'photons' - set(handles.popMenuBioOpt,'Enable','off'); + set(handles.popMenuBioOpt,'Enable','on'); ix = find(strcmp(contentPopUp,'none')); set(handles.popMenuBioOpt,'Value',ix); - set(handles.btnSetTissue,'Enable','off'); + set(handles.btnSetTissue,'Enable','on'); set(handles.btnRunSequencing,'Enable','on'); set(handles.btnRunDAO,'Enable','on'); @@ -1876,10 +1876,10 @@ function UpdateState(handles) set(handles.btnSetTissue,'Enable','on'); elseif strcmp(pln.radiationMode,'protons') set(handles.popMenuBioOpt,'Enable','on'); - set(handles.btnSetTissue,'Enable','off'); + set(handles.btnSetTissue,'Enable','on'); else - set(handles.popMenuBioOpt,'Enable','off'); - set(handles.btnSetTissue,'Enable','off'); + set(handles.popMenuBioOpt,'Enable','on'); + set(handles.btnSetTissue,'Enable','on'); end cMapControls = allchild(handles.uipanel_colormapOptions); @@ -2031,6 +2031,7 @@ function setPln(handles) set(handles.popupRadMode,'Value',find(strcmp(get(handles.popupRadMode,'String'),pln.radiationMode))); set(handles.popUpMachine,'Value',find(strcmp(get(handles.popUpMachine,'String'),pln.machine))); + if ~strcmp(pln.propOpt.bioOptimization,'none') set(handles.popMenuBioOpt,'Enable','on'); contentPopUp = get(handles.popMenuBioOpt,'String'); @@ -2038,9 +2039,10 @@ function setPln(handles) set(handles.popMenuBioOpt,'Value',ix); set(handles.btnSetTissue,'Enable','on'); else - set(handles.popMenuBioOpt,'Enable','off'); - set(handles.btnSetTissue,'Enable','off'); + set(handles.popMenuBioOpt,'Enable','on'); + set(handles.btnSetTissue,'Enable','on'); end + %% enable sequencing and DAO button if radiation mode is set to photons if strcmp(pln.radiationMode,'photons') && pln.propOpt.runSequencing set(handles.btnRunSequencing,'Enable','on'); @@ -2179,14 +2181,19 @@ function getPlnFromGUI(handles) contents = get(handles.popUpMachine,'String'); pln.machine = contents{get(handles.popUpMachine,'Value')}; +%{ if (~strcmp(pln.radiationMode,'photons')) contentBioOpt = get(handles.popMenuBioOpt,'String'); pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:}; -elseif strcmp(pln.propOpt.bioOptimization, 'BED') +elseif isfield(pln, 'propOpt') && strcmp(pln.propOpt.bioOptimization, 'BED') pln.propOpt.bioOptimization = 'BED'; else pln.propOpt.bioOptimization = 'none'; end +%} + +contentBioOpt = get(handles.popMenuBioOpt,'String'); +pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:}; pln.propOpt.runSequencing = logical(get(handles.btnRunSequencing,'Value')); pln.propOpt.runDAO = logical(get(handles.btnRunDAO,'Value')); @@ -2887,6 +2894,10 @@ function SaveResultToGUI(~, ~, ~) if isfield(resultGUI,'RBExDose') resultGUI.(['RBExDose' Suffix]) = resultGUI.RBExDose; end + if isfield(resultGUI,'BED') + resultGUI.(['BED' Suffix]) = resultGUI.BED; + end + if strcmp(pln.radiationMode,'carbon') == 1 if isfield(resultGUI,'effect') From b8b365db752c2c0088a17a161ebac62ef3de98f9 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 30 Jun 2023 11:47:35 +0200 Subject: [PATCH 10/14] warning message BED photon equivalent calculated --- matRad_calcCubes.m | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/matRad_calcCubes.m b/matRad_calcCubes.m index 2d20cd61d..be6e78dac 100644 --- a/matRad_calcCubes.m +++ b/matRad_calcCubes.m @@ -30,6 +30,8 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +matRad_cfg = MatRad_Config.instance(); + if ~exist('sceneNum', 'var') scenNum = 1; end @@ -112,7 +114,7 @@ resultGUI.(['BED', beamInfo(i).suffix])(ix) = full(pln.numOfFractions.*resultGUI.(['effect', beamInfo(i).suffix])(ix) ./alphax(ix)); resultGUI.(['BED', beamInfo(i).suffix]) = reshape(resultGUI.(['BED', beamInfo(i).suffix]), dij.doseGrid.dimensions); end - + matRad_cfg.dispWarning('Photon Equiavlent BED calculated'); else % Get Alpha and Beta Values form dij.ax and dij.bx alphax = reshape(dij.ax, dij.doseGrid.dimensions); @@ -129,6 +131,9 @@ resultGUI.(['BED', beamInfo(i).suffix])(ix) = full(pln.numOfFractions.*effect(ix)./alphax(ix)); resultGUI.(['BED', beamInfo(i).suffix]) = reshape(resultGUI.(['BED', beamInfo(i).suffix]), dij.doseGrid.dimensions); end + if isfield(resultGUI, 'RBExDose') + matRad_cfg.dispWarning('Photon Equiavlent BED calculated'); + end end end % group similar fields together From 9a10d77003eec1c20e64fdab0ad96208078cdcc2 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 30 Jun 2023 11:48:02 +0200 Subject: [PATCH 11/14] renaming LEMIV_BED --- matRad_calcDoseFillDij.m | 4 ++-- matRad_calcParticleDose.m | 10 +++++----- matRad_showDVH.m | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/matRad_calcDoseFillDij.m b/matRad_calcDoseFillDij.m index 87ce86807..d4749ef2a 100644 --- a/matRad_calcDoseFillDij.m +++ b/matRad_calcDoseFillDij.m @@ -13,7 +13,7 @@ end if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... - ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') + ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED') ) && strcmp(pln.radiationMode,'carbon') % score alpha and beta matrices dij.mAlphaDose{1}(:,i) = dij.mAlphaDose{1}(:,i) + stf(i).ray(j).weight(k) * alphaDoseTmpContainer{1,1}; @@ -34,7 +34,7 @@ end if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... - ||isequal(pln.propOpt.bioOptimization,'BED') )&& strcmp(pln.radiationMode,'carbon') + ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED') )&& strcmp(pln.radiationMode,'carbon') dij.mAlphaDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [alphaDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}]; dij.mSqrtBetaDose{1}(:,(ceil(counter/numOfBixelsContainer)-1)*numOfBixelsContainer+1:counter) = [betaDoseTmpContainer{1:mod(counter-1,numOfBixelsContainer)+1,1}]; diff --git a/matRad_calcParticleDose.m b/matRad_calcParticleDose.m index 81ce531a1..ad7f420a0 100644 --- a/matRad_calcParticleDose.m +++ b/matRad_calcParticleDose.m @@ -47,7 +47,7 @@ round2 = @(a,b)round(a*10^b)/10^b; if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... - ||isequal(pln.propOpt.bioOptimization,'BED')) ... + ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED')) ... && strcmp(pln.radiationMode,'carbon') alphaDoseTmpContainer = cell(numOfBixelsContainer,dij.numOfScenarios); @@ -57,7 +57,7 @@ dij.mSqrtBetaDose{i} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1); end -elseif (isequal(pln.propOpt.bioOptimization,'const_RBExD') ||isequal(pln.propOpt.bioOptimization,'BED') )... +elseif (isequal(pln.propOpt.bioOptimization,'const_RBExD') ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED') )... && strcmp(pln.radiationMode,'protons') dij.RBE = 1.1; matRad_cfg.dispInfo('matRad: Using a constant RBE of %g\n',dij.RBE); @@ -93,7 +93,7 @@ [dij.ax,dij.bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1,VdoseGrid); if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')... - ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') + ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED') ) && strcmp(pln.radiationMode,'carbon') if isfield(machine.data,'alphaX') && isfield(machine.data,'betaX') @@ -204,7 +204,7 @@ % just use tissue classes of voxels found by ray tracer if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')... - ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') + ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED') ) && strcmp(pln.radiationMode,'carbon') vTissueIndex_j = vTissueIndex(ix,:); end @@ -366,7 +366,7 @@ if (isequal(pln.propOpt.bioOptimization,'LEMIV_effect') || isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD') ... - ||isequal(pln.propOpt.bioOptimization,'BED') ) && strcmp(pln.radiationMode,'carbon') + ||isequal(pln.propOpt.bioOptimization,'LEMIV_BED') ) && strcmp(pln.radiationMode,'carbon') % calculate alpha and beta values for bixel k on ray j of [bixelAlpha, bixelBeta] = matRad_calcLQParameter(... currRadDepths,... diff --git a/matRad_showDVH.m b/matRad_showDVH.m index c2f7115ff..1a463e60d 100644 --- a/matRad_showDVH.m +++ b/matRad_showDVH.m @@ -96,7 +96,7 @@ function matRad_showDVH(dvh,cst,pln,lineStyleIndicator) if exist('pln','var') && ~isempty(pln) if strcmp(pln.propOpt.bioOptimization,'none') xlabel('Dose [Gy]','FontSize',fontSizeValue); - elseif strcmp(pln.propOpt.bioOptimization,'BED') + elseif strcmp(pln.propOpt.bioOptimization,'LEMIV_BED') xlabel('Dose Gy]','FontSize',fontSizeValue); else xlabel('RBE x Dose [Gy(RBE)]','FontSize',fontSizeValue); From 5e1e1654a31452c086bb98255726f2b281794261 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 30 Jun 2023 11:48:21 +0200 Subject: [PATCH 12/14] no NANs in dose Grad --- optimization/projections/matRad_BEDProjection.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/optimization/projections/matRad_BEDProjection.m b/optimization/projections/matRad_BEDProjection.m index d6e530edd..b8236d3e8 100644 --- a/optimization/projections/matRad_BEDProjection.m +++ b/optimization/projections/matRad_BEDProjection.m @@ -16,8 +16,10 @@ end function wGrad = projectSingleScenarioGradient(obj,dij,doseGrad,scen,w) - doseGrad{scen} = doseGrad{scen}./dij.ax; - wGradEffect = projectSingleScenarioGradient@matRad_EffectProjection(obj,dij,doseGrad,scen,w); + ix = ~(dij.ax==0); + doseGradtmp{scen} = zeros(size(doseGrad{scen})); + doseGradtmp{scen}(ix) = doseGrad{scen}(ix)./dij.ax(ix); + wGradEffect = projectSingleScenarioGradient@matRad_EffectProjection(obj,dij,doseGradtmp,scen,w); wGrad = obj.numOfFractions.*wGradEffect; end end From 6ea309ed459849407ddef24f3326d431ad445f68 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 30 Jun 2023 11:48:38 +0200 Subject: [PATCH 13/14] renaming LEMIV_BED --- matRadGUI.m | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/matRadGUI.m b/matRadGUI.m index e0bc1b530..b8f1d3816 100644 --- a/matRadGUI.m +++ b/matRadGUI.m @@ -184,7 +184,7 @@ %Add Bed as Biogical Optimization method stringBioOpt = get(handles.popMenuBioOpt, 'String'); -stringBioOpt{5} = 'BED'; +stringBioOpt{5} = 'LEMIV_BED'; set(handles.popMenuBioOpt, 'String', stringBioOpt); vChar = get(handles.editGantryAngle,'String'); @@ -2181,16 +2181,6 @@ function getPlnFromGUI(handles) contents = get(handles.popUpMachine,'String'); pln.machine = contents{get(handles.popUpMachine,'Value')}; -%{ -if (~strcmp(pln.radiationMode,'photons')) - contentBioOpt = get(handles.popMenuBioOpt,'String'); - pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:}; -elseif isfield(pln, 'propOpt') && strcmp(pln.propOpt.bioOptimization, 'BED') - pln.propOpt.bioOptimization = 'BED'; -else - pln.propOpt.bioOptimization = 'none'; -end -%} contentBioOpt = get(handles.popMenuBioOpt,'String'); pln.propOpt.bioOptimization = contentBioOpt{get(handles.popMenuBioOpt,'Value'),:}; From 14a50ab68983309fe1d87bce9d847b60cfe62f73 Mon Sep 17 00:00:00 2001 From: JenHardt Date: Fri, 30 Jun 2023 11:50:03 +0200 Subject: [PATCH 14/14] rnaming LEMIV and also check ax, bx parameters for BED, not only effect --- matRad_fluenceOptimization.m | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/matRad_fluenceOptimization.m b/matRad_fluenceOptimization.m index 757e1568b..fa571160e 100644 --- a/matRad_fluenceOptimization.m +++ b/matRad_fluenceOptimization.m @@ -169,7 +169,15 @@ 4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*(dij.physicalDose{1}(V,:)*wOnes))); wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(dij.physicalDose{1}(V,:)*wOnes)))* wOnes; end -elseif strcmp(pln.propOpt.bioOptimization, 'BED') +elseif strcmp(pln.propOpt.bioOptimization, 'LEMIV_BED') + + % retrieve photon LQM parameter + [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1); + if ~isequal(dij.ax(dij.ax~=0),ax(dij.ax~=0)) || ... + ~isequal(dij.bx(dij.bx~=0),bx(dij.bx~=0)) + matRad_cfg.dispError('Inconsistent biological parameter - please recalculate dose influence matrix!\n'); + end + if isfield(dij, 'mAlphaDose') && isfield(dij, 'mSqrtBetaDose') abr = cst{ixTarget,5}.alphaX./cst{ixTarget,5}.betaX; meanBED = mean(pln.numOfFractions.*(dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^2)./cst{ixTarget,5}.alphaX); @@ -206,7 +214,7 @@ backProjection = matRad_ConstantRBEProjection; case 'LEMIV_RBExD' backProjection = matRad_VariableRBEProjection; - case 'BED' + case 'LEMIV_BED' backProjection = matRad_BEDProjection; case 'none' backProjection = matRad_DoseProjection;