Pareto navigation

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)

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)

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;

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();
+
+
+
+
+
+
+

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
+
+