diff --git a/examples/matRad_example16_paretoOptimization.m b/examples/matRad_example16_paretoOptimization.m new file mode 100644 index 000000000..5e735816b --- /dev/null +++ b/examples/matRad_example16_paretoOptimization.m @@ -0,0 +1,121 @@ +%% Example: Photon Treatment Plan +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2017 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% In this example we will show +% (i) how to load patient data into matRad +% (ii) how to setup a photon dose calculation and +% (iii) how to inversely optimize beamlet intensities +% (iv) how to visually and quantitatively evaluate the result +%global timed; +%timed = []; +%% Patient Data Import +% Let's begin with a clear Matlab environment. Then, import the TG119 +% phantom into your workspace. The phantom is comprised of a 'ct' and 'cst' +% structure defining the CT images and the structure set. Make sure the +% matRad root directory with all its subdirectories is added to the Matlab +% search path. + +matRad_rc; %If this throws an error, run it from the parent directory first to set the paths +matRad_rc; +matRad_cfg = MatRad_Config.instance(); +%matRad_cfg.propOpt.defaultMaxIter = 50000; +%% +load('TG119.mat'); +%% +cst{1,6}{1} = struct(DoseObjectives.matRad_EUD(1000,0)); +cst{3,6}{1} = struct(DoseObjectives.matRad_MeanDose(1000,0)); +%% +%cst{2,6}{2} = struct(DoseConstraints.matRad_MinMaxDose(45,55)); + +%pln.radiationMode = 'protons'; +pln.radiationMode = 'photons'; +pln.machine = 'Generic'; + +quantityOpt = 'physicalDose'; +modelName = 'none'; + +%pln.propDoseCalc.calcLET = 0; + +pln.numOfFractions = 30; +pln.propStf.gantryAngles = [0:40:359]; +pln.propStf.couchAngles = zeros(1,numel(pln.propStf.gantryAngles)); +pln.propStf.bixelWidth = 5; + + +pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles); +pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0); + +pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm] +pln.propDoseCalc.doseGrid.resolution.y = 5; % [mm] +pln.propDoseCalc.doseGrid.resolution.z = 5; % [mm] + +%% +% Enable sequencing and disable direct aperture optimization (DAO) for now. +% A DAO optimization is shown in a seperate example. +pln.propOpt.runSequencing = 1; +pln.propOpt.runDAO = 0; + +% retrieve bio model parameters +pln.bioParam = matRad_bioModel(pln.radiationMode,quantityOpt, modelName); + +% retrieve scenarios for dose calculation and optimziation +pln.multScen = matRad_multScen(ct,'nomScen'); + +%% +% and et voila our treatment plan structure is ready. Lets have a look: +display(pln); + + +%% Generate Beam Geometry STF +% The steering file struct comprises the complete beam geometry along with +% ray position, pencil beam positions and energies, source to axis distance (SAD) etc. +stf = matRad_generateStf(ct,cst,pln); + +%% +% Let's display the beam geometry information of the 6th beam +%display(stf(6)); + +%% Dose Calculation +% Let's generate dosimetric information by pre-computing dose influence +% matrices for unit beamlet intensities. Having dose influences available +% allows subsequent inverse optimization. +%dij = matRad_calcParticleDose(ct,stf,pln,cst); +dij = matRad_calcPhotonDose(ct,stf,pln,cst); +%% Inverse Optimization for IMRT +% The goal of the fluence optimization is to find a set of beamlet/pencil +% beam weights which yield the best possible dose distribution according to +% the clinical objectives and constraints underlying the radiation +% treatment. Once the optimization has finished, trigger once the GUI to +% visualize the optimized dose cubes. +%% Paretooptimization +% The goal of this step is to define a grid of penalty values that +% are then evaluated using the matRad_paretoGeneration method +% The VOI and their respective penalties are defined in the following way +% It is possible to have more than one objective function per VOI +% penVal stores the Grid which is then passed on. penGrid contains an +% version easier to visualize, however harder to loop over +%% Add constraints +cst{1,6}{2} = struct(DoseConstraints.matRad_MinMaxDose(0,40)); +cst{3,6}{2} = struct(DoseConstraints.matRad_MinMaxDose(0,45)); +cst{2,6}{2} = struct(DoseConstraints.matRad_MinMaxDose(45,57)); + +%% +[resultGUI,retStruct] = matRad_paretoOptimization(dij,cst,pln,3); + +%% +matRadGUI +%% +matRadParetoGUI +% matRad_UIInterpolation(retStruct,dij,pln,ct,matRad_setOverlapPriorities(cst),retStruct.optiProb) diff --git a/examples/matRad_example17_lexicographicOptimization.m b/examples/matRad_example17_lexicographicOptimization.m new file mode 100644 index 000000000..30aa92950 --- /dev/null +++ b/examples/matRad_example17_lexicographicOptimization.m @@ -0,0 +1,51 @@ +matRad_rc; %If this throws an error, run it from the parent directory first to set the paths +load('TG119.mat'); + + +pln.radiationMode = 'photons'; +pln.machine = 'Generic'; + +quantityOpt = 'physicalDose'; +modelName = 'none'; + +% +pln.numOfFractions = 30; +pln.propStf.gantryAngles = [0:40:359]; +pln.propStf.couchAngles = zeros(1,numel(pln.propStf.gantryAngles)); +pln.propStf.bixelWidth = 5; + +pln.propStf.numOfBeams = numel(pln.propStf.gantryAngles); +pln.propStf.isoCenter = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0); + +pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm] +pln.propDoseCalc.doseGrid.resolution.y = 5; % [mm] +pln.propDoseCalc.doseGrid.resolution.z = 5; % [mm] + +pln.propSeq.runSequencing = 1; +pln.propOpt.runDAO = 0; + +% retrieve bio model parameters +pln.bioParam = matRad_bioModel(pln.radiationMode,quantityOpt, modelName); + +% retrieve scenarios for dose calculation and optimziation +pln.multScen = matRad_multScen(ct,'nomScen'); +%% +stf = matRad_generateStf(ct,cst,pln); +dij = matRad_calcPhotonDose(ct,stf,pln,cst); + +%% +PriorityList1 = matRad_PriorityList1(); +PriorityList1.addObjective(1,DoseObjectives.matRad_SquaredDeviation(100,50),4,2); +PriorityList1.addObjective(2,DoseObjectives.matRad_EUD(100,0),10,1); +PriorityList1.addObjective(3,DoseObjectives.matRad_MeanDose(100,0),5,3); + +%constraints +%PriorityList1.addConstraint(DoseConstraints.matRad_MinMaxDose(0,40),1); +%PriorityList1.addConstraint(DoseConstraints.matRad_MinMaxDose(45,57),2); +%PriorityList1.addConstraint(DoseConstraints.matRad_MinMaxDose(0,45),3); +%% +PriorityList1.printPriorityList(cst) +%% +[resultGUIs1,resultGUIs2,cst1,cst2,PriorityList2] = matRad_2pecOptimization(PriorityList1,dij,cst,pln); +%% +PriorityList1.printPriorityList(cst) \ No newline at end of file diff --git a/matRad/MatRad_Config.m b/matRad/MatRad_Config.m index 183734f63..184585845 100644 --- a/matRad/MatRad_Config.m +++ b/matRad/MatRad_Config.m @@ -210,18 +210,16 @@ function setDefaultProperties(obj) %Optimization Options obj.defaults.propOpt.optimizer = 'IPOPT'; - obj.defaults.propOpt.maxIter = 500; + obj.defaults.propOpt.maxIter = 1000; obj.defaults.propOpt.runDAO = 0; obj.defaults.propOpt.clearUnusedVoxels = false; %Sequencing Options obj.defaults.propSeq.sequencer = 'siochi'; + + obj.defaults.samplingScenarios = 25; - - obj.disableGUI = false; - - obj.defaults.samplingScenarios = 25; obj.devMode = false; obj.eduMode = false; diff --git a/matRad/gui/matRadParetoGUI.m b/matRad/gui/matRadParetoGUI.m new file mode 100644 index 000000000..92a1f0c37 --- /dev/null +++ b/matRad/gui/matRadParetoGUI.m @@ -0,0 +1,47 @@ +function hGUI = matRadParetoGUI(varargin) +% matRad compatability function to call the matRad_MainGUI +% The function checks input parameters and handles the GUI as a +% singleton, so following calls will not create new windows +% +% call +% hGUI = matRadGUI +% matRadGUI +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +persistent hMatRadParetoGUI; + +%Initialize navigation structures -> Check for retStruct etc. +retStruct = evalin('base','retStruct'); +finds = retStruct.finds; +weights = retStruct.weights; + +%Matrix storing objective function values of all calculated plans +assignin('base','finds',finds); + +%Matrix storing the plans' associated weights +assignin('base','weights',weights); + +matRad_ParetoGUI(); + + + + + + + diff --git a/matRad/gui/pareto/matRad_ParetoDVHWidget.m b/matRad/gui/pareto/matRad_ParetoDVHWidget.m new file mode 100644 index 000000000..1d0acc894 --- /dev/null +++ b/matRad/gui/pareto/matRad_ParetoDVHWidget.m @@ -0,0 +1,75 @@ +classdef matRad_ParetoDVHWidget < matRad_Widget + % matRad_ParetoDVHWidget class creates a widget for the Pareto GUI to + % display plan DVH's + % + % + % References + % - + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + DVHPlotAxes + end + + methods + function this = matRad_ParetoDVHWidget(handleParent) + this = this@matRad_Widget(handleParent); + this.DVHPlotAxes = axes(this.widgetHandle,'Position',[0.09 0.1 0.9 0.89]); + this.plotInitialDVH(); + end + + function this = initialize(this) + this.update(); + end + + end + + methods (Access = protected) + function this = createLayout(this) + + parent = this.widgetHandle; + + matRad_cfg = MatRad_Config.instance(); + + + this.createHandles(); + end + + function this = doUpdate(this,evt) + %getFromWorkspace(this); + %updateInWorkspace(this); + end + end + + methods (Access = private) + function plotInitialDVH(this) + ParetoHelperObject = evalin('base','ParetoHelperObject'); + %Shows the DVH for the current plan + resultGUI = matRad_calcCubes(ParetoHelperObject.currentWeights,evalin('base','dij')); + dvh = matRad_calcDVH(evalin('base','cst'),resultGUI.physicalDose,'cum'); + + matRad_showDVH(this.DVHPlotAxes,... + dvh,evalin('base','cst'),evalin('base','pln'),ParetoHelperObject.linestyle); + if ParetoHelperObject.linestyle < 4 + ParetoHelperObject.linestyle = ParetoHelperObject.linestyle + 1; + else + ParetoHelperObject.linestyle = 1; + end + end + end + +end + + diff --git a/matRad/gui/pareto/matRad_ParetoData.m b/matRad/gui/pareto/matRad_ParetoData.m new file mode 100644 index 000000000..2a1d2fbf5 --- /dev/null +++ b/matRad/gui/pareto/matRad_ParetoData.m @@ -0,0 +1,80 @@ +classdef matRad_ParetoData < handle + % matRad_ParetoData implements a class that allows easy storing of + % variables related to the pareto Navigation + % + % References + % - + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + currentWeights + allWeights + currentPoint %objective function values of last calculated plan + allPoints % matrix storing all objective function values + availablePoints % stores only the points not blocked by boundaries + initialUpbound + availableUpbound + upboundSlider + lowboundSliderInit + lowboundSlider + linestyle + end + + methods + function this = matRad_ParetoData(allPoints,allWeights) + %initialize first set + this.allWeights = allWeights; + this.allPoints = allPoints; + this.availablePoints = allPoints; + this.initialUpbound = max(allPoints,[],1); + this.availableUpbound = max(allPoints,[],1); + this.upboundSlider = max(allPoints,[],1); + this.lowboundSliderInit = min(allPoints,[],1); + this.lowboundSlider = min(allPoints,[],1); + this.linestyle = 1; + + %calculate an initial point based on the shortes distance to the + %polyhedral center + center = mean(this.allPoints,1); + shiftedPoints = this.allPoints-center; + distToCenter = sum(shiftedPoints.^2,2); + [~,idx] = min(distToCenter); + this.currentPoint = this.allPoints(idx,:); + this.currentWeights = this.allWeights(:,idx); + + end + + + function [sliderLowBound,sliderUpBound] = restrictObjective(this,i,bound) + %restrict the upper value of a specific slider + this.availableUpbound(i) = bound; + this.availablePoints = this.availablePoints(this.availablePoints(:,i) <= bound,:); + this.upboundSlider= max([this.availablePoints;this.currentPoint'],[],1); + this.upboundSlider(i) = bound; + this.lowboundSlider = min([this.availablePoints;this.currentPoint'],[],1); + sliderLowBound = this.lowboundSlider; + sliderUpBound = this.upboundSlider; + end + + function releaseObjectiveBounds(this) + %reset objective bounds + this.availableUpbound = this.initialUpbound; + this.availablePoints = this.allPoints; + this.upboundSlider = this.initialUpbound; + this.lowboundSlider = this.lowboundSliderInit; + end + + end +end \ No newline at end of file diff --git a/matRad/gui/pareto/matRad_ParetoGUI.m b/matRad/gui/pareto/matRad_ParetoGUI.m new file mode 100644 index 000000000..c46e41237 --- /dev/null +++ b/matRad/gui/pareto/matRad_ParetoGUI.m @@ -0,0 +1,136 @@ +classdef matRad_ParetoGUI < handle +% matRad_ParetoGUI Pareto navigation Graphical User Interface configuration class +% This Class is used to create the Pareto GUI interface where all the widgets +% are called and created. +% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2019 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + properties + guiHandle + lastStoragePath = [] + end + + properties (Access = private) + LogoWidget + ParetoSliceWidget + ParetoSliderWidget + ParetoDVHWidget + ParetoNavigationButtonWidget + end + + methods + function this = matRad_ParetoGUI(varargin) + %Panel for Main Widget + + matRad_cfg = MatRad_Config.instance(); + + this.guiHandle = figure(... + 'Units','normalized',... + 'Position',[0.005 0.04 0.99 0.9],... %approximate fullscreen position + 'Visible','on',... + 'Color',matRad_cfg.gui.backgroundColor,... + 'IntegerHandle','off',... + 'Colormap',[0 0 0.5625;0 0 0.625;0 0 0.6875;0 0 0.75;0 0 0.8125;0 0 0.875;0 0 0.9375;0 0 1;0 0.0625 1;0 0.125 1;0 0.1875 1;0 0.25 1;0 0.3125 1;0 0.375 1;0 0.4375 1;0 0.5 1;0 0.5625 1;0 0.625 1;0 0.6875 1;0 0.75 1;0 0.8125 1;0 0.875 1;0 0.9375 1;0 1 1;0.0625 1 1;0.125 1 0.9375;0.1875 1 0.875;0.25 1 0.8125;0.3125 1 0.75;0.375 1 0.6875;0.4375 1 0.625;0.5 1 0.5625;0.5625 1 0.5;0.625 1 0.4375;0.6875 1 0.375;0.75 1 0.3125;0.8125 1 0.25;0.875 1 0.1875;0.9375 1 0.125;1 1 0.0625;1 1 0;1 0.9375 0;1 0.875 0;1 0.8125 0;1 0.75 0;1 0.6875 0;1 0.625 0;1 0.5625 0;1 0.5 0;1 0.4375 0;1 0.375 0;1 0.3125 0;1 0.25 0;1 0.1875 0;1 0.125 0;1 0.0625 0;1 0 0;0.9375 0 0;0.875 0 0;0.8125 0 0;0.75 0 0;0.6875 0 0;0.625 0 0;0.5625 0 0],... + 'MenuBar','none',... + 'Name','matRad Pareto Navigator (NOT FOR CLINICAL USE!)',... + 'HandleVisibility','callback',... + 'Tag','figure1',... + 'CloseRequestFcn',@(src,hEvent) figure1_CloseRequestFcn(this,src,hEvent)); + + %WindowState not available in all versions + if isprop(this.guiHandle,'WindowState') + set(this.guiHandle,'WindowState','maximized'); + end + + pParetoSlice = uipanel(... + 'Parent',this.guiHandle,... + 'Units','normalized',... + 'Title','Current plan Slice',... + 'BackgroundColor',matRad_cfg.gui.backgroundColor,... + 'Tag','PlanSlicePanel',... + 'Clipping','off',... + 'Position',[0.005 0.5 0.5 0.49],... + 'FontName',matRad_cfg.gui.fontName,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontWeight',matRad_cfg.gui.fontWeight); + + + pSliderWidget = uipanel( ... + 'Parent',this.guiHandle,... + 'Units','normalized',... + 'BackgroundColor',matRad_cfg.gui.backgroundColor,... + 'Tag','SliderPanel',... + 'Clipping','off',... + 'Position',[0.005 0.005 0.5 0.49],... + 'FontName',matRad_cfg.gui.fontName,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontWeight',matRad_cfg.gui.fontWeight); + pSliderWidget.Scrollable = "on"; + + + pButtonWidget = uipanel(... + 'Parent',this.guiHandle,... + 'Units','normalized',... + 'BackgroundColor',matRad_cfg.gui.backgroundColor,... + 'Tag','Additional Button Panel',... + 'Clipping','off',... + 'Position',[0.51 0.395 0.15 0.1],... + 'FontName',matRad_cfg.gui.fontName,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontWeight',matRad_cfg.gui.fontWeight); + + pParetoDVHWidget = uipanel(... + 'Parent',this.guiHandle,... + 'Units','normalized',... + 'Title','DVH ',... + 'BackgroundColor',matRad_cfg.gui.backgroundColor,... + 'Tag','DVHPanel',... + 'Clipping','off',... + 'Position',[0.51 0.5 0.45 0.49],... + 'FontName',matRad_cfg.gui.fontName,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontWeight',matRad_cfg.gui.fontWeight); + + + fAll = evalin('base','retStruct.finds'); + wAll = evalin('base','retStruct.weights'); + assignin('base','ParetoHelperObject',matRad_ParetoData(fAll,wAll)) + %this.optiProb = evalin('base','retStruct.optiProb'); + %this.helperObject = matRad_ParetoData(fAll,wAll); + + + + this.ParetoSliceWidget = matRad_ParetoSliceWidget(pParetoSlice); + this.ParetoSliderWidget = matRad_ParetoSliderWidget(pSliderWidget,this.ParetoSliceWidget); + this.ParetoDVHWidget = matRad_ParetoDVHWidget(pParetoDVHWidget); + this.ParetoNavigationButtonWidget = matRad_ParetoNavigationButtonWidget(pButtonWidget,this.ParetoSliderWidget,this.ParetoDVHWidget); + end + + function figure1_CloseRequestFcn(this,hObject, ~) + matRad_cfg = MatRad_Config.instance(); + set(0,'DefaultUicontrolBackgroundColor',matRad_cfg.gui.backgroundColor); + selection = questdlg('Do you really want to close the matRad ParetoInterface?',... + 'Close matRad',... + 'Yes','No','Yes'); + + switch selection + case 'Yes' + delete(hObject); + delete(this); + case 'No' + return + end + end + end +end \ No newline at end of file diff --git a/matRad/gui/pareto/matRad_ParetoNavigationButtonWidget.m b/matRad/gui/pareto/matRad_ParetoNavigationButtonWidget.m new file mode 100644 index 000000000..23e68623b --- /dev/null +++ b/matRad/gui/pareto/matRad_ParetoNavigationButtonWidget.m @@ -0,0 +1,127 @@ +classdef matRad_ParetoNavigationButtonWidget < matRad_Widget + % matRad_ParetoNavigationButtonWidget class to generate GUI widget + % to display buttons with additional functionality for Pareto surface + % navigation + % + % References + % - + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + sliderWidgetHandle + dvhWidgetHandle + end + + methods + function this = matRad_ParetoNavigationButtonWidget(handleParent,sliderWidgetHandle,dvhWidgetHandle) + + this = this@matRad_Widget(handleParent); + this.sliderWidgetHandle = sliderWidgetHandle; + this.dvhWidgetHandle = dvhWidgetHandle; + end + end + + methods (Access = protected) + function this = createLayout(this) + ButtonPanel= this.widgetHandle; + matRad_cfg = MatRad_Config.instance(); + %About button + DVHButton = uicontrol(... + 'Parent',ButtonPanel,... + 'Units','normalized',... + 'String','Show DVH',... + 'Tooltip', 'Show the DVH for the current plan',... + 'Position',[0.1 0.55 0.35 0.25],... + 'BackgroundColor',matRad_cfg.gui.elementColor,... + 'ForegroundColor',matRad_cfg.gui.textColor,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontName',matRad_cfg.gui.fontName,... + 'FontWeight',matRad_cfg.gui.fontWeight,... + 'Tag','btn_export',... + 'Callback',@this.DVHButton_callback); + + exportButton = uicontrol(... + 'Parent',ButtonPanel,... + 'Units','normalized',... + 'String','Export',... + 'Tooltip', 'Export the current plan to the main GUI',... + 'UserData',[],... + 'Position',[0.55 0.55 0.35 0.25],... + 'BackgroundColor',matRad_cfg.gui.elementColor,... + 'ForegroundColor',matRad_cfg.gui.textColor,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontName',matRad_cfg.gui.fontName,... + 'FontWeight',matRad_cfg.gui.fontWeight,... + 'Tag','btn_export',... + 'Callback',@this.ExportButton_callback); + + resetButton = uicontrol(... + 'Parent',ButtonPanel,... + 'Units','normalized',... + 'String','Reset constraints',... + 'Tooltip', 'Reset all slider constraints',... + 'UserData',[],... + 'Position',[0.25 0.15 0.5 0.25],... + 'BackgroundColor',matRad_cfg.gui.elementColor,... + 'ForegroundColor',matRad_cfg.gui.textColor,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontName',matRad_cfg.gui.fontName,... + 'FontWeight',matRad_cfg.gui.fontWeight,... + 'Tag','btn_export',... + 'Callback',@this.ResetConstraintButton_callback); + end + + + function DVHButton_callback(this,~,~) + ParetoHelperObject = evalin('base','ParetoHelperObject'); + %Shows the DVH for the current plan + resultGUI = matRad_calcCubes(ParetoHelperObject.currentWeights,evalin('base','dij')); + dvh = matRad_calcDVH(evalin('base','cst'),resultGUI.physicalDose,'cum'); + + matRad_showDVH(this.dvhWidgetHandle.DVHPlotAxes,... + dvh,evalin('base','cst'),evalin('base','pln'),ParetoHelperObject.linestyle); + if ParetoHelperObject.linestyle < 4 + ParetoHelperObject.linestyle = ParetoHelperObject.linestyle + 1; + else + ParetoHelperObject.linestyle = 1; + end + % + + end + + function ExportButton_callback(this,~,~) + ParetoHelperObject = evalin('base','ParetoHelperObject'); + %Export the current plan to the matRadGUi for better inspection + resultGUI = matRad_calcCubes(ParetoHelperObject.currentWeights,evalin('base','dij')); + assignin('base',"resultGUI",resultGUI); + matRadGUI; + + end + + function ResetConstraintButton_callback(this,~,~) + ParetoHelperObject = evalin('base','ParetoHelperObject'); + ParetoHelperObject.releaseObjectiveBounds(); + + for i = 1:numel(this.sliderWidgetHandle.objectiveSliders) + set(this.sliderWidgetHandle.objectiveSliders{i},'Min',ParetoHelperObject.lowboundSliderInit(i)); + set(this.sliderWidgetHandle.objectiveSliders{i},'Max',ParetoHelperObject.initialUpbound(i)); + end + end + + + end + + +end diff --git a/matRad/gui/pareto/matRad_ParetoSliceWidget.m b/matRad/gui/pareto/matRad_ParetoSliceWidget.m new file mode 100644 index 000000000..071a11ec9 --- /dev/null +++ b/matRad/gui/pareto/matRad_ParetoSliceWidget.m @@ -0,0 +1,52 @@ +classdef matRad_ParetoSliceWidget < matRad_Widget + % matRad_ParetoSliceWidget class + % + % + % References + % - + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + DosePlotAxes + end + + methods + function this = matRad_ParetoSliceWidget(handleParent) + this = this@matRad_Widget(handleParent); + this.DosePlotAxes = axes(this.widgetHandle,'Position',[0 0 1 1]); + end + + function this = initialize(this) + end + + + end + + methods (Access = protected) + function this = createLayout(this) + + parent = this.widgetHandle; + + matRad_cfg = MatRad_Config.instance(); + + + this.createHandles(); + end + + end + +end + + diff --git a/matRad/gui/pareto/matRad_ParetoSliderWidget.m b/matRad/gui/pareto/matRad_ParetoSliderWidget.m new file mode 100644 index 000000000..07d7333c4 --- /dev/null +++ b/matRad/gui/pareto/matRad_ParetoSliderWidget.m @@ -0,0 +1,466 @@ +classdef matRad_ParetoSliderWidget < matRad_Widget + % matRad_ParetoSliderWidget class to generate GUI widget to move sliders for + % Pareto surface navigation + % + % References + % - + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2020 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + paretoSliceWidgetHandle + objectiveSliders + fixationButtons + helperObject + optiProb + end + + methods + function this = matRad_ParetoSliderWidget(handleParent,paretoSliceWidgetHandle) + + matRad_cfg = MatRad_Config.instance(); + this = this@matRad_Widget(handleParent); + + this.paretoSliceWidgetHandle = paretoSliceWidgetHandle; + this.objectiveSliders = {}; + %TODO: Best way to pass objective function values? + %fAll = evalin('base','retStruct.finds'); + %wAll = evalin('base','retStruct.weights'); + %this.optiProb = evalin('base','retStruct.optiProb'); + %this.helperObject = matRad_ParetoData(fAll,wAll); + this.initialize(); + this.plotSlice(evalin('base','ParetoHelperObject.currentWeights')); + end + + end + + methods (Access = protected) + function this = createLayout(this) + h1 = this.widgetHandle; + + matRad_cfg = MatRad_Config.instance(); + % handle environment + switch matRad_cfg.env + case 'MATLAB' + set(h1,'SizeChangedFcn',@(hObject,eventdata) widget_SizeChangedFcn(this,hObject,eventdata)); + case 'OCTAVE' + % this creates an infinite loop in octave + %set(h1,'SizeChangedFcn',@(hObject,eventdata) widget_SizeChangedFcn(this,hObject,eventdata)); + end + + this.createHandles(); + + end + + function this = doUpdate(this,evt) + doUpdate = true; + + if nargin == 2 + %At pln changes and at cst/cst (for Isocenter and new settings) + %we need to update + doUpdate = this.checkUpdateNecessary({'cst'},evt); + end + + if doUpdate + if evalin('base','exist(''ct'')') && evalin('base','exist(''cst'')') + generateSliderTable(this, evalin('base','cst')); + else + delete(get(this.widgetHandle,'Children')); + end + end + + end + end + + methods(Access = private) + + function cst = generateSliderTable(this,cst) + matRad_cfg = MatRad_Config.instance(); + ParetoHelperObject = evalin('base','ParetoHelperObject'); + %cst = updateStructureTable(this,cst); + handles = this.handles; + SliderPanel = this.widgetHandle; + + SliderPanelPos = get(SliderPanel,'Position'); + SliderPanelPosUnit = get(SliderPanel,'Units'); + set(SliderPanel,'Units','pixels'); + SliderPanelPosPix = get(SliderPanel,'Position'); + set(SliderPanel,'Units',SliderPanelPosUnit); + aspectRatio = SliderPanelPosPix(3) / SliderPanelPosPix(4); + + %Parameters for line height + objHeight = 0.055;% 22; + lineHeight = 0.1; %25; %Height of a table line + yTopSep = 0.12;%40; %Separation of the first line from the top + %tableViewHeight = SliderPanelPos(4) - yTopSep; %Full height of the view + tableViewHeight = 1 - yTopSep; + + %Widths of the fields + buttonW = objHeight*1.5 / aspectRatio; % Make button squared + nameW = 3*buttonW;%60; + sliderW = 5*buttonW; + functionW = 5*buttonW;%120; + paramTitleW = 4*buttonW;%120; + paramW = 1*buttonW;%30; + fieldSep = 0.25*buttonW; %Separation between fields horizontally + + %disp(num2str(sliderPos)); + SliderPanelChildren = get(SliderPanel,'Children'); + cstVertTableScroll = findobj(SliderPanelChildren,'Tag','VerticalSlider'); + if isempty(cstVertTableScroll) + sliderPos = 0; + else + sliderPos = get(cstVertTableScroll,'Max') - get(cstVertTableScroll,'Value'); + end + + ypos = @(c) tableViewHeight - c*lineHeight + sliderPos; + + delete(SliderPanelChildren); + + + %Creates a dummy axis to allow for the use of textboxes instead of uicontrol to be able to use the (la)tex interpreter + tmpAxes = axes('Parent',SliderPanel,'units','normalized','position',[0 0 1 1],'visible','off', 'FontSize',8); + + organTypes = {'OAR', 'TARGET','IGNORED'}; + + %Get all Classes & classNames + classNames = matRad_getObjectivesAndConstraints(); + + numOfObjectives = 0; + for i = 1:size(cst,1) + for j = 1:numel(cst{i,6}) + if isa(cst{i,6}{j},'DoseObjectives.matRad_DoseObjective') + numOfObjectives = numOfObjectives + 1; + end + end + end + %line + % step count + cnt = 0; + + newline = '\n'; + + %Setup Headlines + xPos = 0.01; %5 + + h2 = uicontrol(SliderPanel,'Style','text', ... + 'String','VOI name', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) nameW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'Tooltip','Name of the structure with objective/constraint', ... + 'BackgroundColor',matRad_cfg.gui.backgroundColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor); + + tmp_pos = get(h2,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + h7 = uicontrol(SliderPanel,'Style','text', ... + 'String','Function slider', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) functionW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'Tooltip','Objective/Constraint function type', ... + 'BackgroundColor',matRad_cfg.gui.backgroundColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor); + + tmp_pos = get(h7,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + h3 = uicontrol(SliderPanel,'Style','text', ... + 'String','Fix slider', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) buttonW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'Tooltip','Button to fixate a slider', ... + 'BackgroundColor',matRad_cfg.gui.backgroundColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor); + + tmp_pos = get(h3,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + h6 = uicontrol(SliderPanel,'Style','text', ... + 'String','Function', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) sliderW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'Tooltip','Objective/Constraint function type', ... + 'BackgroundColor',matRad_cfg.gui.backgroundColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor); + + tmp_pos = get(h6,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + h8 = uicontrol(SliderPanel,'Style','text', ... + 'String','| Parameters', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) paramTitleW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'Tooltip','List of parameters', ... + 'HorizontalAlignment','left', ... + 'BackgroundColor',matRad_cfg.gui.backgroundColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor); + + tmp_pos = get(h8,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + cnt = cnt + 1; + %} + %Create Objectives / Constraints controls + for i = 1:size(cst,1) + if strcmp(cst(i,3),'IGNORED')~=1 + %Compatibility Layer for old objective format + if isstruct(cst{i,6}) + cst{i,6} = num2cell(arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct,cst{i,6})); + end + for j=1:numel(cst{i,6}) + %TODO: Add check if constraint or objective + obj = cst{i,6}{j}; + + %Convert to class if not + if ~isa(obj,'matRad_DoseOptimizationFunction') + try + obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); + catch ME + this.showWarning('Objective/Constraint not valid!\n%s',ME.message) + continue; + end + end + + if isa(obj,'DoseObjectives.matRad_DoseObjective') + xPos = 0.01; + h = uicontrol(SliderPanel','Style','edit', ... + 'String',cst{i,2}, ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) nameW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Tooltip','Name',... + 'Enable','inactive',... %Disable editing of name atm + 'UserData',[i,2]); + + tmp_pos = get(h,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + %Slider + this.objectiveSliders{cnt} = uicontrol(SliderPanel,'Units','normalized',... + 'String','Slider',... + 'Tooltip','Adjust objective penalty',... + 'Position',[xPos ypos(cnt) sliderW objHeight], ... + 'Style','slider',... + 'BusyAction','cancel',... + 'Interruptible','off',... + 'BackgroundColor',matRad_cfg.gui.elementColor,... + 'ForegroundColor',matRad_cfg.gui.textColor,... + 'FontSize',matRad_cfg.gui.fontSize,... + 'FontName',matRad_cfg.gui.fontName,... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'Callback',@(hObject,eventdata)ObjectiveSlider_callback(this,hObject,eventdata,cnt)); + + this.objectiveSliders{cnt}.Value = ParetoHelperObject.currentPoint(cnt); + + tmp_pos = get(this.objectiveSliders{cnt},'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + this.fixationButtons{cnt} = uicontrol(SliderPanel,'Style','pushbutton', ... + 'String','fix', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) buttonW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Tooltip','Remove Objective/Constraint',... + 'Callback',@(hObject,eventdata)FixButton_callback(this,hObject,eventdata,cnt)); + tmp_pos = get(this.fixationButtons{cnt},'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + h = uicontrol(SliderPanel,'Style','text', ... + 'String',obj.name, ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) functionW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Tooltip','Objective/Constraint',... + 'UserData',{[i,j],classNames(1,:)}); + + tmp_pos = get(h,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + + + + for p = 1:numel(obj.parameterNames) + % h = text('Parent',tmpAxes,'String',['| ' obj.parameterNames{p} ':'],'VerticalAlignment','middle','Units','normalized','Position',[xPos ypos(cnt)+lineHeight/2],'Interpreter','tex','FontWeight','normal',... + % 'FontSize',get(SliderPanel,'FontSize'),'FontName',get(SliderPanel,'FontName'),'FontUnits',get(SliderPanel,'FontUnits'),'FontWeight','normal');%[xPos ypos(cnt) 100 objHeight]); + % there is no fontsize for SliderPanel + h = text('Parent',tmpAxes, ... + 'String',['| ' obj.parameterNames{p} ':'], ... + 'VerticalAlignment','middle', ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt)+lineHeight/2], ... + 'Interpreter','tex', ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.backgroundColor, ... + 'Color',matRad_cfg.gui.textColor);%[xPos ypos(cnt) 100 objHeight]); + + tmp_pos = get(h,'Extent'); + xPos = xPos + tmp_pos(3) + fieldSep; + %h = annotation(SliderPanel,'textbox','String',obj.parameters{1,p},'Units','pix','Position', [xPos ypos(cnt) 100 objHeight],'Interpreter','Tex'); + + %Check if we have a cell and therefore a parameter list + if iscell(obj.parameterTypes{p}) + h = uicontrol(SliderPanel,'Style','popupmenu', ... + 'String',obj.parameterTypes{p}', ... + 'Value',obj.parameters{p}, ... + 'Tooltip',obj.parameterNames{p}, ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) paramW*2 objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'UserData',[i,j,p]); + else + h = uicontrol(SliderPanel,'Style','text', ... + 'String',num2str(obj.parameters{p}), ... + 'Tag','FunctionValueSlider',... + 'Tooltip',obj.parameterNames{p}, ... + 'Units','normalized', ... + 'Position',[xPos ypos(cnt) paramW objHeight], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'UserData',[i,j,p]); + end + + tmp_pos = get(h,'Position'); + xPos = xPos + tmp_pos(3) + fieldSep; + end + + + cnt = cnt +1; + end + end + end + end + %Calculate Scrollbar + lastPos = ypos(cnt); + firstPos = ypos(0); + tableHeight = abs(firstPos - lastPos); + + exceedFac = tableHeight / tableViewHeight; + if exceedFac > 1 + sliderFac = exceedFac - 1; + uicontrol(SliderPanel,'Style','slider', ... + 'Units','normalized', ... + 'Tag','VerticalSlider',... + 'Position',[0.975 0 0.025 1], ... + 'FontSize',matRad_cfg.gui.fontSize, ... + 'FontName',matRad_cfg.gui.fontName, ... + 'FontWeight',matRad_cfg.gui.fontWeight, ... + 'BackgroundColor',matRad_cfg.gui.elementColor, ... + 'ForegroundColor',matRad_cfg.gui.textColor, ... + 'Min',0,'Max',ceil(sliderFac)*tableViewHeight, ... + 'SliderStep',[lineHeight tableViewHeight] ./ (ceil(sliderFac)*tableViewHeight), ... + 'Value',ceil(sliderFac)*tableViewHeight - sliderPos, ... + 'Callback', @(hObject,eventdata)cstTableSlider_Callback(this,hObject,eventdata)); + end + this.handles = handles; + end + + function ObjectiveSlider_callback(this,slider,~,idx) + % + ParetoHelperObject = evalin('base','ParetoHelperObject'); + v = matRad_paretoSurfaceNavigation(ParetoHelperObject.allPoints',ParetoHelperObject.currentPoint,slider.Value,idx,ParetoHelperObject.availableUpbound); + + if numel(v) == 0 + %navigation algorithm didnt find a new point + for i = 1:numel(this.objectiveSliders) + set(this.objectiveSliders{i},'Value',ParetoHelperObject.currentPoint(i)); + end + else + wNew = ParetoHelperObject.allWeights*v; + + fNew = ParetoHelperObject.allPoints'*v;% evaluate function(strictly necessary?) + %fNew = this.optiProb.normalizeObjectives(fNew'); + + + ParetoHelperObject.currentPoint = fNew; + ParetoHelperObject.currentWeights = wNew; + + %Open: Current + for i = 1:numel(this.objectiveSliders) + set(this.objectiveSliders{i},'Value',fNew(i)); + end + + this.plotSlice(wNew); + end + end + + function plotSlice(this,w) + slice = round(evalin('base','pln.propStf.isoCenter(1,3)./ct.resolution.z')); + cubes = matRad_calcFastCubes(w,evalin('base','dij'),evalin('base','pln')); + + matRad_plotSliceWrapper(this.paretoSliceWidgetHandle.DosePlotAxes,evalin('base','ct'),evalin('base','cst'),1,cubes,3,slice,[],[],[],[],[],[],[],[],[],'LineWidth',2); + end + + function FixButton_callback(this,~,~,idx) + + ParetoHelperObject = evalin('base','ParetoHelperObject'); + [lb,ub] = ParetoHelperObject.restrictObjective(idx,this.objectiveSliders{idx}.Value); %update refObjects bounds + + + for i = 1:numel(this.objectiveSliders) + set(this.objectiveSliders{i},'Min',lb(i)); + set(this.objectiveSliders{i},'Max',ub(i)); + end + end + + + % --- Executes when the widget is resized. + function widget_SizeChangedFcn(this,hObject, eventdata) + % hObject handle to h1 (see GCBO) + % eventdata reserved - to be defined in a future version of MATLAB + % handles structure with handles and user data (see GUIDATA) + try + generateSliderTable(this,evalin('base','cst')); + catch + end + end + + % --- Executes on slider movement. + function cstTableSlider_Callback(this,~,~) + % hObject handle to cstTableSlider (see GCBO) + % eventdata reserved - to be defined in a future version of MATLAB + % handles structure with handles and user data (see GUIDATA) + + % Hints: get(hObject,'Value') returns position of slider + % get(hObject,'Min') and get(hObject,'Max') to determine range of slider + try + generateSliderTable(this,evalin('base','cst')); + catch + end + end + end +end diff --git a/matRad/matRad_fluenceOptimization.m b/matRad/matRad_fluenceOptimization.m old mode 100644 new mode 100755 index e5f9951d9..510a36e30 --- a/matRad/matRad_fluenceOptimization.m +++ b/matRad/matRad_fluenceOptimization.m @@ -34,271 +34,10 @@ matRad_cfg = MatRad_Config.instance(); -% consider VOI priorities -cst = matRad_setOverlapPriorities(cst); - -% check & adjust objectives and constraints internally for fractionation -for i = 1:size(cst,1) - %Compatibility Layer for old objective format - if isstruct(cst{i,6}) - cst{i,6} = arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct,cst{i,6},'UniformOutput',false); - end - for j = 1:numel(cst{i,6}) - - obj = cst{i,6}{j}; - - %In case it is a default saved struct, convert to object - %Also intrinsically checks that we have a valid optimization - %objective or constraint function in the end - if ~isa(obj,'matRad_DoseOptimizationFunction') - try - obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); - catch - matRad_cfg.dispError('cst{%d,6}{%d} is not a valid Objective/constraint! Remove or Replace and try again!',i,j); - end - end - - obj = obj.setDoseParameters(obj.getDoseParameters()/pln.numOfFractions); - - cst{i,6}{j} = obj; - end -end - -% 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); - -% Get rid of voxels that are not interesting for the optimization problem -if ~isfield(pln.propOpt, 'clearUnusedVoxels') - pln.propOpt.clearUnusedVoxels = matRad_cfg.defaults.propOpt.clearUnusedVoxels; -end - -if pln.propOpt.clearUnusedVoxels - dij = matRad_clearUnusedVoxelsFromDij(cst, dij); -end - - -% find target indices and described dose(s) for weight vector -% initialization -V = []; -doseTarget = []; -ixTarget = []; - -for i = 1:size(cst,1) - if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6}) - V = [V;cst{i,4}{1}]; - - %Iterate through objectives/constraints - fDoses = []; - for fObjCell = cst{i,6} - dParams = fObjCell{1}.getDoseParameters(); - %Don't care for Inf constraints - dParams = dParams(isfinite(dParams)); - %Add do dose list - fDoses = [fDoses dParams]; - end - - doseTarget = [doseTarget fDoses]; - ixTarget = [ixTarget i*ones(1,length(fDoses))]; - end -end -[doseTarget,i] = max(doseTarget); -ixTarget = ixTarget(i); -wOnes = ones(dij.totalNumOfBixels,1); - -% calculate initial beam intensities wInit -matRad_cfg.dispInfo('Estimating initial weights... '); - if exist('wInit','var') - %do nothing as wInit was passed to the function - matRad_cfg.dispInfo('chosen provided wInit!\n'); - - % Write ixDose which is needed for the optimizer - if pln.bioParam.bioOpt - dij.ixDose = dij.bx~=0; - - %pre-calculations - dij.gamma = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); - dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose)); - end - -elseif strcmp(pln.bioParam.model,'constRBE') && strcmp(pln.radiationMode,'protons') - % check if a constant RBE is defined - if not use 1.1 - if ~isfield(dij,'RBE') - dij.RBE = 1.1; - end - - doseTmp = dij.physicalDose{1}*wOnes; - bixelWeight = (doseTarget)/(dij.RBE * mean(doseTmp(V))); - wInit = wOnes * bixelWeight; - matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); - -elseif pln.bioParam.bioOpt - % retrieve photon LQM parameter - [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels); - checkAxBx = cellfun(@(ax1,bx1,ax2,bx2) isequal(ax1(ax1~=0),ax2(ax1~=0)) && isequal(bx1(bx1~=0),bx2(bx1~=0)),dij.ax,dij.bx,ax,bx); - if ~all(checkAxBx) - matRad_cfg.dispError('Inconsistent biological parameters in dij.ax and/or dij.bx - please recalculate dose influence matrix before optimization!\n'); - end - - - - for i = 1:size(cst,1) - - for j = 1:size(cst{i,6},2) - % check if prescribed doses are in a valid domain - if any(cst{i,6}{j}.getDoseParameters() > 5) && isequal(cst{i,3},'TARGET') - matRad_cfg.dispError('Reference dose > 5 Gy[RBE] for target. Biological optimization outside the valid domain of the base data. Reduce dose prescription or use more fractions.\n'); - end - - end - end - - for s = 1:numel(dij.bx) - dij.ixDose{s} = dij.bx{s}~=0; - end - - if isequal(pln.bioParam.quantityOpt,'effect') - - effectTarget = cst{ixTarget,5}.alphaX * doseTarget + cst{ixTarget,5}.betaX * doseTarget^2; - aTmp = dij.mAlphaDose{1}*wOnes; - bTmp = dij.mSqrtBetaDose{1} * wOnes; - p = sum(aTmp(V)) / sum(bTmp(V).^2); - q = -(effectTarget * length(V)) / sum(bTmp(V).^2); - - wInit = -(p/2) + sqrt((p^2)/4 -q) * wOnes; - - elseif isequal(pln.bioParam.quantityOpt,'RBExD') - - %pre-calculations - for s = 1:numel(dij.ixDose) - dij.gamma{s} = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); - dij.gamma{s}(dij.ixDose{s}) = dij.ax{s}(dij.ixDose{s})./(2*dij.bx{s}(dij.ixDose{s})); - end - - - % calculate current effect in target - aTmp = dij.mAlphaDose{1}*wOnes; - bTmp = dij.mSqrtBetaDose{1} * wOnes; - doseTmp = dij.physicalDose{1}*wOnes; - - CurrEffectTarget = aTmp(V) + bTmp(V).^2; - % ensure a underestimated biological effective dose - TolEstBio = 1.2; - % calculate maximal RBE in target - 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 - - matRad_cfg.dispInfo('chosen weights adapted to biological dose calculation!\n'); -else - doseTmp = dij.physicalDose{1}*wOnes; - bixelWeight = (doseTarget)/mean(doseTmp(V)); - wInit = wOnes * bixelWeight; - matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); -end - - -%% calculate probabilistic quantities for probabilistic optimization if at least -% one robust objective is defined - -%Check how to use 4D data -if isfield(pln,'propOpt') && isfield(pln.propOpt,'scen4D') - scen4D = pln.propOpt.scen4D; -else - scen4D = 1; %Use only first 4D scenario for optimization -end - -%If "all" provided, use all scenarios -if isequal(scen4D,'all') - scen4D = 1:size(dij.physicalDose,1); -end - -linIxDIJ = find(~cellfun(@isempty,dij.physicalDose(scen4D,:,:)))'; - -%Only select the indexes of the nominal ct Scenarios -linIxDIJ_nominalCT = find(~cellfun(@isempty,dij.physicalDose(scen4D,1,1)))'; - -FLAG_CALC_PROB = false; -FLAG_ROB_OPT = false; - - -for i = 1:size(cst,1) - for j = 1:numel(cst{i,6}) - if strcmp(cst{i,6}{j}.robustness,'PROB') && numel(linIxDIJ) > 1 - FLAG_CALC_PROB = true; - end - if ~strcmp(cst{i,6}{j}.robustness,'none') && numel(linIxDIJ) > 1 - FLAG_ROB_OPT = true; - end - end -end - -if FLAG_CALC_PROB - [dij] = matRad_calculateProbabilisticQuantities(dij,cst,pln); -end - - -% set optimization options -if ~FLAG_ROB_OPT || FLAG_CALC_PROB % if multiple robust objectives are defined for one structure then remove FLAG_CALC_PROB from the if clause - ixForOpt = scen4D; -else - ixForOpt = linIxDIJ; -end - -switch pln.bioParam.quantityOpt - case 'effect' - backProjection = matRad_EffectProjection; - case 'RBExD' - %Capture special case of constant RBE - if strcmp(pln.bioParam.model,'constRBE') - backProjection = matRad_ConstantRBEProjection; - else - backProjection = matRad_VariableRBEProjection; - end - case 'physicalDose' - backProjection = matRad_DoseProjection; - otherwise - warning(['Did not recognize bioloigcal setting ''' pln.probOpt.bioOptimization '''!\nUsing physical dose optimization!']); - backProjection = matRad_DoseProjection; -end - -%Give scenarios used for optimization -backProjection.scenarios = ixForOpt; -backProjection.scenarioProb = pln.multScen.scenProb; -backProjection.nominalCtScenarios = linIxDIJ_nominalCT; -%backProjection.scenDim = pln.multScen - -optiProb = matRad_OptimizationProblem(backProjection); -optiProb.quantityOpt = pln.bioParam.quantityOpt; -if isfield(pln,'propOpt') && isfield(pln.propOpt,'useLogSumExpForRobOpt') - optiProb.useLogSumExpForRobOpt = pln.propOpt.useLogSumExpForRobOpt; -end - -%Get Bounds -if ~isfield(pln.propOpt,'boundMU') - pln.propOpt.boundMU = false; -end - -if pln.propOpt.boundMU - if (isfield(dij,'minMU') || isfield(dij,'maxMU')) && ~isfield(dij,'numParticlesPerMU') - matRad_cfg.dispWarning('Requested MU bounds but number of particles per MU not set! Bounds will not be enforced and standard [0,Inf] will be used instead!'); - elseif ~isfield(dij,'minMU') && ~isfield(dij,'maxMU') - matRad_cfg.dispWarning('Requested MU bounds but machine bounds not defined in dij.minMU & dij.maxMU! Bounds will not be enforced and standard [0,Inf] will be used instead!'); - else - if isfield(dij,'minMU') - optiProb.minimumW = dij.numParticlesPerMU .* dij.minMU / 1e6; - matRad_cfg.dispInfo('Using lower MU bounds provided in dij!\n') - end - - if isfield(dij,'maxMU') - optiProb.maximumW = dij.numParticlesPerMU .* dij.maxMU / 1e6; - matRad_cfg.dispInfo('Using upper MU bounds provided in dij!\n') - end - end + [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln,wInit); else - matRad_cfg.dispInfo('Using standard MU bounds of [0,Inf]!\n') + [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln); end if ~isfield(pln.propOpt,'optimizer') @@ -328,9 +67,10 @@ resultGUI.wUnsequenced = wOpt; resultGUI.usedOptimizer = optimizer; resultGUI.info = info; +resultGUI.optiProb = optiProb; %Robust quantities -if FLAG_ROB_OPT || numel(ixForOpt) > 1 +if FLAG_ROB_OPT || numel(optiProb.BP.scenarios) > 1 Cnt = 1; for i = find(~cellfun(@isempty,dij.physicalDose))' tmpResultGUI = matRad_calcCubes(wOpt,dij,i); @@ -341,3 +81,5 @@ % unblock mex files clear mex + +end \ No newline at end of file diff --git a/matRad/optimization/+DoseConstraints/matRad_DoseConstraintFromObjective.m b/matRad/optimization/+DoseConstraints/matRad_DoseConstraintFromObjective.m index 2a3ee0473..6be2c6472 100644 --- a/matRad/optimization/+DoseConstraints/matRad_DoseConstraintFromObjective.m +++ b/matRad/optimization/+DoseConstraints/matRad_DoseConstraintFromObjective.m @@ -19,17 +19,17 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% properties (Constant) name = 'Objective Constraint'; - parameterTypes = {'objFunc','scalar'}; - parameterNames = {'f^{max}','slackParameter'}; + parameterTypes = {'objFunc'}; + parameterNames = {'f^{max}'}; end properties objective; - parameters = {1e-5, 1e-3}; + parameters = {1e-5}; end methods (Access = public) - function this = matRad_DoseConstraintFromObjective(objective,maxObj,slackParameter) + function this = matRad_DoseConstraintFromObjective(objective,maxObj) %check if objective is a struct and a DoseObjective or Constraint (for init from constraint) if isstruct(objective) && ~isempty(strfind(objective.className,'DoseObjectives')) @@ -48,12 +48,8 @@ if ~initFromStruct - - if nargin == 3 && isscalar(slackParameter) - this.parameters{2} = slackParameter; - end - if nargin >= 2 && isscalar(maxObj) + if nargin == 2 && isscalar(maxObj) this.parameters{1} = maxObj; end @@ -71,9 +67,10 @@ end function cu = upperBounds(this,n) - cu = this.parameters{1}+this.slackParameter; + cu = this.parameters{1}; %cu = [Inf; this.parameters{2}]; end + function cl = lowerBounds(this,n) cl = 0; end diff --git a/matRad/optimization/+DoseConstraints/matRad_MinMaxDose.m b/matRad/optimization/+DoseConstraints/matRad_MinMaxDose.m index b6ba18611..e5d35e613 100644 --- a/matRad/optimization/+DoseConstraints/matRad_MinMaxDose.m +++ b/matRad/optimization/+DoseConstraints/matRad_MinMaxDose.m @@ -173,41 +173,41 @@ function cDose = computeDoseConstraintFunctionLogSumExp(this,dose) dose_min = min(dose); dose_max = max(dose); - + modEpsilon = (this.epsilon/3)*dose_max; %Validate parameters if this.parameters{1} <= 0 && isinf(this.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf) cDose = []; elseif this.parameters{2} == Inf %Only min dose - cDose = dose_min - this.epsilon * log( sum(exp((dose_min - dose)/this.epsilon))); + cDose = dose_min - modEpsilon * log( sum(exp((dose_min - dose)/modEpsilon))); elseif this.parameters{1} <= 0 %Only max dose - cDose = dose_max + this.epsilon * log( sum(exp((dose - dose_max)/this.epsilon))); + cDose = dose_max + modEpsilon * log( sum(exp((dose - dose_max)/modEpsilon))); else %both are set sensible - cDose(2,1) = dose_max + this.epsilon * log( sum(exp((dose - dose_max)/this.epsilon))); - cDose(1,1) = dose_min - this.epsilon * log( sum(exp((dose_min - dose)/this.epsilon))); + cDose(2,1) = dose_max + modEpsilon * log( sum(exp((dose - dose_max)/modEpsilon))); + cDose(1,1) = dose_min - modEpsilon * log( sum(exp((dose_min - dose)/modEpsilon))); end end function cDoseJacob = computeDoseConstraintJacobianLogSumExp(this,dose) %Validate parameters + max_dose = max(dose); + modEpsilon = (this.epsilon/3)*max_dose; + %espilon = this.epsilon; if this.parameters{1} <= 0 && isinf(this.parameters{2}) %Constraint doesn't make sense (min = 0 & max = Inf) cDoseJacob = []; elseif this.parameters{2} == Inf %Only min dose - cDoseJacob(:,1) = exp( (min(dose)-dose)/this.epsilon ); + cDoseJacob(:,1) = exp( (min(dose)-dose)/modEpsilon); cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1)); elseif this.parameters{1} <= 0 %Only max dose - cDoseJacob(:,1) = exp( (dose-max(dose))/this.epsilon ); + cDoseJacob(:,1) = exp( (dose-max_dose)/modEpsilon ); cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1)); else %both are set sensible - cDoseJacob(:,1) = exp( (min(dose)-dose)/this.epsilon ); + cDoseJacob(:,1) = exp( (min(dose)-dose)/modEpsilon ); cDoseJacob(:,1) = cDoseJacob(:,1)/sum(cDoseJacob(:,1)); - - cDoseJacob(:,2) = exp( (dose-max(dose))/this.epsilon ); + cDoseJacob(:,2) = exp( (dose-max_dose)/modEpsilon ); cDoseJacob(:,2) = cDoseJacob(:,2)/sum(cDoseJacob(:,2)); end - - + end - %Exact voxel-wise function cDose = computeDoseConstraintFunctionVoxelwise(this,dose) cDose = dose; diff --git a/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m b/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m index 5322c97b6..6a7301ed8 100644 --- a/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m +++ b/matRad/optimization/+DoseObjectives/matRad_DoseObjective.m @@ -28,6 +28,14 @@ rob = {'none','STOCH','PROB','VWWC','VWWC_INV','COWC','OWC'}; %By default, no robustness is available end end + + %These should be abstract methods, however Octave can't parse them. As soon + %as Octave is able to do this, they should be made abstract again + methods %(Static,Abstract) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + error('Function needs to be implemented'); + end + end %These should be abstract methods, however Octave can't parse them. As soon %as Octave is able to do this, they should be made abstract again diff --git a/matRad/optimization/+DoseObjectives/matRad_EUD.m b/matRad/optimization/+DoseObjectives/matRad_EUD.m index 67f28ab25..0e4caebc6 100644 --- a/matRad/optimization/+DoseObjectives/matRad_EUD.m +++ b/matRad/optimization/+DoseObjectives/matRad_EUD.m @@ -99,6 +99,19 @@ error(['EUD computation failed. Reduce exponent to resolve numerical problems.']); end end + + function constr = turnIntoLexicographicConstraint(obj,goal) + objective = DoseObjectives.matRad_EUD(100,obj.parameters{1,1},obj.parameters{1,2}); + constr = DoseConstraints.matRad_DoseConstraintFromObjective(objective,goal); + end + + end + + methods (Static) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue/numOfFractions; + end + end end diff --git a/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m b/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m index bf505662c..1a58fd5a8 100644 --- a/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m +++ b/matRad/optimization/+DoseObjectives/matRad_MaxDVH.m @@ -95,4 +95,9 @@ end end + methods (Static) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue; + end + end end diff --git a/matRad/optimization/+DoseObjectives/matRad_MeanDose.m b/matRad/optimization/+DoseObjectives/matRad_MeanDose.m index 93feeffe2..07a0bde4c 100644 --- a/matRad/optimization/+DoseObjectives/matRad_MeanDose.m +++ b/matRad/optimization/+DoseObjectives/matRad_MeanDose.m @@ -106,6 +106,12 @@ matRad_cfg.dispError('Invalid setting for %s in Mean Dose Objective!',obj.parameterNames{2}); end end + + function constr = turnIntoLexicographicConstraint(obj,goal) + objective = DoseObjectives.matRad_MeanDose(100,obj.parameters{1},obj.parameters{2}); + constr = DoseConstraints.matRad_DoseConstraintFromObjective(objective,goal); + end + end methods (Access = protected) @@ -125,6 +131,11 @@ fDoseGrad = (1/numel(dose))*sign(dose(:)-obj.parameters{1}); end end - + + methods (Static) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue/numOfFractions; + end + end end diff --git a/matRad/optimization/+DoseObjectives/matRad_MinDVH.m b/matRad/optimization/+DoseObjectives/matRad_MinDVH.m index 6ddd4d4f9..e00bf0f32 100644 --- a/matRad/optimization/+DoseObjectives/matRad_MinDVH.m +++ b/matRad/optimization/+DoseObjectives/matRad_MinDVH.m @@ -95,5 +95,12 @@ fDoseGrad = (2/numel(dose))*deviation; end end + + methods (Static) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue/numOfFractions; + end + + end end \ No newline at end of file diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredDeviation.m b/matRad/optimization/+DoseObjectives/matRad_SquaredDeviation.m index 75c84c8ac..66d837dc3 100644 --- a/matRad/optimization/+DoseObjectives/matRad_SquaredDeviation.m +++ b/matRad/optimization/+DoseObjectives/matRad_SquaredDeviation.m @@ -75,6 +75,15 @@ % calculate delta fDoseGrad = 2 * 1/numel(dose) * deviation; end + + function constr = turnIntoLexicographicConstraint(obj,goal) + if goal < 5e-4 + goal = 5e-4*1.03; + end + objective = DoseObjectives.matRad_SquaredDeviation(100,obj.parameters{1}); + constr = DoseConstraints.matRad_DoseConstraintFromObjective(objective,goal); + end + end methods (Static) @@ -82,6 +91,11 @@ rob = DoseObjectives.matRad_DoseObjective.availableRobustness(); rob{end+1} = 'PROB'; %By default, no robustness is available end + + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue/numOfFractions^2; + end + end end diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredOverdosing.m b/matRad/optimization/+DoseObjectives/matRad_SquaredOverdosing.m index abb0571fc..72ae7b714 100644 --- a/matRad/optimization/+DoseObjectives/matRad_SquaredOverdosing.m +++ b/matRad/optimization/+DoseObjectives/matRad_SquaredOverdosing.m @@ -78,6 +78,20 @@ % calculate delta fDoseGrad = 2 * 1/numel(dose) * overdose; end + + function constr = turnIntoLexicographicConstraint(obj,goal) + if goal < 5e-4 + goal = 5e-4*1.03; + end + objective = DoseObjectives.matRad_SquaredOverdosing(100,obj.parameters{1}); + constr = DoseConstraints.matRad_DoseConstraintFromObjective(objective,goal); + end + end + + methods (Static) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue/numOfFractions^2; + end end end diff --git a/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m b/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m index 006295b9a..d10c66f82 100644 --- a/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m +++ b/matRad/optimization/+DoseObjectives/matRad_SquaredUnderdosing.m @@ -78,6 +78,21 @@ % calculate delta fDoseGrad = 2/numel(dose) * underdose; end + + function constr = turnIntoLexicographicConstraint(obj,goal) + if goal < 5e-4 + goal = 5e-4*1.03; + end + objective = DoseObjectives.matRad_SquaredUnderdosing(100,obj.parameters{1}); + constr = DoseConstraints.matRad_DoseConstraintFromObjective(objective,goal); + end + + end + + methods (Static) + function newGoalValue = adaptGoalToFraction(goalValue,numOfFractions) + newGoalValue = goalValue/numOfFractions^2; + end end end diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m index 529e20f3b..f8ea0e1b7 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_OptimizationProblem.m @@ -24,6 +24,11 @@ properties BP + normalizationScheme = struct('scheme','none'); % used in pareto optimization + objectives = {}; %cell array storing all objectives, has to be initialized at the start + constraints = {}; % + objIdx; + constrIdx; quantityOpt = ''; useMaxApprox = 'logsumexp'; %'pnorm'; %'logsumexp'; %'none'; p = 30; %Can be chosen larger (closer to maximum) or smaller (closer to mean). Only tested 20 >= p >= 1 @@ -31,14 +36,22 @@ minimumW = NaN; maximumW = NaN; end + methods - function obj = matRad_OptimizationProblem(backProjection) - obj.BP = backProjection; + function obj = matRad_OptimizationProblem(backProjection,cst) + + obj.BP = backProjection; %needs to be initalized to have access to setBiologicalDosePrescriptions + if nargin == 2 + obj.extractObjectivesAndConstraintsFromCst(cst); + end end %Objective function declaration fVal = matRad_objectiveFunction(optiProb,w,dij,cst) + + %Objective function declaration + fIndv = matRad_objectiveFunctions(optiProb,w,dij,cst) %Objective gradient declaration fGrad = matRad_objectiveGradient(optiProb,w,dij,cst) @@ -81,6 +94,84 @@ matRad_cfg.dispError('Maximum Bounds for Optimization Problem could not be set!'); end end + + function normalizedfVals = normalizeObjectives(optiProb,fVals) + %function to normalize objectives (used for sandwich + %algorithms) + switch optiProb.normalizationScheme.name + case 'none' + %default case no normalization + normalizedfVals = fVals; + case 'UL' + %used to normalize with respect to min and max values + %maybe check that U and L are defined + normalizedfVals = (fVals - optiProb.normalizationScheme.L)./(optiProb.normalizationScheme.U-optiProb.normalizationScheme.L); %might have to check that U and L work! + otherwise + matRad_cfg.dispError('Normalization scheme not known!'); + end + end + + function normalizedGradient = normalizeGradient(optiProb,Gradient,i) + %function to normalize objectives (used for sandwich + %algorithms) + switch optiProb.normalizationScheme.scheme + + case 'none' + %default case no normalization + normalizedGradient = Gradient; + case 'UL' + %used to normalize with respect to min and max values + %maybe check that U and L are defined + normalizedGradient = Gradient./(optiProb.normalizationScheme.U(i)-optiProb.normalizationScheme.L(i)); %might have to check that U and L work! + otherwise + matRad_cfg.dispError('Normalization scheme not known!'); + end + end + + function updatePenalties(optiProb,newPen) + %TODO: Allow grouping of objectives + if numel(optiProb.objectives) ~= numel(newPen) + matRad_cfg.dispError('Number of objectives in optimization Problem not equal to number of new penalties to be set!'); + end + for i=1:numel(newPen) + optiProb.objectives{i}.penalty = newPen(i)*100; %scaling of penalties with factor 100 + end + end + + function extractObjectivesAndConstraintsFromCst(optiProb,cst) + %used to extract objectives from cst and store in cell array as property of optimization Problem + optiProb.objIdx = []; + optiProb.constrIdx = []; + optiProb.objectives = {}; + optiProb.constraints = {}; + + for i = 1:size(cst,1) % loop over cst + + if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) + + for j = 1:numel(cst{i,6}) + %check whether dose objective or constraint + obj = cst{i,6}{j}; + if isstruct(cst{i,6}{j}) + obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); + end + if contains(class(obj),'DoseObjectives') + optiProb.objIdx = [optiProb.objIdx;i,j]; + obj = optiProb.BP.setBiologicalDosePrescriptions(obj,cst{i,5}.alphaX,cst{i,5}.betaX); + optiProb.objectives(end+1) = {obj}; + + elseif contains(class(obj),'DoseConstraints') + optiProb.constrIdx = [optiProb.constrIdx;i,j]; + obj = optiProb.BP.setBiologicalDosePrescriptions(obj,cst{i,5}.alphaX,cst{i,5}.betaX); + optiProb.constraints(end+1) = {obj}; + end + end + end + end + end + + + end methods (Access = protected) @@ -117,6 +208,9 @@ grad = fac * (tmp ./ pNormVal).^(p-1); end - end -end + end + + + +end diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m index c355e1da7..f70388cc4 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintFunctions.m @@ -32,7 +32,6 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % get current dose / effect / RBExDose vector optiProb.BP.compute(dij,w); d = optiProb.BP.GetResult(); @@ -50,96 +49,77 @@ c = []; % compute objective function for every VOI. -for i = 1:size(cst,1) - - % Only take OAR or target VOI. - if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) - - % loop over the number of constraints for the current VOI - for j = 1:numel(cst{i,6}) +for i = 1:size(optiProb.constrIdx,1) - constraint = cst{i,6}{j}; + constraint = optiProb.constraints{i}; + curConIdx = optiProb.constrIdx(i,1); + + robustness = constraint.robustness; - % only perform computations for constraints - if isa(constraint,'DoseConstraints.matRad_DoseConstraint') + switch robustness + case 'none' % if conventional opt: just sum objectives of nominal dose + d_i = d{1}(cst{curConIdx,4}{1}); + c = [c; constraint.computeDoseConstraintFunction(d_i)]; + + case 'PROB' % if prob opt: sum up expectation value of objectives + + d_i = dExp{1}(cst{curConIdx,4}{1}); + c = [c; constraint.computeDoseConstraintFunction(d_i)]; + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end - % rescale dose parameters to biological optimization quantity if required - constraint = optiProb.BP.setBiologicalDosePrescriptions(constraint,cst{i,5}.alphaX,cst{i,5}.betaX); + d_Scen = d_tmp(cst{curConIdx,4}{contourIx},:); - % retrieve the robustness type - robustness = constraint.robustness; + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); - switch robustness - case 'none' % if conventional opt: just sum objectives of nominal dose - d_i = d{1}(cst{i,4}{1}); - c = [c; constraint.computeDoseConstraintFunction(d_i)]; - - case 'PROB' % if prob opt: sum up expectation value of objectives - - d_i = dExp{1}(cst{i,4}{1}); - c = [c; constraint.computeDoseConstraintFunction(d_i)]; - - case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - - d_max = max(d_Scen,[],2); - d_min = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_max; - elseif isequal(cst{i,3},'TARGET') - d_i = d_min; - end - - c = [c; constraint.computeDoseConstraintFunction(d_i)]; - - case 'VWWC_INV' %inverse voxel-wise conformitiy - takes maximum dose in TARGET and minimum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D inverted VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector - if ~exist('d_tmp','var') - d_tmp = [d{:}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - d_max = max(d_Scen,[],2); - d_min = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_min; - elseif isequal(cst{i,3},'TARGET') - d_i = d_max; - end - - c = [c; constraint.computeDoseConstraintFunction(d_i)]; - otherwise - matRad_cfg.dispError('Robustness setting %s not yet supported!',constraint.robustness); + if isequal(cst{curConIdx,3},'OAR') + d_i = d_max; + elseif isequal(cst{curConIdx,3},'TARGET') + d_i = d_min; end + c = [c; constraint.computeDoseConstraintFunction(d_i)]; - end + case 'VWWC_INV' %inverse voxel-wise conformitiy - takes maximum dose in TARGET and minimum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D inverted VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [d{:}]; + end + + d_Scen = d_tmp(cst{curConIdx,4}{contourIx},:); + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{curConIdx,3},'OAR') + d_i = d_min; + elseif isequal(cst{curConIdx,3},'TARGET') + d_i = d_max; + end + + c = [c; constraint.computeDoseConstraintFunction(d_i)]; + otherwise + matRad_cfg.dispError('Robustness setting %s not yet supported!',constraint.robustness); + end - end % if we are a constraint - - end % over all defined constraints & objectives - -end % if structure not empty and oar or target - -end % over all structures +end % if we are a constraint +end + % if structure not empty and oar or target diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m index a411cc5a2..1118acacc 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_constraintJacobian.m @@ -32,7 +32,6 @@ % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % get current dose / effect / RBExDose vector %d = matRad_backProjection(w,dij,optiProb); %d = optiProb.matRad_backProjection(w,dij); @@ -44,6 +43,7 @@ % initialize projection matrices and id containers DoseProjection{1} = sparse([]); +%DoseProjection{1} = []; mAlphaDoseProjection{1} = sparse([]); mSqrtBetaDoseProjection{1} = sparse([]); voxelID = []; @@ -59,154 +59,138 @@ contourScen = fullScen{1}; % compute objective function for every VOI. -for i = 1:size(cst,1) +for i = 1:size(optiProb.constrIdx,1) - % Only take OAR or target VOI. - if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) - - % loop over the number of constraints for the current VOI - for j = 1:numel(cst{i,6}) + constraint = optiProb.constraints{i}; %Get the Optimization Object + curConIdx = optiProb.constrIdx(i,1); + robustness = constraint.robustness; + switch robustness - constraint = cst{i,6}{j}; %Get the Optimization Object - - % only perform computations for constraints - if isa(constraint,'DoseConstraints.matRad_DoseConstraint') + case 'none' % if conventional opt: just sum objectiveectives of nominal dose + d_i = d{1}(cst{curConIdx,4}{1}); + jacobSub = constraint.computeDoseConstraintJacobian(d_i); + + case 'PROB' % if prob opt: sum up expectation value of objectives - % retrieve the robustness type - robustness = constraint.robustness; + d_i = dExp{1}(cst{curConIdx,4}{1}); + jacobSub = constraint.computeDoseConstraintJacobian(d_i); - % rescale dose parameters to biological optimization quantity if required - constraint = optiProb.BP.setBiologicalDosePrescriptions(constraint,cst{i,5}.alphaX,cst{i,5}.betaX); + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end - switch robustness - - case 'none' % if conventional opt: just sum objectiveectives of nominal dose - d_i = d{1}(cst{i,4}{1}); - jacobSub = constraint.computeDoseConstraintJacobian(d_i); - - case 'PROB' % if prob opt: sum up expectation value of objectives - - d_i = dExp{1}(cst{i,4}{1}); - jacobSub = constraint.computeDoseConstraintJacobian(d_i); - - case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - - d_max = max(d_Scen,[],2); - d_min = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_max; - elseif isequal(cst{i,3},'TARGET') - d_i = d_min; - end - jacobSub = constraint.computeDoseConstraintJacobian(d_i); - - case 'VWWC_INV' %inverse voxel-wise conformitiy - takes maximum dose in TARGET and minimum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D inverted VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - - d_max = max(d_Scen,[],2); - d_min = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_min; - elseif isequal(cst{i,3},'TARGET') - d_i = d_max; - end - jacobSub = constraint.computeDoseConstraintJacobian(d_i); - - otherwise - matRad_cfg.dispError('Robustness setting %s not yet supported!',constraint.robustness); + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end + + d_Scen = d_tmp(cst{curConIdx,4}{contourIx},:); + + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{curConIdx,3},'OAR') + d_i = d_max; + elseif isequal(cst{curConIdx,3},'TARGET') + d_i = d_min; + end + jacobSub = constraint.computeDoseConstraintJacobian(d_i); + + case 'VWWC_INV' %inverse voxel-wise conformitiy - takes maximum dose in TARGET and minimum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D inverted VWWC optimization is currently not supported'); end - nConst = size(jacobSub,2); + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end + + d_Scen = d_tmp(cst{curConIdx,4}{contourIx},:); - %Iterate through columns of the sub-jacobian - if isa(optiProb.BP,'matRad_DoseProjection') && ~isempty(jacobSub) || isa(optiProb.BP,'matRad_ConstantRBEProjection') - - startIx = size(DoseProjection{1},2) + 1; - %First append the Projection matrix with sparse zeros - DoseProjection{1} = [DoseProjection{1},sparse(dij.doseGrid.numOfVoxels,nConst)]; - - %Now directly write the jacobian in there - DoseProjection{1}(cst{i,4}{1},startIx:end) = jacobSub; - - elseif isa(optiProb.BP,'matRad_EffectProjection') && ~isempty(jacobSub) - - if isa(optiProb.BP,'matRad_VariableRBEProjection') - scaledEffect = (dij.gamma(cst{i,4}{1}) + d_i); - jacobSub = jacobSub./(2*dij.bx(cst{i,4}{1}) .* scaledEffect); - end - - startIx = size(mAlphaDoseProjection{1},2) + 1; - - %First append the alphaDose matrix with sparse - %zeros then insert - mAlphaDoseProjection{1} = [mAlphaDoseProjection{1},sparse(dij.doseGrid.numOfVoxels,nConst)]; - mAlphaDoseProjection{1}(cst{i,4}{1},startIx:end) = jacobSub; - - %The betadose has a different structure due to the - %quadratic transformation, but in principle the - %same as above - mSqrtBetaDoseProjection{1} = [mSqrtBetaDoseProjection{1}, sparse(repmat(cst{i,4}{1},nConst,1),repmat(1:numel(cst{i,4}{1}),1,nConst),2*reshape(jacobSub',[],1),dij.doseGrid.numOfVoxels,nConst*numel(cst{i,4}{1}))]; - - if isempty(constraintID) - newID = 1; - else - newID = constraintID(end)+1; - end - - voxelID = [voxelID;repmat(cst{i,4}{1},nConst,1)]; %Keep track of voxels for organizing the sqrt(beta)Dose projection later - constraintID = [constraintID, ... - reshape(ones(numel(cst{i,4}{1}),1)*[newID:newID+nConst-1],[1 nConst*numel(cst{i,4}{1})])]; %Keep track of constraints for organizing the sqrt(beta)Dose projection later + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{curConIdx,3},'OAR') + d_i = d_min; + elseif isequal(cst{curConIdx,3},'TARGET') + d_i = d_max; end + jacobSub = constraint.computeDoseConstraintJacobian(d_i); - + otherwise + matRad_cfg.dispError('Robustness setting %s not yet supported!',constraint.robustness); + end + + nConst = size(jacobSub,2); + + %Iterate through columns of the sub-jacobian + if isa(optiProb.BP,'matRad_DoseProjection') && ~isempty(jacobSub) || isa(optiProb.BP,'matRad_ConstantRBEProjection') + + startIx = size(DoseProjection{1},2) + 1; + %First append the Projection matrix with sparse zeros + DoseProjection{1} = [DoseProjection{1},sparse(dij.doseGrid.numOfVoxels,nConst)]; +% DoseProjection{1} = [DoseProjection{1},zeros(dij.doseGrid.numOfVoxels,nConst)]; + + %Now directly write the jacobian in there + DoseProjection{1}(cst{curConIdx,4}{1},startIx:end) = jacobSub; + + + elseif isa(optiProb.BP,'matRad_EffectProjection') && ~isempty(jacobSub) + + if isa(optiProb.BP,'matRad_VariableRBEProjection') + scaledEffect = (dij.gamma(cst{curConIdx,4}{1}) + d_i); + jacobSub = jacobSub./(2*dij.bx(cst{curConIdx,4}{1}) .* scaledEffect); + end + + startIx = size(mAlphaDoseProjection{1},2) + 1; + + %First append the alphaDose matrix with sparse + %zeros then insert + mAlphaDoseProjection{1} = [mAlphaDoseProjection{1},sparse(dij.doseGrid.numOfVoxels,nConst)]; + mAlphaDoseProjection{1}(cst{curConIdx,4}{1},startIx:end) = jacobSub; + + %The betadose has a different structure due to the + %quadratic transformation, but in principle the + %same as above + mSqrtBetaDoseProjection{1} = [mSqrtBetaDoseProjection{1}, sparse(repmat(cst{curConIdx,4}{1},nConst,1),repmat(1:numel(cst{curConIdx,4}{1}),1,nConst),2*reshape(jacobSub',[],1),dij.doseGrid.numOfVoxels,nConst*numel(cst{curConIdx,4}{1}))]; + + if isempty(constraintID) + newID = 1; + else + newID = constraintID(end)+1; end + voxelID = [voxelID;repmat(cst{curConIdx,4}{1},nConst,1)]; %Keep track of voxels for organizing the sqrt(beta)Dose projection later + constraintID = [constraintID, ... + reshape(ones(numel(cst{curConIdx,4}{1}),1)*[newID:newID+nConst-1],[1 nConst*numel(cst{curConIdx,4}{1})])]; %Keep track of constraints for organizing the sqrt(beta)Dose projection later end - - end - + + end - + + scenario = 1; % enter if statement also for protons using a constant RBE if isa(optiProb.BP,'matRad_DoseProjection') - if ~isempty(DoseProjection{scenario}) - jacob = DoseProjection{scenario}' * dij.physicalDose{scenario}; + jacob = transpose(DoseProjection{scenario}) * dij.physicalDose{scenario}; %for some reason here faster than shorthand notation + end elseif isa(optiProb.BP,'matRad_ConstantRBEProjection') if ~isempty(DoseProjection{scenario}) - jacob = DoseProjection{scenario}' * dij.RBE * dij.physicalDose{scenario}; + %jacob = DoseProjection{scenario}' * dij.RBE * dij.physicalDose{scenario}; + jacob = dij.RBE * (DoseProjection{scenario}' * dij.physicalDose{scenario}); end elseif isa(optiProb.BP,'matRad_EffectProjection') @@ -220,5 +204,4 @@ mSqrtBetaDoseProjection{scenario}' * dij.mSqrtBetaDose{scenario}; end -end - +end \ No newline at end of file diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m index d376ac58e..a07248033 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_getConstraintBounds.m @@ -1,6 +1,5 @@ function [cl,cu] = matRad_getConstraintBounds(optiProb,cst) % matRad IPOPT get constraint bounds wrapper function -% % call % [cl,cu] = matRad_getConstraintBounds(optiProb,cst) % @@ -8,7 +7,8 @@ % cst: matRad cst struct % % output -% cl: lower bounds on constraints +% cl: lower bounds on constfunction [cl,cu] = matRad_getConstraintBounds(optiProb,cst) +% matRad IPOPT get constraint bounds wrapper functioraints % cu: lower bounds on constraints % % References @@ -28,61 +28,20 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -BPtype = class(optiProb.BP); -isEffectBP = strcmp(BPtype,'matRad_EffectProjection'); - % Initialize bounds cl = []; cu = []; % compute objective function for every VOI. -for i = 1:size(cst,1) - - % Only take OAR or target VOI. - if ~any(cellfun(@isempty,cst{i,4})) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) - - % loop over the number of constraints for the current VOI - for j = 1:numel(cst{i,6}) - - optiFunc = cst{i,6}{j}; - - % only perform computations for constraints -%{ - if ~isempty(strfind(cst{i,6}(j).type,'constraint')) - - if isequal(options.quantityOpt,'effect') - param = cst{i,5}.alphaX .* cst{i,6}(j).dose + cst{i,5}.betaX .* cst{i,6}(j).dose.^2; - else - param = cst{i,6}(j).dose; - end - - if strcmp(cst{i,6}(j).robustness,'none') || strcmp(cst{i,6}(j).robustness,'probabilistic') || strcmp(cst{i,6}(j).robustness,'VWWC') ||... - strcmp(cst{i,6}(j).robustness,'COWC') || strcmp(cst{i,6}(j).robustness,'VWWC_CONF')|| strcmp(cst{i,6}(j).robustness,'OWC') - - - [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 - - doses = optiFunc.getDoseParameters(); - effect = cst{i,5}.alphaX*doses + cst{i,5}.betaX*doses.^2; - - optiFunc = optiFunc.setDoseParameters(effect); - 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 +for i = 1:size(optiProb.constrIdx,1) + obj = optiProb.constraints{i}; + curConIdx = optiProb.constrIdx(i,1); - end % if structure not empty and target or oar + cl = [cl;obj.lowerBounds(numel(cst{curConIdx,4}{1}))]; + cu = [cu;obj.upperBounds(numel(cst{curConIdx,4}{1}))]; + +end -end % over all structures - + diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_getJacobianStructure.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_getJacobianStructure.m index 131621544..43f6745ae 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_getJacobianStructure.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_getJacobianStructure.m @@ -19,7 +19,7 @@ % References % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Copyright 2016 the matRad development team. % @@ -31,26 +31,16 @@ % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Initializes constraints +% Initializes constraints jacobStruct = sparse([]); - % compute objective function for every VOI. -for i = 1:size(cst,1) - % Only take OAR or target VOI. - if ~any(cellfun(@isempty,cst{i,4})) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) - % loop over the number of constraints for the current VOI - for j = 1:numel(cst{i,6}) - - obj = cst{i,6}{j}; - - % only perform computations for constraints - if isa(obj,'DoseConstraints.matRad_DoseConstraint') - - % get the jacobian structure depending on dose - jacobDoseStruct = obj.getDoseConstraintJacobianStructure(numel(cst{i,4}{1})); - nRows = size(jacobDoseStruct,2); - jacobStruct = [jacobStruct; repmat(spones(mean(dij.physicalDose{1}(cst{i,4}{1},:),1)),nRows,1)]; - - end - end - end - end \ No newline at end of file +% compute objective function for every VOI. +for i = 1:size(optiProb.constrIdx,1) + obj = optiProb.constraints{i}; + curConIdx = optiProb.constrIdx(i,1); + + % get the jacobian structure depending on dose + jacobDoseStruct = obj.getDoseConstraintJacobianStructure(numel(cst{curConIdx,4}{1})); + nRows = size(jacobDoseStruct,2); + jacobStruct = [jacobStruct; repmat(spones(mean(dij.physicalDose{1}(cst{curConIdx,4}{1},:))),nRows,1)]; + +end \ No newline at end of file diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m index 17d4abab3..b6c265361 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunction.m @@ -1,222 +1,93 @@ function f = matRad_objectiveFunction(optiProb,w,dij,cst) -% matRad IPOPT objective function wrapper -% -% call -% f = matRad_objectiveFuncWrapper(optiProb,w,dij,cst) -% -% input -% optiProb: matRad optimization problem -% w: beamlet/ pencil beam weight vector -% dij: matRad dose influence struct -% cst: matRad cst struct -% scenario: index of dij scenario to consider (optional: default 1) -% -% output -% f: objective function value -% -% References -% - -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Copyright 2016 the matRad development team. -% -% This file is part of the matRad project. It is subject to the license -% terms in the LICENSE file found in the top-level directory of this -% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part -% of the matRad project, including this file, may be copied, modified, -% propagated, or distributed except according to the terms contained in the -% LICENSE file. -% -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -matRad_cfg = MatRad_Config.instance(); - -% get current dose / effect / RBExDose vector -optiProb.BP.compute(dij,w); -d = optiProb.BP.GetResult(); - -% get probabilistic quantities (nearly no overhead if empty) -[dExp,dOmega] = optiProb.BP.GetResultProb(); - -% get the used scenarios -useScen = optiProb.BP.scenarios; -scenProb = optiProb.BP.scenarioProb; -useNominalCtScen = optiProb.BP.nominalCtScenarios; - -% retrieve matching 4D scenarios -fullScen = cell(ndims(d),1); -[fullScen{:}] = ind2sub(size(d),useScen); -contourScen = fullScen{1}; - -% initialize f -f = 0; - -% required for COWC opt -f_COWC = zeros(numel(useScen),1); + % matRad IPOPT objective function wrapper + % + % call + % f = matRad_objectiveFuncWrapper(optiProb,w,dij,cst) + % + % input + % optiProb: matRad optimization problem + % w: beamlet/ pencil beam weight vector + % dij: matRad dose influence struct + % cst: matRad cst struct + % scenario: index of dij scenario to consider (optional: default 1) + % + % output + % f: objective function value + % + % References + % - + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2016 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + fIndv = matRad_objectiveFunctions(optiProb,w,dij,cst); + %fIndv should be an mxn matrix where m is the number of objectives and n should be the number of scenarios used + + % normalization is so far only used for pareto optimization + switch optiProb.normalizationScheme.scheme + case 'UL' % f = (f-L)/(U-L) + fIndv = optiProb.normalizeObjectives(fIndv')';%ISSUES IF COWC is used - + end -% compute objective function for every VOI. -for i = 1:size(cst,1) - - % Only take OAR or target VOI. - if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) + %COWC calculation + + %get indices of objectives and their respective penalties + COWCIdxs = []; + Idxs = []; + penalties = []; + COWCpenalties = []; + f = 0; + + for i = 1:numel(optiProb.objectives) + if strcmp(optiProb.objectives{i}.robustness,'COWC') + COWCIdxs = [COWCIdxs,i]; + COWCpenalties = [COWCpenalties,optiProb.objectives{i}.penalty]; + else + Idxs = [Idxs,i]; + penalties = [penalties,optiProb.objectives{i}.penalty]; + end + end + if isempty(COWCIdxs) + f_COWC = [0]; + else + f_COWCs = fIndv(COWCIdxs,:); + f_COWC = COWCpenalties*f_COWCs; + end%sum over scenarios respect penalties - % loop over the number of constraints for the current VOI - for j = 1:numel(cst{i,6}) - - objective = cst{i,6}{j}; - - % only perform gradient computations for objectiveectives - if isa(objective,'DoseObjectives.matRad_DoseObjective') - - % rescale dose parameters to biological optimization quantity if required - objective = optiProb.BP.setBiologicalDosePrescriptions(objective,cst{i,5}.alphaX,cst{i,5}.betaX); - - % retrieve the robustness type - robustness = objective.robustness; - - switch robustness - case 'none' % if conventional opt: just sum objectives of nominal dose - for ixScen = useNominalCtScen - d_i = d{ixScen}(cst{i,4}{useScen(ixScen)}); - f = f + objective.penalty * objective.computeDoseObjectiveFunction(d_i); - end - - case 'STOCH' % if prob opt: sum up expectation value of objectives - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - d_i = d{ixScen}(cst{i,4}{ixContour}); - - f = f + scenProb(s) * objective.penalty*objective.computeDoseObjectiveFunction(d_i); - - end - - case 'PROB' % if prob opt: sum up expectation value of objectives - - d_i = dExp{1}(cst{i,4}{1}); - - f = f + objective.penalty*objective.computeDoseObjectiveFunction(d_i); - - p = objective.penalty/numel(cst{i,4}{1}); - - % only one variance term per VOI - if j == 1 - f = f + p * w' * dOmega{i,1}; - end - - case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - - d_max = max(d_Scen,[],2); - d_min = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_max; - elseif isequal(cst{i,3},'TARGET') - d_i = d_min; - end - - f = f + objective.penalty*objective.computeDoseObjectiveFunction(d_i); - - case 'VWWC_INV' %inverse voxel-wise conformitiy - consider the maximum and minimum dose in the target and optimize the dose conformity - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispWarning('4D inverted VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - d_max = max(d_Scen,[],2); - d_min = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_min; - elseif isequal(cst{i,3},'TARGET') - d_i = d_max; - end - - f = f + objective.penalty*objective.computeDoseObjectiveFunction(d_i); - - case 'COWC' % composite worst case consideres ovarall the worst objective function value - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - d_i = d{ixScen}(cst{i,4}{ixContour}); - - f_COWC(s) = f_COWC(s) + objective.penalty*objective.computeDoseObjectiveFunction(d_i); - end - - case 'OWC' % objective-wise worst case considers the worst individual objective function value - - f_OWC = zeros(numel(useScen),1); - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - d_i = d{ixScen}(cst{i,4}{ixContour}); - f_OWC(s) = objective.penalty*objective.computeDoseObjectiveFunction(d_i); - end - - % compute the maximum objective function value - switch optiProb.useMaxApprox - case 'logsumexp' - fMax = optiProb.logSumExp(f_OWC); - case 'pnorm' - fMax = optiProb.pNorm(f_OWC,numel(useScen)); - case 'none' - fMax = max(f_OWC); - case 'otherwise' - matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); - fMax = max(f_OWC); - end - f = f + fMax; - - otherwise - matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); - - end %robustness type - end % objective check - end %objective loop - end %empty check -end %cst structure loop - - -%Handling the maximum of the composite worst case part -fMax = max(f_COWC(:)); -if fMax > 0 - switch optiProb.useMaxApprox - case 'logsumexp' - fMax = optiProb.logSumExp(f_COWC); - case 'pnorm' - fMax = optiProb.pNorm(f_COWC,numel(useScen)); - case 'none' - fMax = max(f_COWC); - case 'otherwise' - matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); - fMax = max(f_COWC); + fMax = max(f_COWC(:)); + if fMax > 0 + switch optiProb.useMaxApprox + case 'logsumexp' + fMax = optiProb.logSumExp(f_COWC); + case 'pnorm' + fMax = optiProb.pNorm(f_COWC,numel(useScen)); + case 'none' + fMax = max(f_COWC); + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + fMax = max(f_COWC); + end + end + + f = f + fMax; + + %now sum up over other robustness scenarios + if isempty(Idxs) + f_wo_COWC = [0]; + else + f_wo_COWCs = fIndv(Idxs,:); + f_wo_COWC = sum(penalties*f_wo_COWCs,2); %sum over scenarios? end -end -%Sum up max of composite worst case part -f = f + fMax; + f = f + f_wo_COWC; +end \ No newline at end of file diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunctions.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunctions.m new file mode 100644 index 000000000..5bd7382a2 --- /dev/null +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveFunctions.m @@ -0,0 +1,203 @@ +function fIndv = matRad_objectiveFunctions(optiProb,w,dij,cst) +% matRad IPOPT objective function wrapper +% +% call +% fIndv = matRad_objectiveFunctions(optiProb,w,dij,cst) +% +% input +% optiProb: matRad optimization problem +% w: beamlet/ pencil beam weight vector +% dij: matRad dose influence struct +% cst: matRad cst struct +% scenario: index of dij scenario to consider (optional: default 1) +% +% output +% fIndv: m x n matrix that stores all objective function values +% m: number of objectives +% n: number of scenarios +% there are some inconsistencies when using robust optimization +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2016 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + matRad_cfg = MatRad_Config.instance(); + + % get current dose / effect / RBExDose vector + optiProb.BP.compute(dij,w); + d = optiProb.BP.GetResult(); + + % get probabilistic quantities (nearly no overhead if empty) + [dExp,dOmega] = optiProb.BP.GetResultProb(); + + % get the used scenarios + useScen = optiProb.BP.scenarios; + scenProb = optiProb.BP.scenarioProb; + useNominalCtScen = optiProb.BP.nominalCtScenarios; + + % retrieve matching 4D scenarios + fullScen = cell(ndims(d),1); + [fullScen{:}] = ind2sub(size(d),useScen); + contourScen = fullScen{1}; + + % initialize matrix with objective function values + fIndv = [zeros(numel(optiProb.objectives),numel(useScen))]; + + + %individual objective functions + for i = 1:size(optiProb.objIdx,1) %loop over objectives + + objective = optiProb.objectives{i}; + curObjIdx = optiProb.objIdx(i,1); + + %calculation differs based on robustness + robustness = objective.robustness; + + + switch robustness + case 'none' % if conventional opt: just sum objectives of nominal dose + + for ixScen = useNominalCtScen + d_i = d{ixScen}(cst{curObjIdx,4}{useScen(1)}); + fInd = objective.computeDoseObjectiveFunction(d_i); + fIndv(i,ixScen) = fInd; + end + + case 'STOCH' % if prob opt: sum up expectation value of objectives + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + + fInd = scenProb(s) * objective.computeDoseObjectiveFunction(d_i); + fIndv(i,s) = fInd; + + end + + case 'PROB' % if prob opt: sum up expectation value of objectives TODO: CHECK FOR VALUE TO APPEND + + %NOT UPDATED YET! + d_i = dExp{1}(cst{curObjIdx,4}{1}); + + fInd = objective.computeDoseObjectiveFunction(d_i); + fIndv = [fIndv fInd] + + p = 1/numel(cst{curObjIdx,4}{1}); + + % only one variance term per VOI + if j == 1 + %need to add + fInd = p * w' * dOmega{i,1}; + fIndv(end) = fIndv(end) + fInd + end + + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end + + d_Scen = d_tmp(cst{curObjIdx,4}{contourIx},:); + + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{curObjIdx,3},'OAR') + d_i = d_max; + elseif isequal(cst{curObjIdx,3},'TARGET') + d_i = d_min; + end + + fInd = objective.computeDoseObjectiveFunction(d_i); + fIndv(i,1) = fInd; %no different scenarios? + + case 'VWWC_INV' %inverse voxel-wise conformitiy - consider the maximum and minimum dose in the target and optimize the dose conformity + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispWarning('4D inverted VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end + + d_Scen = d_tmp(cst{curObjIdx,4}{contourIx},:); + d_max = max(d_Scen,[],2); + d_min = min(d_Scen,[],2); + + if isequal(cst{curObjIdx,3},'OAR') + d_i = d_min; + elseif isequal(cst{curObjIdx,3},'TARGET') + d_i = d_max; + end + + fInd = objective.computeDoseObjectiveFunction(d_i); + fIndv(i,1) = fInd; + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + fIndv(i,s) = objective.computeDoseObjectiveFunction(d_i); + end + + case 'OWC' % objective-wise worst case considers the worst individual objective function value + + f_OWC = zeros(numel(useScen),1); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + f_OWC(s) = objective.computeDoseObjectiveFunction(d_i); + end + + % compute the maximum objective function value + switch optiProb.useMaxApprox + case 'logsumexp' + fMax = optiProb.logSumExp(f_OWC); + case 'pnorm' + fMax = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + fMax = max(f_OWC); + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + fMax = max(f_OWC); + end + fIndv(i,1) = fMax; %does this make sense? + + otherwise + matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); + + end %robustness type + end +end \ No newline at end of file diff --git a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m index 1d92ff56c..2a2c7de4b 100644 --- a/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m +++ b/matRad/optimization/@matRad_OptimizationProblem/matRad_objectiveGradient.m @@ -33,7 +33,6 @@ % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - matRad_cfg = MatRad_Config.instance(); % get current dose / effect / RBExDose vector @@ -56,6 +55,7 @@ doseGradient = cell(size(dij.physicalDose)); doseGradient(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + %For probabilistic optimization vOmega = 0; @@ -63,213 +63,199 @@ f_COWC = zeros(size(dij.physicalDose)); % compute objective function for every VOI. -for i = 1:size(cst,1) - - % Only take OAR or target VOI. - if ~isempty(cst{i,4}{1}) && ( isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET') ) +for i = 1:size(optiProb.objIdx,1) + objective = optiProb.objectives{i}; + curObjIdx = optiProb.objIdx(i,1); + + % retrieve the robustness type + robustness = objective.robustness; + + + switch robustness + case 'none' % if conventional opt: just sum objectives of nominal dose + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + %add to dose gradient + doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) = doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) + objective.penalty * optiProb.normalizeGradient(objective.computeDoseObjectiveGradient(d_i),i); + end + case 'STOCH' % perform stochastic optimization with weighted / random scenarios + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + + doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) = doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) + ... + (objective.penalty*objective.computeDoseObjectiveGradient(d_i) * scenProb(s)); + + end + + case 'PROB' % use the expectation value and the integral variance influence matrix + %First check the speficic cache for probabilistic + if ~exist('doseGradientExp','var') + doseGradientExp{1} = zeros(dij.doseGrid.numOfVoxels,1); + end + + d_i = dExp{1}(cst{curObjIdx,4}{1}); + + doseGradientExp{1}(cst{curObjIdx,4}{1}) = doseGradientExp{1}(cst{curObjIdx,4}{1}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); + + p = objective.penalty/numel(cst{curObjIdx,4}{1}); + + vOmega = vOmega + p * dOmega{i,1}; - % loop over the number of constraints and objectives for the current VOI - for j = 1:numel(cst{i,6}) + case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end + + d_Scen = d_tmp(cst{curObjIdx,4}{contourIx},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{curObjIdx,3},'OAR') + d_i = d_max; + elseif isequal(cst{curObjIdx,3},'TARGET') + d_i = d_min; + end - %Get current optimization function - objective = cst{i,6}{j}; + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradient.',numel(isnan(d_i))); + end - % only perform gradient computations for objectives - if isa(objective,'DoseObjectives.matRad_DoseObjective') + deltaTmp = objective.penalty*objective.computeDoseObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); - % retrieve the robustness type - robustness = objective.robustness; + if isequal(cst{curObjIdx,3},'OAR') + currWcIx = double(max_ix == s); + elseif isequal(cst{curObjIdx,3},'TARGET') + currWcIx = double(min_ix == s); + end - % rescale dose parameters to biological optimization quantity if required - objective = optiProb.BP.setBiologicalDosePrescriptions(objective,cst{i,5}.alphaX,cst{i,5}.betaX); + doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) = doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'VWWC_INV' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR + contourIx = unique(contourScen); + if ~isscalar(contourIx) + % voxels need to be tracked through the 4D CT, + % not yet implemented + matRad_cfg.dispError('4D VWWC optimization is currently not supported'); + end + + % prepare min/max dose vector for voxel-wise worst case + if ~exist('d_tmp','var') + d_tmp = [d{useScen}]; + end + + d_Scen = d_tmp(cst{curObjIdx,4}{1},:); + [d_max,max_ix] = max(d_Scen,[],2); + [d_min,min_ix] = min(d_Scen,[],2); + + if isequal(cst{curObjIdx,3},'OAR') + d_i = d_min; + elseif isequal(cst{curObjIdx,3},'TARGET') + d_i = d_max; + end + + if any(isnan(d_i)) + matRad_cfg.dispWarning('%d NaN values in gradFuncWrapper.',numel(isnan(d_i))); + end + + deltaTmp = objective.penalty*objective.computeDoseObjectiveGradient(d_i); + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + if isequal(cst{curObjIdx,3},'OAR') + currWcIx = double(min_ix == s); + elseif isequal(cst{curObjIdx,3},'TARGET') + currWcIx = double(max_ix == s); + end + + doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) = doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) + deltaTmp.*currWcIx; + end + + case 'COWC' % composite worst case consideres ovarall the worst objective function value + %First check the speficic cache for COWC + if ~exist('delta_COWC','var') + delta_COWC = cell(size(doseGradient)); + delta_COWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + + f_COWC(ixScen) = f_COWC(ixScen) + objective.penalty*objective.computeDoseObjectiveFunction(d_i); + delta_COWC{ixScen}(cst{curObjIdx,4}{ixContour}) = delta_COWC{ixScen}(cst{curObjIdx,4}{ixContour}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); + end + + case 'OWC' % objective-wise worst case consideres the worst individual objective function value + %First check the speficic cache for COWC + f_OWC = zeros(size(doseGradient)); + + if ~exist('delta_OWC','var') + delta_OWC = cell(size(doseGradient)); + delta_OWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + + d_i = d{ixScen}(cst{curObjIdx,4}{ixContour}); + + f_OWC(ixScen) = objective.penalty*objective.computeDoseObjectiveFunction(d_i); + + delta_OWC{ixScen}(cst{curObjIdx,4}{ixContour}) = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - switch robustness - case 'none' % if conventional opt: just sum objectiveectives of nominal dose - for s = useNominalCtScen - ixScen = useScen(s); - ixContour = contourScen(s); - d_i = d{ixScen}(cst{i,4}{ixContour}); - %add to dose gradient - doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); - end - case 'STOCH' % perform stochastic optimization with weighted / random scenarios - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - d_i = d{ixScen}(cst{i,4}{ixContour}); - - doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + ... - (objective.penalty*objective.computeDoseObjectiveGradient(d_i) * scenProb(s)); - - end - - case 'PROB' % use the expectation value and the integral variance influence matrix - %First check the speficic cache for probabilistic - if ~exist('doseGradientExp','var') - doseGradientExp{1} = zeros(dij.doseGrid.numOfVoxels,1); - end - - d_i = dExp{1}(cst{i,4}{1}); - - doseGradientExp{1}(cst{i,4}{1}) = doseGradientExp{1}(cst{i,4}{1}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); - - p = objective.penalty/numel(cst{i,4}{1}); - - vOmega = vOmega + p * dOmega{i,1}; - - case 'VWWC' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector for voxel-wise worst case - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{contourIx},:); - [d_max,max_ix] = max(d_Scen,[],2); - [d_min,min_ix] = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_max; - elseif isequal(cst{i,3},'TARGET') - d_i = d_min; - end - - if any(isnan(d_i)) - matRad_cfg.dispWarning('%d NaN values in gradient.',numel(isnan(d_i))); - end - - deltaTmp = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - if isequal(cst{i,3},'OAR') - currWcIx = double(max_ix == s); - elseif isequal(cst{i,3},'TARGET') - currWcIx = double(min_ix == s); - end - - doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; - end - - case 'VWWC_INV' % voxel-wise worst case - takes minimum dose in TARGET and maximum in OAR - contourIx = unique(contourScen); - if ~isscalar(contourIx) - % voxels need to be tracked through the 4D CT, - % not yet implemented - matRad_cfg.dispError('4D VWWC optimization is currently not supported'); - end - - % prepare min/max dose vector for voxel-wise worst case - if ~exist('d_tmp','var') - d_tmp = [d{useScen}]; - end - - d_Scen = d_tmp(cst{i,4}{1},:); - [d_max,max_ix] = max(d_Scen,[],2); - [d_min,min_ix] = min(d_Scen,[],2); - - if isequal(cst{i,3},'OAR') - d_i = d_min; - elseif isequal(cst{i,3},'TARGET') - d_i = d_max; - end - - if any(isnan(d_i)) - matRad_cfg.dispWarning('%d NaN values in gradFuncWrapper.',numel(isnan(d_i))); - end - - deltaTmp = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - if isequal(cst{i,3},'OAR') - currWcIx = double(min_ix == s); - elseif isequal(cst{i,3},'TARGET') - currWcIx = double(max_ix == s); - end - - doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + deltaTmp.*currWcIx; - end - - case 'COWC' % composite worst case consideres ovarall the worst objective function value - %First check the speficic cache for COWC - if ~exist('delta_COWC','var') - delta_COWC = cell(size(doseGradient)); - delta_COWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; - end - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - d_i = d{ixScen}(cst{i,4}{ixContour}); - - f_COWC(ixScen) = f_COWC(ixScen) + objective.penalty*objective.computeDoseObjectiveFunction(d_i); - delta_COWC{ixScen}(cst{i,4}{ixContour}) = delta_COWC{ixScen}(cst{i,4}{ixContour}) + objective.penalty*objective.computeDoseObjectiveGradient(d_i); - end - - case 'OWC' % objective-wise worst case consideres the worst individual objective function value - %First check the speficic cache for COWC - f_OWC = zeros(size(doseGradient)); - - if ~exist('delta_OWC','var') - delta_OWC = cell(size(doseGradient)); - delta_OWC(useScen) = {zeros(dij.doseGrid.numOfVoxels,1)}; - end - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - - d_i = d{ixScen}(cst{i,4}{ixContour}); - - f_OWC(ixScen) = objective.penalty*objective.computeDoseObjectiveFunction(d_i); - - delta_OWC{ixScen}(cst{i,4}{ixContour}) = objective.penalty*objective.computeDoseObjectiveGradient(d_i); - - end - - switch optiProb.useMaxApprox - case 'logsumexp' - [~,fGrad] = optiProb.logSumExp(f_OWC); - case 'pnorm' - [~,fGrad] = optiProb.pNorm(f_OWC,numel(useScen)); - case 'none' - [~,ix] = max(f_OWC(:)); - fGrad = zeros(size(f_OWC)); - fGrad(ix) = 1; - case 'otherwise' - matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); - [~,ix] = max(f_OWC(:)); - fGrad = zeros(size(f_OWC)); - fGrad(ix) = 1; - end - - for s = 1:numel(useScen) - ixScen = useScen(s); - ixContour = contourScen(s); - if fGrad(ixScen ) ~= 0 - doseGradient{ixScen}(cst{i,4}{ixContour}) = doseGradient{ixScen}(cst{i,4}{ixContour}) + fGrad(ixScen)*delta_OWC{ixScen}(cst{i,4}{ixContour}); - end - end - - otherwise - matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); - - end %robustness type - end % objective check - end %objective loop + end + + switch optiProb.useMaxApprox + case 'logsumexp' + [~,fGrad] = optiProb.logSumExp(f_OWC); + case 'pnorm' + [~,fGrad] = optiProb.pNorm(f_OWC,numel(useScen)); + case 'none' + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + case 'otherwise' + matRad_cfg.dispWarning('Unknown maximum approximation desired. Using ''none'' instead.'); + [~,ix] = max(f_OWC(:)); + fGrad = zeros(size(f_OWC)); + fGrad(ix) = 1; + end + + for s = 1:numel(useScen) + ixScen = useScen(s); + ixContour = contourScen(s); + if fGrad(ixScen ) ~= 0 + doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) = doseGradient{ixScen}(cst{curObjIdx,4}{ixContour}) + fGrad(ixScen)*delta_OWC{ixScen}(cst{curObjIdx,4}{ixContour}); + end + end + + otherwise + matRad_cfg.dispError('Robustness setting %s not supported!',objective.robustness); + + %objective loop end %empty check end %cst structure loop @@ -304,7 +290,7 @@ g = optiProb.BP.GetGradient(); for s = 1:numel(useScen) - weightGradient = weightGradient + g{useScen(s)}; + weightGradient = weightGradient + g{useScen(s)}; end if vOmega ~= 0 @@ -313,4 +299,4 @@ %Only implemented for first scenario now weightGradient = weightGradient + gProb{1}; -end +end \ No newline at end of file diff --git a/matRad/optimization/lexicographic/matRad_2pecOptimization.m b/matRad/optimization/lexicographic/matRad_2pecOptimization.m new file mode 100644 index 000000000..902889f00 --- /dev/null +++ b/matRad/optimization/lexicographic/matRad_2pecOptimization.m @@ -0,0 +1,151 @@ +function [resultGUIs,resultGUIs2,cst1,cst2,PriorityList2]= matRad_2pecOptimization(PriorityList1,dij,cst,pln,wInit) + % Lexicographic optimization using the 2pec method as introducded by + % Breedveld et al + % + % call + % + % input + % PriorityList1: Initial structure to store all objectives + % dij: matRad dij struct + % cst: modified matRad cst struct. + % pln: matRad pln struct + % wInit: (optional) custom weights to initialize problems + % + % output + % resultGUIs: Plans created in the first phase + % resultGUIs2: Plans created in teh second phase + % cst1: Modified cst after the first phase + % cst2: Modified cst as after the second phase + % PriorityList2: Final priority list for lexicographic optimization + % + % References + % https://iopscience.iop.org/article/10.1088/0031-9155/54/23/011 + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + matRad_cfg = MatRad_Config.instance(); + + PriorityList1.adaptToFractionSize(pln.numOfFractions); + + if exist('wInit','var') + [dij,cst,pln,wInit,optiProb] = matRad_initOptimization(dij,cst,pln,wInit); + else + [dij,cst,pln,wInit,optiProb] = matRad_initOptimization(dij,cst,pln); + end + + if ~isfield(pln.propOpt,'optimizer') + pln.propOpt.optimizer = 'IPOPT'; + end + + switch pln.propOpt.optimizer + case 'IPOPT' + optimizer = matRad_OptimizerIPOPT; + case 'fmincon' + optimizer = matRad_OptimizerFmincon; + otherwise + warning(['Optimizer ''' pln.propOpt.optimizer ''' not known! Fallback to IPOPT!']); + optimizer = matRad_OptimizerIPOPT; + end + + if ~optimizer.IsAvailable() + matRad_cfg.dispError(['Optimizer ''' pln.propOpt.optimizer ''' not available!']); + end + + resultGUIs = {}; % resultGUIs for first phase + resultGUIs2 = {}; % resultGUIs for second phase + PriorityList2 = matRad_PriorityList2(); % create empty priority list + PriorityList2.ConstraintList = PriorityList1.ConstraintList; %same constraints + + %initial update for cst with all objectives and constraints for overlap + cst = PriorityList1.generateFullCst(cst); + + %generate the actual csts that get modified throughout the algorithm + cst1 = PriorityList1.generateConstraintCst(cst); + cst2 = cst1; + + matRad_cfg.propOpt.defaultAccChangeTol = 1e-6; + %add first objective(s) to cst + cst1 = PriorityList1.modifyCst(cst1); + while PriorityList1.numOfObj <= numel(PriorityList1.GoalList) + + optiProb.extractObjectivesAndConstraintsFromCst(cst1); + + + %check if objective can be skipped + [objectives,FastCalc] = PriorityList1.fastObjectiveCalc(dij,cst1,optiProb,wInit); + + if ~FastCalc % can the objective be skipped? + %If not optimize the current objective + optimizer = optimizer.optimize(wInit,optiProb,dij,cst1); + + wOpt = optimizer.wResult; + info = optimizer.resultInfo; + + resultGUI = matRad_calcCubes(wOpt,dij); + resultGUI.wUnsequenced = wOpt; + resultGUI.usedOptimizer = optimizer; + resultGUI.info = info; + resultGUI.optiProb = optiProb; + objectives = matRad_objectiveFunctions(optiProb,wOpt,dij,cst1); + + wInit = wOpt; + + else %objectives can be met + resultGUI = matRad_calcCubes(wInit,dij); + resultGUI.wUnsequenced = wInit; + resultGUI.optiProb = optiProb; + objectives = matRad_objectiveFunctions(optiProb,wInit,dij,cst1); + end + + resultGUI.objectives = objectives; %add objectives + + resultGUIs{end+1} = resultGUI; + % exchange previously optimized objectives to constraints and remove from priorityList + [cst1,cst2,PriorityList2] = PriorityList1.updateStep(cst1,cst2,PriorityList2,objectives); + + %add next objective + if PriorityList1.numOfObj <= numel(PriorityList1.GoalList) + cst1 = PriorityList1.modifyCst(cst1); + end + end + + %step2 + matRad_cfg.propOpt.defaultAccChangeTol = 1e-6; + cst2 = PriorityList2.modifyCst(cst2); %should set objective in appropriate spot -> use VOIIdx + while PriorityList2.numOfObj <= numel(PriorityList2.GoalList) + optiProb.extractObjectivesAndConstraintsFromCst(cst2); + optimizer = optimizer.optimize(wInit,optiProb,dij,cst2); + + wOpt = optimizer.wResult; + info = optimizer.resultInfo; + + resultGUI = matRad_calcCubes(wOpt,dij); + resultGUI.wUnsequenced = wOpt; + resultGUI.usedOptimizer = optimizer; + resultGUI.info = info; + resultGUI.optiProb = optiProb; + objectives = matRad_objectiveFunctions(optiProb,wOpt,dij,cst2); + resultGUI.objectives = objectives; + wInit = wOpt; + + + resultGUIs2{end+1} = resultGUI; + [cst2] = PriorityList2.updateStep(cst2,resultGUI.objectives); + + if PriorityList2.numOfObj <= numel(PriorityList2.GoalList) + cst2 = PriorityList2.modifyCst(cst2); + end + + end +end \ No newline at end of file diff --git a/matRad/optimization/lexicographic/matRad_PriorityClass.m b/matRad/optimization/lexicographic/matRad_PriorityClass.m new file mode 100644 index 000000000..cc445493b --- /dev/null +++ b/matRad/optimization/lexicographic/matRad_PriorityClass.m @@ -0,0 +1,146 @@ +classdef(Abstract) matRad_PriorityClass < handle +% matRad_PriorityClass creates the super class for PriorityLists used in +% the lexicographic approach +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2024 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + properties + Priorities; + ConstraintList; + GoalList; + slackVariable = 1.03; + numOfObj; + end + + + methods + function obj = matRad_PriorityClass() + %add optional slack Variable adjustment + obj.Priorities = []; % list storing all priorities (e.g [1,2,2,3]) + obj.GoalList = {}; + obj.ConstraintList = {}; + obj.numOfObj = 1; + + end + + function addObjective(obj,Priority,objective,goal,cstIdx) + % function to add an objective to a priorityList + % + % input + % Priority: priority of the objective in the priority list + % objective: corresponding objective (either struct or class based) + % goal: goal value associated with the objetive + % cstIdx: Index of the VOI in the cst that this objective + % belongs to + + obj.Priorities = [obj.Priorities,Priority]; %add priority + + if ~isa(objective,'matRad_DoseOptimizationFunction') + objective = matRad_DoseOptimizationFunction.createInstanceFromStruct(objective); + end + + obj.GoalList{end+1} = matRad_PriorityListObjective(objective,goal,cstIdx); + %sort by priority + [obj.Priorities,I] = sort(obj.Priorities); + obj.GoalList = obj.GoalList(I); + end + + function addConstraint(obj,constraint,cstIdx) + % function to add an constraint to a priorityList + % + % input + % constraint: Constraint to be added (either struct or class based) + % cstIdx: Index of the VOI in the cst that this objective + % belongs to + % + if ~isa(constraint,'matRad_DoseOptimizationFunction') + constraint = matRad_DoseOptimizationFunction.createInstanceFromStruct(constraint); + end + + obj.ConstraintList{end+1} = matRad_PriorityListConstraint(constraint,cstIdx); + end + + function adaptToFractionSize(obj,numOfFractions) + % function to adapt all constraints and objectives to the fraction + % size level + %adapt constraints to fraction size + for i = 1:numel(obj.ConstraintList) + constraint = obj.ConstraintList{i}.constraint; + obj.ConstraintList{i}.constraint = constraint.setDoseParameters(constraint.getDoseParameters()/numOfFractions); + end + %adapt objectives dose and goal values to fraction size + for i= 1:numel(obj.GoalList) + objective = obj.GoalList{i}.objective; + obj.GoalList{i}.objective = objective.setDoseParameters(objective.getDoseParameters()/numOfFractions); + obj.GoalList{i}.goalValue = objective.adaptGoalToFraction(obj.GoalList{i}.goalValue,numOfFractions); + end + end + %} + + function cst = generateConstraintCst(obj,cst) + %generate a bare bone cst struct for prioritized optimization containing only the + % hard constraints + cst(:,6) = cell(1); + for i = 1:numel(obj.ConstraintList) + cst{obj.ConstraintList{i}.cstIdx,6}{end+1} = obj.ConstraintList{i}.constraint; + end + end + + function [Priority,LexiObjectives] = nextPriority(obj) %get next element(s) in prioList + Priority = obj.Priorities(obj.numOfObj); %get next highest priority + LexiObjectives = {}; %cell that is filled with the objectives of next highest prio + i = obj.numOfObj; + % get all objectives for this Priority + while i <= numel(obj.Priorities) && obj.Priorities(i) == Priority + LexiObjectives{end+1} = obj.GoalList{i}; + i = i + 1; + end + end + + function printPriorityList(obj,cst) + %%% + % Function to print out the current PriorityList in table form + %%% + %create data + goals = []; + VOINames= {}; + Objectives = {}; + AchievedValues = []; + AchievedValues2 = []; + for i = 1:numel(obj.GoalList) + idx = obj.GoalList{i}.cstIdx; + goals = [goals;obj.GoalList{i}.goalValue]; + VOINames = [VOINames;cst{idx,2}]; + Objectives = [Objectives;obj.GoalList{i}.objective.name]; + AchievedValues = [AchievedValues;obj.GoalList{i}.achievedValue]; + AchievedValues2 = [AchievedValues2;obj.GoalList{i}.achievedValue2]; + end + Priority = obj.Priorities'; + T = table(Priority,VOINames,Objectives,goals,AchievedValues,AchievedValues2); + fig = uifigure(); + uitable(fig,"Data",T); + end + + + end + + methods (Abstract) + + modifyCst + updateStep + end + +end \ No newline at end of file diff --git a/matRad/optimization/lexicographic/matRad_PriorityList1.m b/matRad/optimization/lexicographic/matRad_PriorityList1.m new file mode 100644 index 000000000..f0c6306fe --- /dev/null +++ b/matRad/optimization/lexicographic/matRad_PriorityList1.m @@ -0,0 +1,96 @@ +classdef matRad_PriorityList1 < matRad_PriorityClass +% matRad_PriorityList1 implements a class based wish list for the first +% phase of the ptwo method +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2024 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + methods + + function obj = matRad_PriorityList1() + obj = obj@matRad_PriorityClass(); + end + + function [cst,Priority] = modifyCst(obj,cst) + %get the next objective(s) from the priority list + + [Priority,nextObjectives] = obj.nextPriority(); % get next Priority + for i = 1:numel(nextObjectives) + cst{nextObjectives{i}.cstIdx,6}{end+1} = nextObjectives{i}.objective; + %set number where objective is stored in cst + nextObjectives{i}.setVOIIdx(numel(cst{nextObjectives{i}.cstIdx,6})); + end + end + + function [cst1,cst2,PriorityList2] = updateStep(obj,cst1,cst2,PriorityList2,fopt) + %Change run objectives to constraints and add them to the second priority list + + [Priority,LexiObjectives] = obj.nextPriority(); + %remove from current PriorityList + obj.numOfObj = obj.numOfObj + numel(LexiObjectives); + + %change objectives in cst to constraints and add them to PriorityList2 + for i = 1:numel(LexiObjectives) + LexiObjective = LexiObjectives{i}; + LexiObjective.achievedValue = fopt(i); + bound = LexiObjective.goalValue; + if bound < fopt(i) && fopt(i) > 1e-4 %not met + bound = fopt(i)*obj.slackVariable; + else + %objectives are used in PriorityList2 + PriorityList2.GoalList{end+1} = LexiObjective; + PriorityList2.Priorities = [PriorityList2.Priorities,Priority]; + end + %update both csts + cst1{LexiObjective.cstIdx,6}{LexiObjective.getVOIIdx()} = LexiObjective.objective.turnIntoLexicographicConstraint(bound); + cst2{LexiObjective.cstIdx,6}{LexiObjective.getVOIIdx()} = LexiObjective.objective.turnIntoLexicographicConstraint(bound); + end + end + + function [objectives,skip] = fastObjectiveCalc(obj,dij,cst,optiProb,wInit) + %fast check if objectives are already met for next iteration + %and do not have to be recalculated + + [Priority,LexiObjectives] = obj.nextPriority(); + + %preassigning + bounds = zeros(numel(LexiObjectives),1); + objidxs = []; + objectives = {}; + for i = 1:numel(LexiObjectives) + bounds(i) = LexiObjectives{i}.goalValue; + objidxs = [objidxs;LexiObjectives{i}.cstIdx,LexiObjectives{i}.getVOIIdx()]; + objectives{end+1} = LexiObjectives{i}.objective; + end + + %modify optimizationProblem + optiProb.objIdx = objidxs; + optiProb.objectives = objectives; + objectives = matRad_objectiveFunctions(optiProb,wInit,dij,cst); + + skip = all(objectives < bounds); + end + + function cst = generateFullCst(obj,cst) + cst = obj.generateConstraintCst(cst); %add constraints + for i = 1:numel(obj.GoalList) + cst{obj.GoalList{i}.cstIdx,6}{end+1} = obj.GoalList{i}.objective; + end + end + + end +end \ No newline at end of file diff --git a/matRad/optimization/lexicographic/matRad_PriorityList2.m b/matRad/optimization/lexicographic/matRad_PriorityList2.m new file mode 100644 index 000000000..ceb0d5117 --- /dev/null +++ b/matRad/optimization/lexicographic/matRad_PriorityList2.m @@ -0,0 +1,53 @@ +classdef matRad_PriorityList2 < matRad_PriorityClass +% matRad_PriorityList2 implements a class based wish list for the first +% phase of the ptwo method +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + methods + + function obj = matRad_PriorityList2() + %constructor + obj = obj@matRad_PriorityClass(); + end + + function [cst,Priority]= modifyCst(obj,cst) + %get the next objective(s) from the priority list and add to cst + [Priority,nextObjectives] = obj.nextPriority(); + for i = 1:numel(nextObjectives) + %use previously allocated position + cst{nextObjectives{i}.cstIdx,6}{nextObjectives{i}.getVOIIdx} = nextObjectives{i}.objective; + end + + end + + function [cst2] = updateStep(obj,cst2,fopt) + %change the objective to constraint + [Priority,LexiObjectives] = obj.nextPriority(); + %remove objective(s) from current PriorityList + + obj.numOfObj = obj.numOfObj + numel(LexiObjectives); + + for i = 1:numel(LexiObjectives) + LexiObjective = LexiObjectives{i}; + LexiObjective.achievedValue2 = fopt(i); + bound = fopt(i)*obj.slackVariable; + cst2{LexiObjective.cstIdx,6}{LexiObjective.getVOIIdx()} = LexiObjective.objective.turnIntoLexicographicConstraint(bound); + end + end + end +end \ No newline at end of file diff --git a/matRad/optimization/lexicographic/matRad_PriorityListConstraint.m b/matRad/optimization/lexicographic/matRad_PriorityListConstraint.m new file mode 100644 index 000000000..863dc0713 --- /dev/null +++ b/matRad/optimization/lexicographic/matRad_PriorityListConstraint.m @@ -0,0 +1,33 @@ +classdef matRad_PriorityListConstraint < handle +% matRad_PriorityListConstraint implements a helper class for the priority list +% object +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + constraint; + cstIdx; % idx of the volume in the cst + end + + methods + function obj = matRad_PriorityListConstraint(constraint,cstIdx) + %constructor + obj.constraint = constraint; + obj.cstIdx = cstIdx; + end + end +end \ No newline at end of file diff --git a/matRad/optimization/lexicographic/matRad_PriorityListObjective.m b/matRad/optimization/lexicographic/matRad_PriorityListObjective.m new file mode 100644 index 000000000..90bc34b23 --- /dev/null +++ b/matRad/optimization/lexicographic/matRad_PriorityListObjective.m @@ -0,0 +1,49 @@ +classdef matRad_PriorityListObjective < handle +% matRad_PriorityListObjective implements a helper class for a priority list +% object +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2020 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + objective; + goalValue; % aspired value for the objective + cstIdx; % idx of the volume in the cst + achievedValue = -1; %not assigned yet + achievedValue2 = -1; %not assigned yet + end + + properties (Access = private) + VOIIdx; % idx of the objective for the VOI + end + + methods + function obj = matRad_PriorityListObjective(objective,goal,cstIdx) + %constructor + obj.objective = objective; + obj.goalValue = goal; + obj.cstIdx = cstIdx; + end + + function setVOIIdx(obj,idx) + obj.VOIIdx = idx; + end + + function idx = getVOIIdx(obj) + idx = obj.VOIIdx; + end + end +end \ No newline at end of file diff --git a/matRad/optimization/matRad_initOptimization.m b/matRad/optimization/matRad_initOptimization.m new file mode 100644 index 000000000..a24d3f94b --- /dev/null +++ b/matRad/optimization/matRad_initOptimization.m @@ -0,0 +1,302 @@ +function [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln,wInit) +% matRad inverse planning wrapper function +% +% call +% [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln) +% [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln,wInit) +% +% input +% dij: matRad dij struct +% cst: matRad cst struct +% pln: matRad pln struct +% wInit: (optional) custom weights to initialize problems +% +% output +% resultGUI: struct containing optimized fluence vector, dose, and (for +% biological optimization) RBE-weighted dose etc. +% optimizer: Used Optimizer Object +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2024 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +matRad_cfg = MatRad_Config.instance(); +% consider VOI priorities +cst = matRad_setOverlapPriorities(cst); + +% check & adjust objectives and constraints internally for fractionation +for i = 1:size(cst,1) + %Compatibility Layer for old objective format + if isstruct(cst{i,6}) + cst{i,6} = arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct,cst{i,6},'UniformOutput',false); + end + for j = 1:numel(cst{i,6}) + + obj = cst{i,6}{j}; + + %In case it is a default saved struct, convert to object + %Also intrinsically checks that we have a valid optimization + %objective or constraint function in the end + if ~isa(obj,'matRad_DoseOptimizationFunction') + try + obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); + catch + matRad_cfg.dispError('cst{%d,6}{%d} is not a valid Objective/constraint! Remove or Replace and try again!',i,j); + end + end + + obj = obj.setDoseParameters(obj.getDoseParameters()/pln.numOfFractions); + + cst{i,6}{j} = obj; + end +end + + + +% 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); + +% Get rid of voxels that are not interesting for the optimization problem +if ~isfield(pln.propOpt, 'clearUnusedVoxels') + pln.propOpt.clearUnusedVoxels = matRad_cfg.propOpt.defaultClearUnusedVoxels; +end + +if pln.propOpt.clearUnusedVoxels + dij = matRad_clearUnusedVoxelsFromDij(cst, dij); +end + + +% find target indices and described dose(s) for weight vector +% initialization +V = []; +doseTarget = []; +ixTarget = []; + +for i = 1:size(cst,1) + if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6}) + V = [V;cst{i,4}{1}]; + + %Iterate through objectives/constraints + fDoses = []; + for fObjCell = cst{i,6} + dParams = fObjCell{1}.getDoseParameters(); + %Don't care for Inf constraints + dParams = dParams(isfinite(dParams)); + %Add do dose list + fDoses = [fDoses dParams]; + end + + doseTarget = [doseTarget fDoses]; + ixTarget = [ixTarget i*ones(1,length(fDoses))]; + end +end + +[doseTarget,i] = max(doseTarget); +ixTarget = ixTarget(i); +wOnes = ones(dij.totalNumOfBixels,1); + +% calculate initial beam intensities wInit +matRad_cfg.dispInfo('Estimating initial weights... '); + +if exist('wInit','var') + %do nothing as wInit was passed to the function + matRad_cfg.dispInfo('chosen provided wInit!\n'); + + % Write ixDose which is needed for the optimizer + if pln.bioParam.bioOpt + for s = 1:numel(dij.bx) + dij.ixDose{s} = dij.bx{s}~=0; + end + %pre-calculations + for s = 1:numel(dij.ixDose) + dij.gamma{s} = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); + dij.gamma{s}(dij.ixDose{s}) = dij.ax{s}(dij.ixDose{s})./(2*dij.bx{s}(dij.ixDose{s})); + end + end + +elseif strcmp(pln.bioParam.model,'constRBE') && strcmp(pln.radiationMode,'protons') + % check if a constant RBE is defined - if not use 1.1 + if ~isfield(dij,'RBE') + dij.RBE = 1.1; + end + + doseTmp = dij.physicalDose{1}*wOnes; + bixelWeight = (doseTarget)/(dij.RBE * mean(doseTmp(V))); + wInit = wOnes * bixelWeight; + matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); + +elseif pln.bioParam.bioOpt + % retrieve photon LQM parameter + [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels); + checkAxBx = cellfun(@(ax1,bx1,ax2,bx2) isequal(ax1(ax1~=0),ax2(ax1~=0)) && isequal(bx1(bx1~=0),bx2(bx1~=0)),dij.ax,dij.bx,ax,bx); + if ~all(checkAxBx) + matRad_cfg.dispError('Inconsistent biological parameters in dij.ax and/or dij.bx - please recalculate dose influence matrix before optimization!\n'); + end + + for i = 1:size(cst,1) + for j = 1:size(cst{i,6},2) + % check if prescribed doses are in a valid domain + if any(cst{i,6}{j}.getDoseParameters() > 5) && isequal(cst{i,3},'TARGET') + matRad_cfg.dispError('Reference dose > 5 Gy[RBE] for target. Biological optimization outside the valid domain of the base data. Reduce dose prescription or use more fractions.\n'); + end + end + end + + for s = 1:numel(dij.bx) + dij.ixDose{s} = dij.bx{s}~=0; + end + + if isequal(pln.bioParam.quantityOpt,'effect') + + effectTarget = cst{ixTarget,5}.alphaX * doseTarget + cst{ixTarget,5}.betaX * doseTarget^2; + aTmp = dij.mAlphaDose{1}*wOnes; + bTmp = dij.mSqrtBetaDose{1} * wOnes; + p = sum(aTmp(V)) / sum(bTmp(V).^2); + q = -(effectTarget * length(V)) / sum(bTmp(V).^2); + + wInit = -(p/2) + sqrt((p^2)/4 -q) * wOnes; + + elseif isequal(pln.bioParam.quantityOpt,'RBExD') + + %pre-calculations + for s = 1:numel(dij.ixDose) + dij.gamma{s} = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); + dij.gamma{s}(dij.ixDose{s}) = dij.ax{s}(dij.ixDose{s})./(2*dij.bx{s}(dij.ixDose{s})); + end + + + % calculate current effect in target + aTmp = dij.mAlphaDose{1}*wOnes; + bTmp = dij.mSqrtBetaDose{1} * wOnes; + doseTmp = dij.physicalDose{1}*wOnes; + + CurrEffectTarget = aTmp(V) + bTmp(V).^2; + % ensure a underestimated biological effective dose + TolEstBio = 1.2; + % calculate maximal RBE in target + 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 + + matRad_cfg.dispInfo('chosen weights adapted to biological dose calculation!\n'); +else + doseTmp = dij.physicalDose{1}*wOnes; + bixelWeight = (doseTarget)/mean(doseTmp(V)); + wInit = wOnes * bixelWeight; + matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); +end + +%% calculate probabilistic quantities for probabilistic optimization if at least +% one robust objective is defined + +%Check how to use 4D data +if isfield(pln,'propOpt') && isfield(pln.propOpt,'scen4D') + scen4D = pln.propOpt.scen4D; +else + scen4D = 1; %Use only first 4D scenario for optimization +end + +%If "all" provided, use all scenarios +if isequal(scen4D,'all') + scen4D = 1:size(dij.physicalDose,1); +end + +linIxDIJ = find(~cellfun(@isempty,dij.physicalDose(scen4D,:,:)))'; + +%Only select the indexes of the nominal ct Scenarios +linIxDIJ_nominalCT = find(~cellfun(@isempty,dij.physicalDose(scen4D,1,1)))'; + +FLAG_CALC_PROB = false; +FLAG_ROB_OPT = false; + + +for i = 1:size(cst,1) + for j = 1:numel(cst{i,6}) + if strcmp(cst{i,6}{j}.robustness,'PROB') && numel(linIxDIJ) > 1 + FLAG_CALC_PROB = true; + end + if ~strcmp(cst{i,6}{j}.robustness,'none') && numel(linIxDIJ) > 1 + FLAG_ROB_OPT = true; + end + end +end + +if FLAG_CALC_PROB + [dij] = matRad_calculateProbabilisticQuantities(dij,cst,pln); +end + + +% set optimization options +if ~FLAG_ROB_OPT || FLAG_CALC_PROB % if multiple robust objectives are defined for one structure then remove FLAG_CALC_PROB from the if clause + ixForOpt = scen4D; +else + ixForOpt = linIxDIJ; +end + +switch pln.bioParam.quantityOpt + case 'effect' + backProjection = matRad_EffectProjection; + case 'RBExD' + %Capture special case of constant RBE + if strcmp(pln.bioParam.model,'constRBE') + backProjection = matRad_ConstantRBEProjection; + else + backProjection = matRad_VariableRBEProjection; + end + case 'physicalDose' + backProjection = matRad_DoseProjection; + otherwise + warning(['Did not recognize bioloigcal setting ''' pln.probOpt.bioOptimization '''!\nUsing physical dose optimization!']); + backProjection = matRad_DoseProjection; +end + +%Give scenarios used for optimization +backProjection.scenarios = ixForOpt; +backProjection.scenarioProb = pln.multScen.scenProb; +backProjection.nominalCtScenarios = linIxDIJ_nominalCT; + +optiProb = matRad_OptimizationProblem(backProjection,cst); +optiProb.quantityOpt = pln.bioParam.quantityOpt; +if isfield(pln,'propOpt') && isfield(pln.propOpt,'useLogSumExpForRobOpt') + optiProb.useLogSumExpForRobOpt = pln.propOpt.useLogSumExpForRobOpt; +end + +%Get Bounds +if ~isfield(pln.propOpt,'boundMU') + pln.propOpt.boundMU = false; +end + +if pln.propOpt.boundMU + if (isfield(dij,'minMU') || isfield(dij,'maxMU')) && ~isfield(dij,'numParticlesPerMU') + matRad_cfg.dispWarning('Requested MU bounds but number of particles per MU not set! Bounds will not be enforced and standard [0,Inf] will be used instead!'); + elseif ~isfield(dij,'minMU') && ~isfield(dij,'maxMU') + matRad_cfg.dispWarning('Requested MU bounds but machine bounds not defined in dij.minMU & dij.maxMU! Bounds will not be enforced and standard [0,Inf] will be used instead!'); + else + if isfield(dij,'minMU') + optiProb.minimumW = dij.numParticlesPerMU .* dij.minMU / 1e6; + matRad_cfg.dispInfo('Using lower MU bounds provided in dij!\n') + end + + if isfield(dij,'maxMU') + optiProb.maximumW = dij.numParticlesPerMU .* dij.maxMU / 1e6; + matRad_cfg.dispInfo('Using upper MU bounds provided in dij!\n') + end + end +else + matRad_cfg.dispInfo('Using standard MU bounds of [0,Inf]!\n') +end + diff --git a/matRad/optimization/optimizer/matRad_OptimizerIPOPT.m b/matRad/optimization/optimizer/matRad_OptimizerIPOPT.m index d1c6d6316..1a1c6ac67 100644 --- a/matRad/optimization/optimizer/matRad_OptimizerIPOPT.m +++ b/matRad/optimization/optimizer/matRad_OptimizerIPOPT.m @@ -25,13 +25,13 @@ showPlot = true; end - properties (SetAccess = protected) + properties(GetAccess = public, SetAccess = protected) + allObjectiveFunctionValues wResult resultInfo end properties (Access = private) - allObjectiveFunctionValues axesHandle plotHandle abortRequested @@ -80,7 +80,7 @@ obj.options.acceptable_constr_viol_tol = 1e-2; % (Acc3) obj.options.acceptable_dual_inf_tol = 1e10; % (Acc4) obj.options.acceptable_compl_inf_tol = 1e10; % (Acc5) - obj.options.acceptable_obj_change_tol = 1e-4; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled + obj.options.acceptable_obj_change_tol = 1e-5; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled obj.options.max_iter = matRad_cfg.defaults.propOpt.maxIter; obj.options.max_cpu_time = 7200; @@ -120,6 +120,7 @@ end function obj = optimize(obj,w0,optiProb,dij,cst) + obj.allObjectiveFunctionValues = []; matRad_cfg = MatRad_Config.instance(); % set optimization options @@ -190,7 +191,6 @@ obj.abortRequested = false; % Empty the array of stored function values - obj.allObjectiveFunctionValues = []; end function [statusmsg,statusflag] = GetStatus(obj) diff --git a/matRad/optimization/pareto/matRad_paretoOptimization.m b/matRad/optimization/pareto/matRad_paretoOptimization.m new file mode 100644 index 000000000..6156dee2f --- /dev/null +++ b/matRad/optimization/pareto/matRad_paretoOptimization.m @@ -0,0 +1,342 @@ +function [resultGUI,returnStruct] = matRad_paretoOptimization(dij,cst,pln,nIter,wInit) +% matRad inverse pareto planning wrapper function +% +% call +% [returnStruct] = matRad_paretoOptimization(dij,cst,pln,nIter) +% [returnstruct] = matRad_paretoOptimization(dij,cst,pln,nIter,wInit) +% +% input +% dij: matRad dij struct +% cst: matRad cst struct +% pln: matRad pln struct +% nIter: number of iterations to run/points to calculate after anchor point generation +% wInit: (optional) custom weights to initialize problems +% +% output +% retStruct: Structure containing the weights of the final plans and other important information +% +% References +% - https://dx.doi.org/10.2139/ssrn.1427721 +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2016 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +matRad_cfg = MatRad_Config.instance(); +if exist('wInit','var') + [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln,wInit); +else + [dij,cst,pln,wInit,optiProb,FLAG_ROB_OPT] = matRad_initOptimization(dij,cst,pln); +end + +if ~isfield(pln.propOpt,'optimizer') + pln.propOpt.optimizer = 'IPOPT'; +end + +switch pln.propOpt.optimizer + case 'IPOPT' + optimizer = matRad_OptimizerIPOPT; + case 'fmincon' + optimizer = matRad_OptimizerFmincon; + otherwise + warning(['Optimizer ''' pln.propOpt.optimizer ''' not known! Fallback to IPOPT!']); + optimizer = matRad_OptimizerIPOPT; +end + +if ~optimizer.IsAvailable() + matRad_cfg.dispError(['Optimizer ''' pln.propOpt.optimizer ''' not available!']); +end + + +% PARETO PART +if optimizer.options.acceptable_obj_change_tol > 1e-5 + warning(['Pareto Optimization requires more accurate results and therefore small objective change tolerance!']); +end +%get number of objectives at the start + +objcount = numel(optiProb.objectives); + +%% generaete Anchor Points for optimization +penGrid = matRad_generateParetoAnchorPoints(objcount); + +%initialize matrices for weight vectors and associated objective function values +weights = zeros(numel(wInit),size(penGrid,1)); +fInd = zeros(size(penGrid,1),objcount); +objectiveFunctionVals = {}; +% loop over all anchor points + +for i = 1:size(penGrid,1) + + optiProb.updatePenalties(penGrid(i,:)); + + + optimizer = optimizer.optimize(wInit,optiProb,dij,cst); + wOpt = optimizer.wResult; + info = optimizer.resultInfo; + weights(:,i) = wOpt; + + %set values for warm start + %optimizer.optionsWarmStart.use = true; + %optimizer.optionsWarmStart.zl = info.zl; + %optimizer.optionsWarmStart.zu = info.zu; + %optimizer.optionsWarmStart.lambda = info.lambda; + fInd(i,:) = matRad_objectiveFunctions(optiProb,wOpt,dij,cst)'; + objectiveFunctionVals(end + 1) = {optimizer.allObjectiveFunctionValues}; + +end + + + +%%normalize objectives + +optiProb.normalizationScheme.name = 'UL'; +optiProb.normalizationScheme.U = max(fInd,[],1); +optiProb.normalizationScheme.L = min(fInd,[],1); + +fInd = (fInd-optiProb.normalizationScheme.L)./(optiProb.normalizationScheme.U-optiProb.normalizationScheme.L); + + +%% Maybe: Grouping based on correlation of objectives +%% Generaete further points + +%initialize array storing OPS boundaries +OPSA = []; +OPSb = []; +np = size(penGrid,1); + +% first on: "Balanced point" after anchor points (slightly different calculation for normal) +% + +[a,b,firstNormal] = matRad_normalFromFacet(fInd,1:np,1); + +penGrid(np+1,:) = abs(firstNormal); +newPen = abs(firstNormal); + +%update cst values +optiProb.updatePenalties(newPen); % change for groupings! + +%calculate first new point + +optimizer = optimizer.optimize(wInit,optiProb,dij,cst); +wOpt = optimizer.wResult; +weights = [weights,wOpt]; +objectiveFunctionVals(end + 1) = {optimizer.allObjectiveFunctionValues}; +%wInit = wOpt; +fIndv = matRad_objectiveFunctions(optiProb,wOpt,dij,cst)'; +fIndv = optiProb.normalizeObjectives(fIndv); + +fInd = [fInd;fIndv]; + +%initialize some variables that might be used if points are removed later on (should only happen with too high stopping criteria) +removedfInd = []; +removedpenGrid = []; +removedweights = []; +removedOPSA = []; +removedOPSb = []; + + +OPSA = [OPSA;-penGrid(np+1,:)]; +OPSb = [OPSb;-1*fInd(np+1,:)*penGrid(np+1,:)']; +errors = []; +allErrors = {}; + +initSize = size(penGrid,1); + +%calculating for distance measure +L = min(fInd,[],1); +U = max(fInd,[],1); +eps = U - L; + +%% calculate remaining facets + +for i = 1:nIter + fprintf('Now in iteration %i',i) + %Step 1 calculate convex Hull -> Inner approximation (IPS) and gives facets + %Rennen Algorithm + fVals = fInd; + + fValsMod = matRad_generateParetoDummyPoints(fVals,U); %generate dummy points + % + [kmod,vol] = convhulln(fValsMod); + %[kred,vol] = convhulln(fVals); + + %check for relevant facets (those that contain points of the original + %fVals set) + IPSidxs = 1:size(fVals,1); + relFacetidxs = []; + for j = 1:size(kmod,1) + if any(ismember(kmod(j,:),IPSidxs)) + relFacetidxs = [relFacetidxs,j]; + end + end + facetMods = kmod(relFacetidxs,:); + facetErrors = zeros(size(facetMods,1),1); + normals = zeros(size(facetMods)); + + %Loop over facets and calculate error + for j = 1:size(relFacetidxs,2) + [facetPoints,refPoint,normal] = matRad_normalFromFacet(fValsMod,facetMods,j); + + + %check for sign of normals (might be redundant) + if all(normal<0) + continue + end + + %now check for OPS point for facet + lb = min(fVals,[],1); + ub = max(fVals,[],1); + z = linprog(normal,OPSA,OPSb,[],[],lb,ub); + + %hyperplane distance + b = refPoint*normal; + + %calculate error for each facet + + facetErrors(j) = (b-z'*normal)/(eps*normal); + normals(j,:) = normal; + + end + allErrors(end+1) = {facetErrors}; + [A,I] = sort(facetErrors,'descend'); + + %%check for next facet to run + found = false; + facetNum= 1; + accuracy = 4; + + while ~found && facetNum <= numel(I) %loop over facets + idx = I(facetNum); + norm = normals(idx,:); + + %check if facet has already been run + if ~any(ismember(round(penGrid,accuracy),round(norm,accuracy),'rows')) + + %update weights in cst + optiProb.updatePenalties(norm); + paretoOptimal = false; + calcPointsForFacet = 0; + + while ~paretoOptimal && calcPointsForFacet < 4 %more than 4 points? + if calcPointsForFacet < 3 + optimizer = optimizer.optimize(weights(:,end-calcPointsForFacet),optiProb,dij,cst); + else + optimizer = optimizer.optimize(wInit,optiProb,dij,cst); + end + + %really necessary? + if optimizer.resultInfo.iter < 15 %sometimes optimizer will terminate fast (points tend to be not optimal!) + calcPointsForFacet = calcPointsForFacet + 1; + continue + end + %check number of iterations needed for optimzation + wOpt = optimizer.wResult; + + fIndv = matRad_objectiveFunctions(optiProb,wOpt,dij,cst)'; + fIndv = optiProb.normalizeObjectives(fIndv); + %how does the newly generated point influence the pareto + %surface? + + dominatedPoints = matRad_paretoCheckPoint(fInd,fIndv); + % if no point is dominated -> returns [] + % if the newly generated point is not optimal -> Returns array with a 0 inside + % if (a) previously generated point(s) is/are dominated by the + % newly generated point -> returns index in fInd -> Need to + % update penGrid, weights, fInd,OPSa,OPSb + + if isempty(dominatedPoints) + paretoOptimal = true; + + elseif dominatedPoints == 0 + calcPointsForFacet = calcPointsForFacet + 1; + continue %% found point is dominated by previously calculated point -> has to be recalculated + + else %newly calculated point dominates points on the previous pareto surface + paretoOptimal = true; + %remove dominated points and related equations %might + %have to check if this fully works + removedfInd = [removedfInd;fInd(dominatedPoints,:)]; + removedpenGrid = [removedpenGrid;penGrid(dominatedPoints,:)]; + removedweights = [removedweights,weights(:,dominatedPoints)]; + + + fInd(dominatedPoints,:) = []; + penGrid(dominatedPoints,:) = []; + weights(:,dominatedPoints) = []; + %Issue if + if ~(dominatedPoints-objcount <=0) + removedOPSA = [removedOPSA;OPSA(dominatedPoints-objcount,:)]; + removedOPSb = [removedOPSb;OPSb(dominatedPoints-objcount,:)]; + OPSA(dominatedPoints-objcount,:) = []; % shifted index + OPSb(dominatedPoints-objcount,:) = []; + + end + %remove dominated points + + end + end + objectiveFunctionVals(end + 1) = {optimizer.allObjectiveFunctionValues}; + errors = [errors,facetErrors(idx)]; + newPen = norm; + %awInit = wOpt; + info = optimizer.resultInfo; + weights = [weights,wOpt]; + penGrid = [penGrid;newPen]; + + + %fIndv = cellfun(@(group)sum(fIndv(:,group)),groups) + fInd = [fInd;fIndv]; + + found = true; + end + facetNum = facetNum +1; + end + + % when final point is found: Update OPSA and OPSb + OPSA = [OPSA;-newPen]; %add normal vector of facet that was run + OPSb = [OPSb;-fIndv*newPen']; + +end + + + + + +%defining a structure that contains all relevant information +returnStruct.OPSA = OPSA; +returnStruct.OPSb = OPSb; +returnStruct.weights = weights; +returnStruct.finds = fInd; +returnStruct.penGrid = penGrid; +returnStruct.errors = errors; +returnStruct.allErrors = allErrors; +returnStruct.removed = {removedfInd,removedpenGrid,removedweights,removedOPSA,removedOPSb}; + +returnStruct.weights = weights; +returnStruct.finds = fInd; +returnStruct.penGrid = penGrid; + +returnStruct.optiProb = optiProb; +returnStruct.allObj = objectiveFunctionVals; +returnStruct.modcst = cst; + +%calculate a single plan +resultGUI = matRad_calcCubes(wOpt,dij); +resultGUI.wUnsequenced = wOpt; +resultGUI.usedOptimizer = optimizer; +resultGUI.info = info; +resultGUI.optiProb = optiProb; + + + +% unblock mex files +clear mex diff --git a/matRad/optimization/pareto/navigation/matRad_ParetoSurfFromFacets.m b/matRad/optimization/pareto/navigation/matRad_ParetoSurfFromFacets.m new file mode 100644 index 000000000..e4f07bc65 --- /dev/null +++ b/matRad/optimization/pareto/navigation/matRad_ParetoSurfFromFacets.m @@ -0,0 +1,74 @@ +function [k,facets] = matRad_ParetoSurfFromFacets(fVals) +% matRad_ParetoSurfFromFacets allows to calculate the facets of the pareto surface +% from its vertices. The facets are calculated by calculating the convex hull. +% +% call +% [k,facets] = matRad_ParetoSurfFromFacets(fVals) +% +% input +% fVals: Matrix storing the objective function values of the pareto surface vertices (objective function values of the calculated points) +% +% output +% k: Matrix storing the vertices belonging to the facets of the convex hull +% facets: Matrix storing only the vertices that belong to the Pareto surface (The ones associated to facets with positive normals) +% +% References +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2023 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %calculate convex hull + [k,vol] = convhulln(fVals); + + %% + %calculate normals for each facet(normal vectors are perpendicular to their respective facet) + %Done by solving P*n= 0 where P is a matrix with the vectors spanning the hyperplane and n the + %normal vector to be calculated + + %initializing some objects that are returned by the function + facets= zeros(size(k)); + + j = 0; + %% loop over all facets of convex hull + for i = 1:size(k,1) + %% Step2: Calculate upper bounds from convex hull + %vertices of facets + ps = fVals(k(i,:),:); + %calculate vectors spanning hyperplane through refPoint + refPoint = ps(1,:); + f = ps-refPoint; % vectors spanning hyperplan through ps(1,:) + zw = f(2:end,:); % remove reference Point + % + normal = null(zw); %solve P*n=0 P Matrix + + %% check orientation of calculated normal vector + %get vertex not in facet and calculate orientation of facet + idxs = (1:size(fVals,1)); + diffs = setdiff(idxs,k(i,:)); % find vertex not in facet + + %check if orientation and normal vector face in same direction + orientationVector = fVals(diffs(1),:)-refPoint; + orientation =(orientationVector*normal>0); + + %flip orientation of normal vector if it goes in the opposite direction + normal = normal*(2*orientation-1); + normal = round(normal,4); + + %reject facet if normal has all negative components + + if any(round(normal,3)<0) + continue + end + facets(i,:) = k(i,:); + end +end \ No newline at end of file diff --git a/matRad/optimization/pareto/navigation/matRad_calcFastCubes.m b/matRad/optimization/pareto/navigation/matRad_calcFastCubes.m new file mode 100644 index 000000000..23493d7ed --- /dev/null +++ b/matRad/optimization/pareto/navigation/matRad_calcFastCubes.m @@ -0,0 +1,50 @@ +function cubes = matRad_calcFastCubes(w,dij,pln) +% matRad computation of cube for plan optimization quantity +% +% call +% resultGUI = matRad_calcFastCubes(w,dij,pln) +% +% input +% w: bixel weight vector +% dij: dose influence matrix +% pln: matRad pln struct +% +% output +% cubes: Array storing cubes of optimization quantity +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2024 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + switch pln.bioParam.quantityOpt + case 'physicalDose' + cubes = reshape(dij.physicalDose{1}*w,dij.doseGrid.dimensions); + case 'RBExD' + if isfield(dij,'mAlphaDose') && isfield(dij,'mSqrtBetaDose') + %TODO + effect = (dij.mAlphaDose{1} * w + (dij.mSqrtBetaDose{1} * w).^2); + ix = dij.bx{1}~=0 & effect(:) > 0; + cubes = zeros(dij.doseGrid.dimensions); + cubes(ix) = (sqrt(dij.ax{1}(ix).^2 + 4 .* dij.bx{1}(ix).* effect(ix)) - dij.ax{1}(ix))./(2.*dij.bx{1}(ix)); + + else + cubes = reshape(dij.physicalDose{1}*w,dij.doseGrid.dimensions)*dij.RBE; + end + case 'effect' + cubes = reshape(full(dij.mAlphaDose{1} * w + (dij.mSqrtBetaDose{1} * w).^2),dij.doseGrid.dimensions); + end + cubes = matRad_interp3(dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z, ... + cubes, ... + dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z,'linear',0); +end \ No newline at end of file diff --git a/matRad/optimization/pareto/navigation/matRad_paretoSurfaceNavigation.m b/matRad/optimization/pareto/navigation/matRad_paretoSurfaceNavigation.m new file mode 100644 index 000000000..87430d9b7 --- /dev/null +++ b/matRad/optimization/pareto/navigation/matRad_paretoSurfaceNavigation.m @@ -0,0 +1,71 @@ +function v = matRad_paretoSurfaceNavigation(Y,yr,tau,j,b) + %matRad_navigation problem allows the exploration of plans that are not in the set + %of precalculated plans. + % + % input + % Y: Set of precalculated Pareto optimal plans + % yr: Reference point that is used as a starting point to calculate the new plna + % tau: Value that objective j is fixed for calculation of next point + % j: Index of the objectives + % b: Array containing the upper bounds for each objective + % + % output + % facetPoints: vertices of the hyperplane in objective space + % refPoint: Reference Point, used to construct the hyperplane + % normal: Normal of the facet + % + % + % References + % [1] DOI: 10.1088/0031-9155/53/4/011 + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2020 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + [k,m] = size(Y); + + + f = zeros(m+k+1,1); + f(end) = 1; + + %Inequality constraints + A = zeros(k,m+k+1); + A(:,1:m) = Y; + + + + %Equality constraints + Aeq = zeros(k+1,m+k+1); + Aeq(1:k,1:m) = Y; + Aeq(k+1,1:m) = 1; + Aeq(1:k,m+1:m+k) = 1; + Aeq(1:k,end) = -1; + Aeq(j,m+1:end) = 0; + + + + beq = zeros(k+1,1); + beq(1:k,1) = yr; + beq(j,1) = tau; + beq(end) = 1; + + %Bounds + lb = zeros(m+k,1); + ub = ones(1,m); + + options = optimoptions('linprog','Display','iter'); + %options = optimoptions('linprog','Algorithm','interior-point','Display','iter'); + v = linprog(f,A,b,Aeq,beq,lb,ub,options); + if numel(v) > 0 + v = v(1:m); + end +end \ No newline at end of file diff --git a/matRad/optimization/pareto/navigation/matRad_plotPS2D.m b/matRad/optimization/pareto/navigation/matRad_plotPS2D.m new file mode 100644 index 000000000..bfd1e7ee7 --- /dev/null +++ b/matRad/optimization/pareto/navigation/matRad_plotPS2D.m @@ -0,0 +1,59 @@ +function matRad_plotPS2D(retStruct,idx) + % matRad_plotParetoFront3D implements a function to visualize a 3d Pareto + % surface + % + % + % input + % retStruct: Structure returned from Pareto optimization + % idx(optional): + % output + % resultGUI: struct containing optimized fluence vector, dose, and (for + % biological optimization) RBE-weighted dose etc. + % optimizer: Used Optimizer Object + + % + % References + % + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + ps = retStruct.finds; + if size(ps,2) < 2 + ME = MException('Matlab:WrongSize','Function requires more than one data point!'); + throw(ME); + end + if ~exist('idx','var') + if size(ps,2)>2 + ME = MException('Matlab:VariableNotFound','No indices provided!'); + throw(ME); + else %ps has exactly two objectives + idx = [1,2]; + end + else %idx variable provided + if numel(idx) ~= 2 %how many indices are provided + ME = MException('Matlab:WrongSize','Function needs exactly two indices!'); + throw(ME); + end + + if any(idx<1) || any(idx>size(ps,2)) %do the indices have the correct dimension? + ME = MException('Matlab:WrongIdx','Indices out of bounds!'); + throw(ME); + end + end + + + figure + scatter(ps(:,idx(1)),ps(:,idx(2)),'black') + xlabel('f_1[a.u]') + ylabel('f_2[a.u]') +end \ No newline at end of file diff --git a/matRad/optimization/pareto/navigation/matRad_plotParetoFront3D.m b/matRad/optimization/pareto/navigation/matRad_plotParetoFront3D.m new file mode 100644 index 000000000..38e4a826d --- /dev/null +++ b/matRad/optimization/pareto/navigation/matRad_plotParetoFront3D.m @@ -0,0 +1,40 @@ +function matRad_plotParetoFront3D(retStruct) + % matRad_plotParetoFront3D implements a function to visualize a 3d Pareto + % surface + % + % + % References + % + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2024 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + ps = retStruct.finds; + if size(ps,2) ~= 3 + ME = MException('Matlab:WrongSize','Function only supports plotting of 3D surfaces!'); + throw(ME); + end + + + [k,facets] = matRad_ParetoSurfFromFacets(ps); + figure + trisurf(facets(all(facets,2),:),ps(:,1),ps(:,2),ps(:,3),'FaceColor',[0.8 0.8 0.8]) + hold on + scatter3(ps(:,1),ps(:,2),ps(:,3),'MarkerEdgeColor','black',... + 'MarkerFaceColor',[0 0 0]) + legend('approx. surface','Location','northwest') + xlabel('f_1[a.u]') + ylabel('f_2[a.u]') + zlabel('f_3[a.u]') + view(-45,5) +end \ No newline at end of file diff --git a/matRad/optimization/pareto/tools/matRad_generateParetoAnchorPoints.m b/matRad/optimization/pareto/tools/matRad_generateParetoAnchorPoints.m new file mode 100644 index 000000000..f34aca1c5 --- /dev/null +++ b/matRad/optimization/pareto/tools/matRad_generateParetoAnchorPoints.m @@ -0,0 +1,30 @@ +function penPoints= matRad_generateParetoAnchorPoints(numObj) +% matRad function that generates the anchor points for pareto analysis +% call +% penPoints = matRad_generateSphericalPenaltyGrid(numObj) +% +% input +% numObj: Number of objectives in optimization +% +% output +% penPoints: Matrix storing the penalty vectors +% e.g in 3D +% [[1,0,0] +% [0,1,0] +% [0,0,1]] +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2023 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +penPoints = eye(numObj); + diff --git a/matRad/optimization/pareto/tools/matRad_generateParetoDummyPoints.m b/matRad/optimization/pareto/tools/matRad_generateParetoDummyPoints.m new file mode 100644 index 000000000..272536344 --- /dev/null +++ b/matRad/optimization/pareto/tools/matRad_generateParetoDummyPoints.m @@ -0,0 +1,44 @@ +function fValsMod = matRad_generateParetoDummyPoints(fVals,U) + % matRad Function that generates Dummy Points for a set of points + % + % input + % fVals: The so-far determined pareto optimal points + % U: Array storing highest objective values (should be an + % array of ones if normalized). Should be calculated after anchorpoint calculation and then kept constant to avoid issues from optimization complications + % + % output + % fValsMod: Modified pareto set to remove mixed normal facets (see + % Rennen) + % + % References + % - https://dx.doi.org/10.2139/ssrn.1427721 + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2023 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +theta = 0.1; +m = size(fVals,2); +n = size(fVals,1); + +U = U*m+theta; %modifying upper bound + +fValsMod = zeros(size(fVals,1)*size(fVals,2),size(fVals,2)); %initialize + +for i = 1:n %loop over and generate a dummy point for each dimension + temp = repmat(fVals(i,:),m,1); + for j = 1:numel(U) + temp(j,j) = U(j); + end + fValsMod((i-1)*m+1:i*m,:) = temp; +end + +fValsMod = [fVals;fValsMod]; \ No newline at end of file diff --git a/matRad/optimization/pareto/tools/matRad_normalFromFacet.m b/matRad/optimization/pareto/tools/matRad_normalFromFacet.m new file mode 100644 index 000000000..765396f9a --- /dev/null +++ b/matRad/optimization/pareto/tools/matRad_normalFromFacet.m @@ -0,0 +1,52 @@ +function [facetPoints,refPoint,normal] = matRad_normalFromFacet(fVals,k,i) + % matRad helper function that calculates the normal for a facet + % + % input + % fVals: All input values of the convexHull + % k: Array with the vertex indices of the facets + % i: Index of the facet + % + % output + % facetPoints: vertices of the hyperplane in objective space + % refPoint: Reference Point, used to construct the hyperplane + % normal: Normal of the facet + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % + % Copyright 2023 the matRad development team. + % + % This file is part of the matRad project. It is subject to the license + % terms in the LICENSE file found in the top-level directory of this + % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part + % of the matRad project, including this file, may be copied, modified, + % propagated, or distributed except according to the terms contained in the + % LICENSE file. + % + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + facetPoints = fVals(k(i,:),:); %points that span the given facet + %choose a reference point of facet that hyperplane is "build on" + refPoint = facetPoints(1,:); + hyperplaneVectors = facetPoints-refPoint; % calculate difference of reference point to all other points on facet + spanningVectors = hyperplaneVectors(2:end,:); %reference point not needed for calculation of normal + + %calculate the normal of the hyperplane by solving V*n = 0 where V is + %the matrix with the vectors spanning the hyperplane + normal = null(spanningVectors); + + %we want to check if the components of the normal contain negative + %components. For that one has to check the orientation first + if rank(spanningVectors) ~= size(spanningVectors,1) + normal = -1*ones(size(spanningVectors,2),1); %simply return all -1 + return + end + %choose a point of the PS that is not part of the hyperplane currently + %investigated + idxs = (1:size(fVals,1)); + + %calculate orientation vector + orientationVector = mean(fVals(unique(k),:))-refPoint; + orientation = (orientationVector*normal>0); %check orientation of facet (either 0 or 1) + normal = normal*(2*orientation-1); %flip normal vector if it faces in the wrong direction + \ No newline at end of file diff --git a/matRad/optimization/pareto/tools/matRad_paretoCheckPoint.m b/matRad/optimization/pareto/tools/matRad_paretoCheckPoint.m new file mode 100644 index 000000000..b29e3e765 --- /dev/null +++ b/matRad/optimization/pareto/tools/matRad_paretoCheckPoint.m @@ -0,0 +1,40 @@ +function dominatedPoints = matRad_paretoCheckPoint(fVals,newfVal) +% matRad helper function that checks if a newly generated point dominates +% another point or is dominated by a point in the previous sets. should +% be called after calculating the first initial points +% +% +% input +% fVals: The so-far determined pareto optimal points +% newfVal: New point to check against fVals +% +% output +% dominatedPoints: 0 if new point dominated by point in previous set +% empty if new point neither dominates a point in old set or is dominated by point in old set +% otherwise array refering to the indices of points +% in the old set dominated by the newly calculated +% point +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2023 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +dominatedPoints = []; +for i = 1:size(fVals,1) + if all(newfVal > fVals(i,:)) % newly generated point is not pareto optimal + dominatedPoints = [dominatedPoints,0]; + break + elseif all(newfVal < fVals(i,:)) % point i is dominated by newly generated point + dominatedPoints = [dominatedPoints,i]; + end +end \ No newline at end of file diff --git a/matRad/optimization/pareto/tools/matRad_paretoSurfFromFacets.m b/matRad/optimization/pareto/tools/matRad_paretoSurfFromFacets.m new file mode 100644 index 000000000..28cf12379 --- /dev/null +++ b/matRad/optimization/pareto/tools/matRad_paretoSurfFromFacets.m @@ -0,0 +1,74 @@ +function [k,facets] = matRad_paretoSurfFromFacets(fVals) +% matRad_paretoSurfFromFacets allows to calculate the facets of the pareto surface +% from its vertices. The facets are calculated by calculating the convex hull. +% +% call +% [k,facets] = matRad_paretoSurfFromFacets(fVals) +% +% input +% fVals: Matrix storing the objective function values of the pareto surface vertices (objective function values of the calculated points) +% +% output +% k: Matrix storing the vertices belonging to the facets of the convex hull +% facets: Matrix storing only the vertices that belong to the Pareto surface (The ones associated to facets with positive normals) +% +% References +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2023 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %calculate convex hull + [k,vol] = convhulln(fVals); + + %% + %calculate normals for each facet(normal vectors are perpendicular to their respective facet) + %Done by solving P*n= 0 where P is a matrix with the vectors spanning the hyperplane and n the + %normal vector to be calculated + + %initializing some objects that are returned by the function + facets= zeros(size(k)); + + j = 0; + %% loop over all facets of convex hull + for i = 1:size(k,1) + %% Step2: Calculate upper bounds from convex hull + %vertices of facets + ps = fVals(k(i,:),:); + %calculate vectors spanning hyperplane through refPoint + refPoint = ps(1,:); + f = ps-refPoint; % vectors spanning hyperplan through ps(1,:) + zw = f(2:end,:); % remove reference Point + % + normal = null(zw); %solve P*n=0 P Matrix + + %% check orientation of calculated normal vector + %get vertex not in facet and calculate orientation of facet + idxs = (1:size(fVals,1)); + diffs = setdiff(idxs,k(i,:)); % find vertex not in facet + + %check if orientation and normal vector face in same direction + orientationVector = fVals(diffs(1),:)-refPoint; + orientation =(orientationVector*normal>0); + + %flip orientation of normal vector if it goes in the opposite direction + normal = normal*(2*orientation-1); + normal = round(normal,4); + + %reject facet if normal has all negative components + + if any(round(normal,3)<0) + continue + end + facets(i,:) = k(i,:); + end +end \ No newline at end of file diff --git a/matRad/planAnalysis/matRad_showDVH.m b/matRad/planAnalysis/matRad_showDVH.m index f61094a5e..f08bdefde 100644 --- a/matRad/planAnalysis/matRad_showDVH.m +++ b/matRad/planAnalysis/matRad_showDVH.m @@ -96,7 +96,7 @@ function matRad_showDVH(axesHandle,dvh,cst,pln,lineStyleIndicator) if exist('pln','var') && ~isempty(pln) if strcmp(pln.bioParam.model,'none') - xlabel('Dose [Gy]','FontSize',fontSizeValue); + xlabel(axesHandle,'Dose [Gy]','FontSize',fontSizeValue); else xlabel(axesHandle,'RBE x Dose [Gy(RBE)]','FontSize',fontSizeValue); end diff --git a/matRad/plotting/matRad_plotParetoSurface.m b/matRad/plotting/matRad_plotParetoSurface.m new file mode 100644 index 000000000..33d9e4ecb --- /dev/null +++ b/matRad/plotting/matRad_plotParetoSurface.m @@ -0,0 +1,44 @@ +function matRad_plotParetoSurface(retStruct) +% matRad function that plots a color coded Pareto Surface. Colors are based +% on penalty values of the data points in 3 dimensions. +% Red means higher penalties of objective plotted along x-axis. +% Green means higher penalties of objective plotted along y-axis. +% Blue means higher penalties of objective plotted along z-axis. +% +% call +% cBarHandle = matRad_plotColorbar(axesHandle,cMap,window,varargin) +% +% input +% axesHandle handle to axes the colorbar will be attached to +% cMap corresponding colormap +% window colormap window (corresponds to clim) +% varargin additional key-value pairs that will be forwarded to the +% MATLAB colorbar(__) call +% +% output +% cBarHandle handle of the colorbar object +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +fig2 = figure; +ps = retStruct.finds; +[k,facets] = matRad_ParetoSurfFromFacets(ps); + +trisurf(facets(all(facets,2),:),ps(:,1),ps(:,2),ps(:,3),'FaceColor',[0.8 0.8 0.8]) +hold on +scatter3(ps(:,1),ps(:,2),ps(:,3),'MarkerEdgeColor','black',... + 'MarkerFaceColor',[0 0 0]) diff --git a/matRad/plotting/matRad_plotPenaltyPoints.m b/matRad/plotting/matRad_plotPenaltyPoints.m new file mode 100644 index 000000000..732a0786b --- /dev/null +++ b/matRad/plotting/matRad_plotPenaltyPoints.m @@ -0,0 +1,48 @@ +function matRad_plotPenaltyPoints(penPoints) +% matRad function to plot the penalty Vectors - works only for 2, 3 or 4 +% dimensional objectives (latter case utilizes the normalization of the objectives) +% +% +% input +% penPoints Matrix storing the penalty vectors +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2023 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +switch size(penPoints,2) + case 2 + figure + scatter(penPoints(:,1),penPoints(:,2),'filled','MarkerFaceColor',[0.4660 0.6740 0.1880]) + xlabel('Penalty 1'); + ylabel('Penalty 2'); + case 3 + figure + scatter3(penPoints(:,1),penPoints(:,2),penPoints(:,3),[], penPoints(:,3),'filled') + colormap(gca,"summer") + xlabel('Penalty 1'); + ylabel('Penalty 2'); + zlabel('Penalty 3'); + case 4 + figure + scatter3(penPoints(:,1),penPoints(:,2),penPoints(:,3),[], penPoints(:,4),'filled') + h = colorbar() + colormap(gca,"summer") + xlabel('Penalty 1'); + ylabel('Penalty 2'); + zlabel('Penalty 3'); + set(get(h,'title'),'string','Penalty 4'); + otherwise + warning(['Number of objectives for Pareto Analysis not suited for Plot!']); +end diff --git a/matRad/tools/matRad_plotSliceWrapper.m b/matRad/tools/matRad_plotSliceWrapper.m index f166574ee..2657b2326 100644 --- a/matRad/tools/matRad_plotSliceWrapper.m +++ b/matRad/tools/matRad_plotSliceWrapper.m @@ -102,7 +102,6 @@ set(axesHandle,'YDir','Reverse'); % plot ct slice hCt = matRad_plotCtSlice(axesHandle,ct.cubeHU,cubeIdx,plane,slice); -hold on; % plot dose if ~isempty(doseWindow) && doseWindow(2) - doseWindow(1) <= 0 diff --git a/matRad/util/matRad_calcCubes.m b/matRad/util/matRad_calcCubes.m index d63b5e12a..56dcbc3f2 100644 --- a/matRad/util/matRad_calcCubes.m +++ b/matRad/util/matRad_calcCubes.m @@ -29,148 +29,12 @@ % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - if nargin < 3 scenNum = 1; end -resultGUI.w = w; - -if isfield(dij,'numParticlesPerMU') - resultGUI.MU = (w.*1e6) ./ dij.numParticlesPerMU; -end - -% get bixel - beam correspondence -for i = 1:dij.numOfBeams - beamInfo(i).suffix = ['_beam', num2str(i)]; - beamInfo(i).logIx = (dij.beamNum == i); -end -beamInfo(dij.numOfBeams+1).suffix = ''; -beamInfo(dij.numOfBeams+1).logIx = true(size(resultGUI.w,1),1); - -[ctScen,~] = ind2sub(size(dij.physicalDose),scenNum); - - -%% Physical Dose -doseFields = {'physicalDose','doseToWater'}; -doseQuantities = {'','_std','_batchStd'}; -% compute physical dose for all beams individually and together -for j = 1:length(doseFields) - for k = 1:length(doseQuantities) - % Check if combination is a field in dij, otherwise skip - if isfield(dij,[doseFields{j} doseQuantities{k}]) - % Handle standard deviation fields and add quadratically - if ~isempty(strfind(lower(doseQuantities{1}),'std')) - for i = 1:length(beamInfo) - resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix]) = sqrt(reshape(full(dij.([doseFields{j} doseQuantities{k}]){scenNum}.^2 * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions)); - resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix])(isnan(resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix]))) = 0; - end - % Handle normal fields as usual - else - for i = 1:length(beamInfo) - resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix]) = reshape(full(dij.([doseFields{j} doseQuantities{k}]){scenNum} * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); - end - end - end - end -end - - -%% LET -% consider LET -if isfield(dij,'mLETDose') - for i = 1:length(beamInfo) - LETDoseCube = reshape(full(dij.mLETDose{scenNum} * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); - resultGUI.(['LET', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); - ix = resultGUI.(['physicalDose', beamInfo(i).suffix]) > 0; - resultGUI.(['LET', beamInfo(i).suffix])(ix) = LETDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix); - end -end - - -%% RBE weighted dose -% consider RBE for protons and skip varRBE calculation -if isfield(dij,'RBE') && isscalar(dij.RBE) - for i = 1:length(beamInfo) - resultGUI.(['RBExD', beamInfo(i).suffix]) = resultGUI.(['physicalDose', beamInfo(i).suffix]) * dij.RBE; - end -elseif any(cellfun(@(teststr) ~isempty(strfind(lower(teststr),'alpha')), fieldnames(dij))) - % Load RBE models if MonteCarlo was calculated for multiple models - if isfield(dij,'RBE_models') - RBE_model = cell(1,length(dij.RBE_models)); - for i = 1:length(dij.RBE_models) - RBE_model{i} = ['_' dij.RBE_models{i}]; - end - else - RBE_model = {''}; - end - - % Loop through RBE models - for j = 1:length(RBE_model) - % Check if combination is a field in dij, otherwise skip - if isfield(dij,['mAlphaDose' RBE_model{j}]) - for i = 1:length(beamInfo) - % Get weights of current beam - wBeam = (resultGUI.w .* beamInfo(i).logIx); - - % consider biological optimization - ix = dij.bx{ctScen}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0; - - % Calculate effect from alpha- and sqrtBetaDose - resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]) = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam + (dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam).^2); - resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]) = reshape(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]),dij.doseGrid.dimensions); - - % Calculate RBExD from the effect - resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix]) = zeros(size(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]))); - resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix) = (sqrt(dij.ax{ctScen}(ix).^2 + 4 .* dij.bx{ctScen}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{ctScen}(ix))./(2.*dij.bx{ctScen}(ix)); - - % Divide RBExD with the physicalDose to get the plain RBE cube - resultGUI.(['RBE', RBE_model{j}, beamInfo(i).suffix]) = resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])./resultGUI.(['physicalDose', beamInfo(i).suffix]); - - % Initialize alpha/beta cubes - resultGUI.(['alpha', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); - resultGUI.(['beta', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); - resultGUI.(['alphaDoseCube', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); - resultGUI.(['SqrtBetaDoseCube', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); - - % Calculate alpha and weighted alphaDose - AlphaDoseCube = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam); - resultGUI.(['alpha', RBE_model{j}, beamInfo(i).suffix])(ix) = AlphaDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix); - resultGUI.(['alphaDoseCube', RBE_model{j}, beamInfo(i).suffix])(ix) = AlphaDoseCube(ix); - - % Calculate beta and weighted sqrtBetaDose - SqrtBetaDoseCube = full(dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam); - resultGUI.(['beta', RBE_model{j}, beamInfo(i).suffix])(ix) = (SqrtBetaDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix)).^2; - resultGUI.(['SqrtBetaDoseCube', RBE_model{j}, beamInfo(i).suffix])(ix) = SqrtBetaDoseCube(ix); - end - end - end -end - -%add some dij meta -if isfield(dij,'meta') && isstruct(dij.meta) - resultGUI.meta = dij.meta; -end - - -%% Final processing -% Remove suffix for RBExD if there's only one available -if any(cellfun(@(teststr) ~isempty(strfind(lower(teststr),'alpha')), fieldnames(dij))) && isfield(dij,'RBE_models') && length(dij.RBE_models) == 1 - % Get fieldnames that include the specified RBE model - fnames = fieldnames(resultGUI); - fnames = fnames(cellfun(@(teststr) ~isempty(strfind(lower(teststr),lower(dij.RBE_models{1}))), fnames)); - - % Rename fields and remove model specifier if there's only one - for f = 1:length(fnames) - resultGUI.(erase(fnames{f},['_',dij.RBE_models{1}])) = resultGUI.(fnames{f}); - end - - % Remove old fields - resultGUI = rmfield(resultGUI,fnames); -end - -% group similar fields together -resultGUI = orderfields(resultGUI); +% First calculate on dose grid +resultGUI = matRad_calcCubesDoseGrid(w,dij,scenNum); % interpolation if dose grid does not match ct grid if isfield(dij,'ctGrid') && any(dij.ctGrid.dimensions~=dij.doseGrid.dimensions) diff --git a/matRad/util/matRad_calcCubesDoseGrid.m b/matRad/util/matRad_calcCubesDoseGrid.m new file mode 100644 index 000000000..ce51341df --- /dev/null +++ b/matRad/util/matRad_calcCubesDoseGrid.m @@ -0,0 +1,172 @@ +function resultGUI = matRad_calcCubesDoseGrid(w,dij,scenNum) +% matRad computation of all cubes for the resultGUI struct +% which is used as result container and for visualization in matRad's GUI +% +% call +% resultGUI = matRad_calcCubes(w,dij) +% resultGUI = matRad_calcCubes(w,dij,scenNum) +% +% input +% w: bixel weight vector +% dij: dose influence matrix +% scenNum: optional: number of scenario to calculated (default 1) +% +% output +% resultGUI: matRad result struct +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2015 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargin < 3 + scenNum = 1; +end + +resultGUI.w = w; + +if isfield(dij,'numParticlesPerMU') + resultGUI.MU = (w.*1e6) ./ dij.numParticlesPerMU; +end + +% get bixel - beam correspondence +for i = 1:dij.numOfBeams + beamInfo(i).suffix = ['_beam', num2str(i)]; + beamInfo(i).logIx = (dij.beamNum == i); +end +beamInfo(dij.numOfBeams+1).suffix = ''; +beamInfo(dij.numOfBeams+1).logIx = true(size(resultGUI.w,1),1); + +[ctScen,~] = ind2sub(size(dij.physicalDose),scenNum); + + +%% Physical Dose +doseFields = {'physicalDose','doseToWater'}; +doseQuantities = {'','_std','_batchStd'}; +% compute physical dose for all beams individually and together +for j = 1:length(doseFields) + for k = 1:length(doseQuantities) + % Check if combination is a field in dij, otherwise skip + if isfield(dij,[doseFields{j} doseQuantities{k}]) + % Handle standard deviation fields and add quadratically + if ~isempty(strfind(lower(doseQuantities{1}),'std')) + for i = 1:length(beamInfo) + resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix]) = sqrt(reshape(full(dij.([doseFields{j} doseQuantities{k}]){scenNum}.^2 * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions)); + resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix])(isnan(resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix]))) = 0; + end + % Handle normal fields as usual + else + for i = 1:length(beamInfo) + resultGUI.([doseFields{j}, doseQuantities{k}, beamInfo(i).suffix]) = reshape(full(dij.([doseFields{j} doseQuantities{k}]){scenNum} * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); + end + end + end + end +end + + +%% LET +% consider LET +if isfield(dij,'mLETDose') + for i = 1:length(beamInfo) + LETDoseCube = reshape(full(dij.mLETDose{scenNum} * (resultGUI.w .* beamInfo(i).logIx)),dij.doseGrid.dimensions); + resultGUI.(['LET', beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + ix = resultGUI.(['physicalDose', beamInfo(i).suffix]) > 0; + resultGUI.(['LET', beamInfo(i).suffix])(ix) = LETDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix); + end +end + + +%% RBE weighted dose +% consider RBE for protons and skip varRBE calculation +if isfield(dij,'RBE') && isscalar(dij.RBE) + for i = 1:length(beamInfo) + resultGUI.(['RBExD', beamInfo(i).suffix]) = resultGUI.(['physicalDose', beamInfo(i).suffix]) * dij.RBE; + end +elseif any(cellfun(@(teststr) ~isempty(strfind(lower(teststr),'alpha')), fieldnames(dij))) + % Load RBE models if MonteCarlo was calculated for multiple models + if isfield(dij,'RBE_models') + RBE_model = cell(1,length(dij.RBE_models)); + for i = 1:length(dij.RBE_models) + RBE_model{i} = ['_' dij.RBE_models{i}]; + end + else + RBE_model = {''}; + end + + % Loop through RBE models + for j = 1:length(RBE_model) + % Check if combination is a field in dij, otherwise skip + if isfield(dij,['mAlphaDose' RBE_model{j}]) + for i = 1:length(beamInfo) + % Get weights of current beam + wBeam = (resultGUI.w .* beamInfo(i).logIx); + + % consider biological optimization + ix = dij.bx{ctScen}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0; + + % Calculate effect from alpha- and sqrtBetaDose + resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]) = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam + (dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam).^2); + resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]) = reshape(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]),dij.doseGrid.dimensions); + + % Calculate RBExD from the effect + resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix]) = zeros(size(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]))); + resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix) = (sqrt(dij.ax{ctScen}(ix).^2 + 4 .* dij.bx{ctScen}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{ctScen}(ix))./(2.*dij.bx{ctScen}(ix)); + + % Divide RBExD with the physicalDose to get the plain RBE cube + resultGUI.(['RBE', RBE_model{j}, beamInfo(i).suffix]) = resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])./resultGUI.(['physicalDose', beamInfo(i).suffix]); + + % Initialize alpha/beta cubes + resultGUI.(['alpha', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + resultGUI.(['beta', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + resultGUI.(['alphaDoseCube', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + resultGUI.(['SqrtBetaDoseCube', RBE_model{j}, beamInfo(i).suffix]) = zeros(dij.doseGrid.dimensions); + + % Calculate alpha and weighted alphaDose + AlphaDoseCube = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam); + resultGUI.(['alpha', RBE_model{j}, beamInfo(i).suffix])(ix) = AlphaDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix); + resultGUI.(['alphaDoseCube', RBE_model{j}, beamInfo(i).suffix])(ix) = AlphaDoseCube(ix); + + % Calculate beta and weighted sqrtBetaDose + SqrtBetaDoseCube = full(dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam); + resultGUI.(['beta', RBE_model{j}, beamInfo(i).suffix])(ix) = (SqrtBetaDoseCube(ix)./resultGUI.(['physicalDose', beamInfo(i).suffix])(ix)).^2; + resultGUI.(['SqrtBetaDoseCube', RBE_model{j}, beamInfo(i).suffix])(ix) = SqrtBetaDoseCube(ix); + end + end + end +end + +%% Final processing +% Remove suffix for RBExD if there's only one available +if any(cellfun(@(teststr) ~isempty(strfind(lower(teststr),'alpha')), fieldnames(dij))) && isfield(dij,'RBE_models') && length(dij.RBE_models) == 1 + % Get fieldnames that include the specified RBE model + fnames = fieldnames(resultGUI); + fnames = fnames(cellfun(@(teststr) ~isempty(strfind(lower(teststr),lower(dij.RBE_models{1}))), fnames)); + + % Rename fields and remove model specifier if there's only one + for f = 1:length(fnames) + resultGUI.(erase(fnames{f},['_',dij.RBE_models{1}])) = resultGUI.(fnames{f}); + end + + % Remove old fields + resultGUI = rmfield(resultGUI,fnames); +end + +%add some dij meta +if isfield(dij,'meta') && isstruct(dij.meta) + resultGUI.meta = dij.meta; +end + +% group similar fields together +resultGUI = orderfields(resultGUI); \ No newline at end of file