e0404 · wahln · Oct 27, 2016 · Oct 27, 2016 · Oct 27, 2016 · Oct 27, 2016
diff --git a/basedata/photons_Generic.mat b/basedata/photons_Generic.mat
diff --git a/dicom/matRad_importDicomSteeringPhotons.m b/dicom/matRad_importDicomSteeringPhotons.m
@@ -48,7 +48,15 @@
     % set necessary steering information 
     stf(i).gantryAngle  = UniqueComb(i,1);
     stf(i).couchAngle   = UniqueComb(i,2);
-    stf(i).isoCenter    = pln.propStf.isoCenter(i,:);
+
+    %Handle possibility of multiple isocenters
+    if size(pln.propStf.isoCenter,1) == 1
+        stf(i).isoCenter    = pln.propStf.isoCenter;
+    elseif size(pln.propStf.isoCenter,1) == pln.propStf.numOfBeams
+        stf(i).isoCenter = pln.propStf.isoCenter(i,:);
+    else
+        error('Invalid number of isocenters - should either be one or as many as beams!');
+    end
 
     % bixelWidth = 'field' as keyword for whole field dose calc
     stf(i).bixelWidth    = 'field';

diff --git a/examples/matRad_example3_photonsDAO.m b/examples/matRad_example3_photonsDAO.m
@@ -66,6 +66,7 @@
 % Enable sequencing and direct aperture optimization (DAO).
 pln.propOpt.runSequencing = 1;
 pln.propOpt.runDAO        = 1;
+pln.propOpt.numLevels = 5;
 
 %% Generate Beam Geometry STF
 stf = matRad_generateStf(ct,cst,pln);
@@ -90,7 +91,7 @@
 % order to modulate the intensity of the beams with multiple static 
 % segments, so that translates each intensity map into a set of deliverable 
 % aperture shapes.
-resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,5);
+resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,pln);
 
 %% DAO - Direct Aperture Optimization
 % The Direct Aperture Optimization is an optimization approach where we 

diff --git a/examples/matRad_example8_photonsVMAT.m b/examples/matRad_example8_photonsVMAT.m
@@ -0,0 +1,124 @@
+%% Example Photon Treatment Plan with VMAT direct aperture optimization
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% 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 input necessary parameters in the pln structure
+% (iii) how to setup a photon dose calculation
+% (iv) how to inversely optimize fluence directly from command window in MatLab.
+% (v) how to apply a sequencing algorithm
+% (vi) how to run a VMAT direct aperture optimization
+% (vii) how to visually and quantitatively evaluate the result
+
+%% Patient Data Import
+% Let's begin with a clear Matlab environment and import the head &
+% neck patient into your workspace.
+matRad_rc
+
+load('HEAD_AND_NECK.mat');
+
+%% Treatment Plan
+% The next step is to define your treatment plan labeled as 'pln'. This 
+% structure requires input from the treatment planner and defines 
+% the most important cornerstones of your treatment plan.
+
+pln.radiationMode   = 'photons';   % either photons / protons / carbon
+pln.machine         = 'Generic';
+
+pln.numOfFractions  = 30;
+
+% beam geometry settings
+pln.propStf.bixelWidth = 5;
+
+% optimization settings
+pln.propOpt.bioOptimization = 'none';
+pln.propOpt.runVMAT         = true;
+pln.propOpt.runDAO          = true;
+pln.propOpt.runSequencing   = true;
+pln.propOpt.preconditioner  = true;
+pln.propOpt.numLevels       = 7;
+
+pln.propOpt.VMAToptions.machineConstraintFile = [pln.radiationMode '_' pln.machine];
+
+pln.propOpt.VMAToptions.maxGantryAngleSpacing    = 2;      % Max gantry angle spacing for dose calculation
+pln.propOpt.VMAToptions.maxDAOGantryAngleSpacing = 4;      % Max gantry angle spacing for DAO
+pln.propOpt.VMAToptions.maxFMOGantryAngleSpacing = 28;     % Max gantry angle spacing for FMO
+
+pln.propOpt.VMAToptions.startingAngle = 0;
+pln.propOpt.VMAToptions.finishingAngle = 359;
+pln.propOpt.VMAToptions.continuousAperture = 0;
+
+pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]
+pln.propDoseCalc.doseGrid.resolution.y = 5; % [mm]
+pln.propDoseCalc.doseGrid.resolution.z = 5; % [mm]
+
+%%
+% Generate dose calculation, DAO, and FMO angles from the parameters input
+% above. FMO is performed only on the initGantryAngles set. In the DAO
+% step, weights and leaf positions are optimized at the angles in the
+% optGantryAngles set. Weights and leaf positions are interpolated at the
+% angles in the gantryAngles set to increase the accuracy of the dose
+% calculation (each iteration).
+
+% FMO: optimize fluence on coarse subset of gantry angles
+% Sequencing: select subset of apertures and spread to finer angles
+% DAO: constrain for leaf speed, gantry rotation speed and MU rate
+
+pln = matRad_VMATGantryAngles(pln,cst,ct);
+
+%% Generate Beam Geometry STF
+stf = matRad_generateStf(ct,cst,pln);
+
+%% Dose Calculation
+% Lets generate dosimetric information by pre-computing dose influence 
+% matrices for unit beamlet intensities. Having dose influences available 
+% allows for subsequent inverse optimization.
+dij = matRad_calcPhotonDose(ct,stf,pln,cst);
+
+%% Inverse Planning for IMRT
+% The goal of the fluence optimization is to find a set of beamlet weights 
+% which yield the best possible dose distribution according to the 
+% predefined clinical objectives and constraints underlying the radiation 
+% treatment. In VMAT, FMO is done only at the angles in the
+% initGantryAngles set. Once the optimization has finished, trigger once the GUI to
+% visualize the optimized dose cubes.
+resultGUI = matRad_fluenceOptimization(dij,cst,pln,stf);
+matRadGUI;
+
+%% Sequencing
+% This is a multileaf collimator leaf sequencing algorithm that is used in 
+% order to modulate the intensity of the beams with multiple static 
+% segments, so that translates each intensity map into a set of deliverable 
+% aperture shapes. The fluence map at each angle in the initGantryAngles
+% set is sequenced, with the resulting apertures spread to neighbouring
+% angles from the optGantryAngles set.
+resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,pln,0);
+
+%% DAO - Direct Aperture Optimization
+% The Direct Aperture Optimization is an optimization approach where we 
+% directly optimize aperture shapes and weights at the angles in the
+% optGantryAngles set.  The gantry angle speed, leaf speed, and MU rate are
+% constrained by the min and max values specified by the user.
+resultGUI = matRad_directApertureOptimization(dij,cst,resultGUI.apertureInfo,resultGUI,pln);
+
+%% Aperture visualization
+% Use a matrad function to visualize the resulting aperture shapes
+matRad_visApertureInfo(resultGUI.apertureInfo);
+
+%% Indicator Calculation and display of DVH and QI
+[dvh,qi] = matRad_indicatorWrapper(cst,pln,resultGUI);
+matRad_showDVH(dvh,cst,pln);
+
diff --git a/matRad.m b/matRad.m
@@ -32,22 +32,26 @@
 
 % beam geometry settings
 pln.propStf.bixelWidth      = 5; % [mm] / also corresponds to lateral spot spacing for particles
-pln.propStf.gantryAngles    = [0:72:359]; % [?]
+pln.propStf.gantryAngles    = 0:72:359; % [?]
 pln.propStf.couchAngles     = [0 0 0 0 0]; % [?]
 pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
 pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
 
-% dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]
 pln.propDoseCalc.doseGrid.resolution.y = 5; % [mm]
 pln.propDoseCalc.doseGrid.resolution.z = 5; % [mm]
 
+
 % optimization settings
 pln.propOpt.optimizer       = 'IPOPT';
 pln.propOpt.bioOptimization = 'none'; % none: physical optimization;             const_RBExD; constant RBE of 1.1;
-                                      % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of RBE-weighted dose
+                                            % LEMIV_effect: effect-based optimization; LEMIV_RBExD: optimization of RBE-weighted dose
 pln.propOpt.runDAO          = false;  % 1/true: run DAO, 0/false: don't / will be ignored for particles
-pln.propOpt.runSequencing   = false;  % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below
+pln.propOpt.runSequencing   = true;  % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below
+pln.propOpt.runVMAT         = false;  % 1/true: run VMAT, 0/false: don't
+pln.propOpt.preconditioner  = false;
+pln.propOpt.numLevels = 12;
+
 
 %% initial visualization and change objective function settings if desired
 matRadGUI
@@ -58,27 +62,33 @@
 %% dose calculation
 if strcmp(pln.radiationMode,'photons')
     dij = matRad_calcPhotonDose(ct,stf,pln,cst);
-    %dij = matRad_calcPhotonDoseVmc(ct,stf,pln,cst);
 elseif strcmp(pln.radiationMode,'protons') || strcmp(pln.radiationMode,'carbon')
     dij = matRad_calcParticleDose(ct,stf,pln,cst);
 end
 
+
 %% inverse planning for imrt
 resultGUI = matRad_fluenceOptimization(dij,cst,pln);
 
 %% sequencing
 if strcmp(pln.radiationMode,'photons') && (pln.propOpt.runSequencing || pln.propOpt.runDAO)
     %resultGUI = matRad_xiaLeafSequencing(resultGUI,stf,dij,5);
-    %resultGUI = matRad_engelLeafSequencing(resultGUI,stf,dij,5);
-    resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,5);
+    %resultGUI = matRad_engelLeafSequencing(resultGUI,stf,dij,pln,5);
+    resultGUI = matRad_siochiLeafSequencing(resultGUI,stf,dij,pln,0);
+    %resultGUI = matRad_svenssonLeafSequencing(resultGUI,stf,dij,pln,0);
 end
 
 %% DAO
 if strcmp(pln.radiationMode,'photons') && pln.propOpt.runDAO
-   resultGUI = matRad_directApertureOptimization(dij,cst,resultGUI.apertureInfo,resultGUI,pln);
-   matRad_visApertureInfo(resultGUI.apertureInfo);
+    resultGUI = matRad_directApertureOptimization(dij,cst,resultGUI.apertureInfo,resultGUI,pln);
+
+    if pln.propOpt.runVMAT
+        resultGUI = matRad_calcDeliveryMetrics(resultGUI,pln,stf);
+    end
+    matRad_visApertureInfo(resultGUI.apertureInfo);
 end
 
+
 %% start gui for visualization of result
 matRadGUI
 

diff --git a/matRadGUI.m b/matRadGUI.m
@@ -1645,11 +1645,11 @@ function btnOptimize_Callback(hObject, eventdata, handles)
     %% sequencing
     if strcmp(pln.radiationMode,'photons') && (pln.propOpt.runSequencing || pln.propOpt.runDAO)
     %   resultGUI = matRad_xiaLeafSequencing(resultGUI,evalin('base','stf'),evalin('base','dij')...
-    %       ,get(handles.editSequencingLevel,'Value'));
+    %       ,pln);
     %   resultGUI = matRad_engelLeafSequencing(resultGUI,evalin('base','stf'),evalin('base','dij')...
-    %       ,str2double(get(handles.editSequencingLevel,'String')));
+    %       ,pln);
         resultGUI = matRad_siochiLeafSequencing(resultGUI,evalin('base','stf'),evalin('base','dij')...
-            ,str2double(get(handles.editSequencingLevel,'String')));
+            ,pln);
 
         assignin('base','resultGUI',resultGUI);
     end
@@ -2062,6 +2062,9 @@ function setPln(handles)
     set(handles.txtSequencing,'Enable','on');
     set(handles.radiobutton3Dconf,'Enable','on');
     set(handles.editSequencingLevel,'Enable','on');
+    if isfield(pln.propOpt,'numLevels')
+        set(handles.editSequencingLevel,'Value',pln.propOpt.numLevels);
+    end
 else
     set(handles.txtSequencing,'Enable','off');
     set(handles.radiobutton3Dconf,'Enable','off');
@@ -2389,6 +2392,15 @@ function editFraction_Callback(hObject, ~, handles)
 % text box: stratification levels
 function editSequencingLevel_Callback(~, ~, ~)
 
+pln = evalin('base','pln');
+levels = str2double(get(hObject,'String'));
+
+pln.propOpt.numLevels = levels;
+
+assignin('base','pln',pln);
+guidata(hObject,handles);
+
+
 % text box: isoCenter in [mm]
 function editIsoCenter_Callback(hObject, ~, handles)
 

diff --git a/matRad_DijSampling.m b/matRad_DijSampling.m
@@ -77,12 +77,12 @@
     ixCore      = bixelDose > relDoseThreshold * max(bixelDose); % get voxels indices having a greater dose than the thresholdDose
 end
 
-bixelDoseCore       = bixelDose(ixCore);                         % save dose values that are not affected by sampling
+bixelDoseCore = bixelDose(ixCore);                         % save dose values that are not affected by sampling
 
 if all(ixCore)
     %% all bixels are in the core
     %exit function with core dose only
-    ixNew = ix;
+    ixNew        = ix;
     bixelDoseNew = bixelDoseCore;
 else
     logIxTail           = ~ixCore;                                   % get voxels indices beyond r0
@@ -126,6 +126,9 @@
     % cut new vectors and add inner core values
     ixNew        = [ix(ixCore);    ixNew(1:IxCnt-1)];
     bixelDoseNew = [bixelDoseCore; bixelDoseNew(1:IxCnt-1)];
+
+    [ixNew, ind] = sort(ixNew);
+    bixelDoseNew = bixelDoseNew(ind);
 end