Class based Sequencing

examples/matRad_example10_4DphotonRobust.m

@@ -177,8 +177,6 @@
 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.propOpt.runDAO        = 0;
-pln.propOpt.runSequencing = 0;
 
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]
@@ -231,7 +229,8 @@
 totalPhaseMatrix = ones(dij.totalNumOfBixels,ct.numOfCtScen)/ct.numOfCtScen;  % the total phase matrix determines a mapping what fluence will be delivered in the which phase
 totalPhaseMatrix = bsxfun(@times,totalPhaseMatrix,resultGUIrobust.w);         % equally distribute the fluence over all fluences
 
-[resultGUIrobust4D, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUIrobust,totalPhaseMatrix); 
+resultGUIrobust4D = matRad_calc4dDose(dij, pln, stf,resultGUIrobust,totalPhaseMatrix); 
+resultGUIrobust4D = matRad_acc4dDose( dij, pln, ct, cst,resultGUIrobust4D, 'DDM');
 
 %% Visualize results
 

examples/matRad_example11_helium.m

@@ -55,8 +55,6 @@
 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.propOpt.runDAO        = 0;
-pln.propOpt.runSequencing = 0;
 
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]

examples/matRad_example13_fitAnalyticalParticleBaseData.m

@@ -151,8 +151,6 @@
 pln.propOpt.optimizer       = 'IPOPT';
 pln.propOpt.bioOptimization = 'none'; % none: physical optimization;             const_RBExDose; constant RBE of 1.1;
                                       % LEMIV_effect: effect-based optimization; LEMIV_RBExDose: 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
 
 % retrieve scenarios for dose calculation and optimziation
 pln.multScen = matRad_multScen(ct,'nomScen');

examples/matRad_example14_spotRemoval.m

@@ -54,8 +54,6 @@
 pln.propStf.bixelWidth    = 3;
 pln.propStf.numOfBeams    = numel(pln.propStf.gantryAngles);
 pln.propStf.isoCenter     = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
-pln.propOpt.runDAO        = 0;
-pln.propOpt.runSequencing = 0;
 
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]

examples/matRad_example16_photonMC_MLC.m

@@ -42,9 +42,6 @@
 pln.propStf.bixelWidth      = 10;
 pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
 pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
-% Enable sequencing and direct aperture optimization (DAO).
-pln.propOpt.runSequencing   = 1;
-pln.propOpt.runDAO          = 1;
 
 quantityOpt    = 'physicalDose';                                     
 modelName      = 'none';  
@@ -74,11 +71,11 @@
 % 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,1);
-[pln,stf] = matRad_aperture2collimation(pln,stf,resultGUI.sequencing,resultGUI.apertureInfo);
+resultGUI = matRad_sequencing(resultGUI,stf,pln, dij);
+
 %% Aperture visualization
 % Use a matrad function to visualize the resulting aperture shapes
-matRad_visApertureInfo(resultGUI.apertureInfo)
+matRad_visApertureInfo(resultGUI.sequencing.apertureInfo)
 %% Plot the Resulting Dose Slice
 % Just let's plot the transversal iso-center dose slice
 slice = matRad_world2cubeIndex(pln.propStf.isoCenter(1,:),ct);

examples/matRad_example17_biologicalModels.m

@@ -15,8 +15,6 @@
 pln.propStf.numOfBeams    = numel(pln.propStf.gantryAngles);
 
 pln.propStf.isoCenter     = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
-pln.propOpt.runDAO        = 0;
-pln.propSeq.runSequencing = 0;
 
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 8;

examples/matRad_example18_FREDMC.m

@@ -34,8 +34,6 @@
 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.propOpt.runDAO        = 0;
-pln.propSeq.runSequencing = 0;
 
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]

examples/matRad_example19_CT_sCT_DVH_difference_photons.m

@@ -215,12 +215,6 @@
 pln.propDoseCalc.doseGrid.resolution.y = 7; % [mm]
 pln.propDoseCalc.doseGrid.resolution.z = 7; % [mm]
 
-%%
-% Enable sequencing and disable direct aperture optimization (DAO) for now.
-% A DAO optimization is shown in a seperate example.
-pln.propSeq.runSequencing = 1;
-pln.propOpt.runDAO        = 0;
-
 %% 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.

examples/matRad_example1_phantom.m

@@ -110,10 +110,6 @@
 pln.propStf.numOfBeams    = numel(pln.propStf.gantryAngles);
 pln.propStf.isoCenter     = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
 
-% Settings for Optimization
-pln.propOpt.runDAO        = 0;
-pln.propOpt.runSequencing = 0;
-
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
 pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]

examples/matRad_example2_photons.m

@@ -132,12 +132,6 @@
 pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
 pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
 
-%%
-% Enable sequencing and disable direct aperture optimization (DAO) for now.
-% A DAO optimization is shown in a seperate example.
-pln.propSeq.runSequencing = 1;
-pln.propOpt.runDAO        = 0;
-
 
 %%
 % and et voila our treatment plan structure is ready. Lets have a look:

examples/matRad_example3_photonsDAO.m

@@ -2,105 +2,110 @@
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-% 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/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 
+% 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/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.
 %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-%% In this example we will show 
+%% In this example we will show
 % (i) how to load patient data into matRad
-% (ii) how to setup a photon dose calculation and 
+% (ii) how to setup a photon dose calculation and
 % (iii) how to inversely optimize directly from command window in MatLab.
 % (iv) how to apply a sequencing algorithm
 % (v) how to run a direct aperture optimization
 % (iv) how to visually and quantitatively evaluate the result
 
 %% set matRad runtime configuration
-matRad_rc; %If this throws an error, run it from the parent directory first to set the paths
+matRad_rc; % If this throws an error, run it from the parent directory first to set the paths
 
 %% Patient Data Import
 % import the head & neck patient into your workspace.
 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 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;
- 
-pln.propStf.gantryAngles    = [0:72:359];
-pln.propStf.couchAngles     = [0 0 0 0 0];
+
+pln.propStf.gantryAngles    = [0:90:359];
+pln.propStf.couchAngles     = [0 0 0 0];
 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.propStf.isoCenter       = ones(pln.propStf.numOfBeams, 1) * matRad_getIsoCenter(cst, ct, 0);
 
-pln.bioModel = 'none'; 
+pln.bioModel = 'none';
 pln.multScen = 'nomScen';
 
 % dose calculation settings
-pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
-pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
-pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]
+pln.propDoseCalc.doseGrid.resolution.x = 5; % [mm]
+pln.propDoseCalc.doseGrid.resolution.y = 5; % [mm]
+pln.propDoseCalc.doseGrid.resolution.z = 5; % [mm]
 
-% We can also use other solver for optimization than IPOPT. matRad 
+% We can also use other solver for optimization than IPOPT. matRad
 % currently supports fmincon from the MATLAB Optimization Toolbox. First we
 % check if the fmincon-Solver is available, and if it es, we set in in the
 % pln.propOpt.optimizer vairable. Otherwise wie set to the default
 % optimizer 'IPOPT'
 if matRad_OptimizerFmincon.IsAvailable()
-    pln.propOpt.optimizer = 'fmincon';   
+    pln.propOpt.optimizer = 'fmincon';
 else
     pln.propOpt.optimizer = 'IPOPT';
 end
-pln.propOpt.quantityOpt = 'physicalDose';  
-
-%%
-% Enable sequencing and direct aperture optimization (DAO).
-pln.propSeq.runSequencing = true;
-pln.propOpt.runDAO        = true;
+pln.propOpt.quantityOpt = 'physicalDose';
 
 %% Generate Beam Geometry STF
-stf = matRad_generateStf(ct,cst,pln);
+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 
+% 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_calcDoseInfluence(ct,cst,stf,pln);
+dij = matRad_calcDoseInfluence(ct, cst, stf, pln);
 
 %% 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 
+% 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. Once the optimization has finished, trigger once the GUI to
 % visualize the optimized dose cubes.
-resultGUI = matRad_fluenceOptimization(dij,cst,pln);
-matRadGUI;
+resultGUI = matRad_fluenceOptimization(dij, cst, pln);
+% 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 
+% 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.
-resultGUI = matRad_sequencing(resultGUI,stf,dij,pln);
+
+%% some testing of sequencing
+pln.propSeq.sequencer = 'siochi';
+pln.propSeq.sequencingLevel = 10;
+resultGUI_SIOCHI = matRad_sequencing(resultGUI, stf, pln, dij);
+
+pln.propSeq.sequencer = 'xia';
+resultGUI_XIA = matRad_sequencing(resultGUI, stf, pln, dij);
+
+pln.propSeq.sequencer = 'engel';
+resultGUI_ENGEL = matRad_sequencing(resultGUI, stf, pln, dij);
 
 %% DAO - Direct Aperture Optimization
-% The Direct Aperture Optimization is an optimization approach where we 
+% The Direct Aperture Optimization is an optimization approach where we
 % directly optimize aperture shapes and weights.
-resultGUI = matRad_directApertureOptimization(dij,cst,resultGUI.apertureInfo,resultGUI,pln);
+resultGUI_SIOCHI_DAO = matRad_directApertureOptimization(dij, cst, resultGUI_SIOCHI.sequencing.apertureInfo, pln);
 
 %% Aperture visualization
 % Use a matrad function to visualize the resulting aperture shapes
-matRad_visApertureInfo(resultGUI.apertureInfo);
+matRad_visApertureInfo(resultGUI_SIOCHI_DAO.sequencing.apertureInfo);
 
 %% Indicator Calculation and display of DVH and QI
-resultGUI = matRad_planAnalysis(resultGUI,ct,cst,stf,pln);
+resultGUI = matRad_planAnalysis(resultGUI, ct, cst, stf, pln);

examples/matRad_example4_photonsMC.m

@@ -43,8 +43,6 @@
 pln.propStf.bixelWidth      = 10;
 pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
 pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
-pln.propSeq.runSequencing   = 0;
-pln.propOpt.runDAO          = 0;
 
 % dose calculation settings
 % We can choose a different dose calculation engine, here "ompMC", by

examples/matRad_example5_protons.m

@@ -63,8 +63,6 @@
 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.propOpt.runDAO        = 0;
-pln.propSeq.runSequencing = 0;
 
 % dose calculation settings
 pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]

examples/matRad_example9_4DDoseCalcMinimal.m

@@ -71,7 +71,8 @@
 %%
 % calc 4D dose
 % make sure that the correct pln, dij and stf are loeaded in the workspace
-[resultGUI, timeSequence] = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI); 
+resultGUI = matRad_calc4dDose(dij, pln,stf,resultGUI);
+resultGUI = matRad_acc4dDose( dij, pln, ct, cst,resultGUI, 'DDM');
 
 % plot the result in comparison to the static dose
 slice = matRad_world2cubeIndex(pln.propStf.isoCenter(1,:),ct);

matRad/4D/matRad_acc4dDose.m

@@ -0,0 +1,76 @@
+function resultGUI = matRad_acc4dDose( dij, pln, ct, cst,resultGUI, accType)
+% wrapper for the whole 4D dose calculation pipeline and calculated dose
+% accumulation
+%
+% call
+%   ct = matRad_calc4dDose(ct, pln, dij, stf, cst, resultGUI)
+%
+% input
+%   ct :             ct cube
+%   pln:             matRad plan meta information struct
+%   dij:             matRad dij struct
+%   stf:             matRad steering information struct
+%   cst:             matRad cst struct
+%   resultGUI:       struct containing optimized fluence vector
+%   totalPhaseMatrix optional intput for totalPhaseMatrix
+%   accType:         witch algorithim for dose accumulation
+% output
+%   resultGUI:      structure containing phase dose, RBE weighted dose, etc
+%   timeSequence:   timing information about the irradiation
+%
+% References
+%   -
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright 2018 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();
+
+if ~isa(pln.bioModel,'matRad_BiologicalModel')
+    pln.bioModel = matRad_BiologicalModel.validate(pln.bioModel,pln.radiationMode);
+end
+
+% accumulation
+resultGUI.accPhysicalDose = matRad_doseAcc(ct,resultGUI.phaseDose, cst, accType);
+if isa(pln.bioModel,'matRad_ConstantRBE')
+
+    resultGUI.accRBExDose = matRad_doseAcc(ct,resultGUI.phaseRBExDose, cst, accType);
+
+elseif isa(pln.bioModel,'matRad_LQBasedModel')
+
+    resultGUI.accAlphaDose    = matRad_doseAcc(ct,resultGUI.phaseAlphaDose, cst,accType);
+    resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType);
+
+    % only compute where we have biologically defined tissue
+    ix = (ax{1} ~= 0);
+
+    resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2;
+
+    resultGUI.accRBExDose     = zeros(ct.cubeDim);
+    resultGUI.accRBExDose(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix)));
+end
+
+for beamIx = 1:dij.numOfBeams
+    resultGUI.(['accPhysicalDose_beam', num2str(beamIx)])= matRad_doseAcc(ct,resultGUI.(['phaseDose_beam', num2str(beamIx)]), cst, accType);
+    if isa(pln.bioModel,'matRad_ConstantRBE')
+        resultGUI.(['accRBExDose_beam', num2str(beamIx)]) = matRad_doseAcc(ct,resultGUI.(['phaseRBExDose_beam', num2str(beamIx)]), cst, accType);
+   elseif isa(pln.bioModel,'matRad_LQBasedModel')
+        resultGUI.(['accAlphaDose_beam', num2str(beamIx)])      = matRad_doseAcc(ct,resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]), cst, accType);
+        resultGUI.(['accSqrtBetaDose_beam', num2str(beamIx)])   = matRad_doseAcc(ct,resultGUI.(['phaseAlphaDose_beam', num2str(beamIx)]), cst, accType);
+        resultGUI.(['accEffect_beam', num2str(beamIx)])         = resultGUI.(['accAlphaDose_beam', num2str(beamIx)]) + resultGUI.(['accSqrtBetaDose_beam', num2str(beamIx)]).^2;
+        resultGUI.(['accRBExDose_beam', num2str(beamIx)]){i}       = zeros(ct.cubeDim);
+        resultGUI.(['accRBExDose_beam', num2str(beamIx)]){i}       =  ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.(['accEffect_beam', num2str(beamIx)]){i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix)));    
+    end
+end
+
+
+end