Class based scenarios and simulations

MatRadMultScen.m

@@ -0,0 +1,135 @@
+classdef MatRadMultScen < handle 
+  properties               
+
+      % scenarios
+      numOfCtScen;            % number of imported ct scenarios                        
+      % shift scenarios
+      shiftSize;              % 3x1 vector to define maximal shift in [mm]  % (e.g. abdominal cases 5mm otherwise 3mm)
+      maxAbsRangeShift;       % maximum absolute over and undershoot in mm   
+      maxRelRangeShift;       % maximum relative over and undershoot in % 
+
+
+      % combination
+      rangeCombType;          % 'permuted':    create every possible range combination
+                              % 'combined':    combine absolute and relative range scenario
+      shiftCombType;          % 'permuted':    create every possible shift combination; 
+      scenCombType;           % 'permuted':    create every possible combination of range and setup scenarios
+
+      samplingMethod;           % 'random' / 'grid'
+      includeNomScen;         % boolean to determine if the nominal scenario should be included
+
+      totNumScen;             % total number of scenarios (might be changed later)
+      numOfFrac;              % number of fractions
+
+      shiftSD_sys     = [0 0 0];   % given in [mm]
+      shiftSD_rand    = [2 2 2];   % given in [mm]
+
+      rangeRelSD_sys  = 1.8;       % given in %
+      rangeRelSD_rand = 0;         % given in %
+      rangeAbsSD_sys  = 0.8;       % given in [mm]
+      rangeAbsSD_rand = 0;         % given in [mm]
+
+  end
+
+  % private properties which can only be changed inside matRad_multScen
+  properties (SetAccess = private)
+      fraction_index % determines which dose j belongs to which treatment scenario i
+      scenForProb;            % matrix for probability calculation - each row denotes one scenario  
+      % these parameters will be filled according to the choosen scenario type
+      isoShift;
+      relRangeShift;
+      absRangeShift;
+      scenMask;
+
+      DEFAULT_TYPE      = 'nomScen'; 
+      DEFAULT_probDist  = 'normDist'; % 'normDist': normal probability distrubtion or 'equalProb' for uniform probability distribution  
+      DEFAULT_TypeProp  = 'probBins'; % 'probBins': cumulative prob in bin or 'pointwise' for point probability  
+  end
+
+  % constant properties which are visible outside of matRad_multScen
+  properties (Constant = true)
+      AvailableScenCreationTYPE = {'nomScen','wcScen','impScen','rndScen'};
+  end
+
+  properties (Dependent = true)
+      SD_sys_all;
+      SD_rand_all;
+      scenProb;               % probability of each scenario stored in a vector
+      numberOfTreatmentScen;
+  end
+
+  methods
+    function this = MatRadMultScen()
+
+    end
+
+    function randomSampling(this, totNumScen, numOfFrac, shiftSD_sys, shiftSD_rand, rangeRelSD_sys, rangeRelSD_rand, rangeAbsSD_sys, rangeAbsSD_rand)
+        this.totNumScen = totNumScen;
+        this.numOfFrac = numOfFrac;
+        this.shiftSD_sys = shiftSD_sys;
+        this.shiftSD_rand = shiftSD_rand;
+        this.rangeRelSD_sys = rangeRelSD_sys;
+        this.rangeRelSD_rand = rangeRelSD_rand;
+        this.rangeAbsSD_sys = rangeAbsSD_sys;
+        this.rangeAbsSD_rand = rangeAbsSD_rand;
+
+        this.createRandomSampledScenarios();
+    end
+
+    function w = get.scenProb(this)
+        w = this.computeScenProb();
+    end
+
+    function numberOfTreatmentScen = get.numberOfTreatmentScen(this)
+        numberOfTreatmentScen = ceil(this.totNumScen / this.numOfFrac);
+    end
+    function SD_sys_all = get.SD_sys_all(this)
+        SD_sys_all = [this.shiftSD_sys this.rangeAbsSD_sys this.rangeRelSD_sys];
+    end
+
+    function SD_rand_all = get.SD_rand_all(this)
+        SD_rand_all = [this.shiftSD_rand this.rangeAbsSD_rand this.rangeRelSD_rand];
+    end
+
+  end % eof public methods
+
+  methods (Access = private)
+
+    function createRandomSampledScenarios(this)
+        this.totNumScen = this.numberOfTreatmentScen * this.numOfFrac;
+
+        this.scenForProb = NaN * ones(this.totNumScen, numel(this.SD_sys_all));
+        this.fraction_index = NaN * ones(this.totNumScen, 1);
+
+        runIx = 1;
+        for i = 1:this.numberOfTreatmentScen
+          mu = this.sampleFromNormalDist(0, this.SD_sys_all, size(this.SD_sys_all));
+          for j = 1:this.numOfFrac
+            this.scenForProb(runIx, :) = this.sampleFromNormalDist(mu, this.SD_rand_all, size(this.SD_rand_all));
+            this.fraction_index(runIx) = i;
+            runIx = runIx + 1;
+          end
+        end
+        this.isoShift = this.scenForProb(:,1:3);
+        this.absRangeShift = this.scenForProb(:,4);
+        this.relRangeShift = this.scenForProb(:,5);
+        this.scenMask = ones(this.numOfCtScen, size(this.isoShift, 1), numel(this.absRangeShift));
+    end
+
+    function scenProb = computeScenProb(this)
+        if strcmp(this.samplingMethod, 'random')
+            scenProb = 1 / this.totNumScen * ones(this.totNumScen, 1);
+        else
+            error('Do that later');
+        end
+    end
+
+  end % eof private methods
+
+  methods (Static)
+    function v = sampleFromNormalDist(mu,sigma, dim)
+      v = sigma .* randn(dim(1),dim(2)) + mu;
+    end
+  end
+
+end

matRad.m

@@ -97,6 +97,3 @@
 structSel = {}; % structSel = {'PTV','OAR1'};
 [caSampRes, mSampDose, plnSamp, resultGUInomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,[],[]);
 [cstStat, resultGUIStat, param]                   = matRad_samplingAnalysis(ct,cst,plnSamp,caSampRes, mSampDose, resultGUInomScen,[]);
-
-
-

matRad_calcDoseDirect.m

@@ -50,7 +50,7 @@
 end
 
 % copy bixel weight vector into stf struct
-if exist('w','var')
+if exist('w','var') && ~isempty(w)
     if sum([stf.totalNumOfBixels]) ~= numel(w)
         error('weighting does not match steering information')
     end

matRad_calcParticleDose.m

@@ -62,7 +62,7 @@
 ct = matRad_calcWaterEqD(ct, pln, param);
 
 % meta information for dij
-dij.numOfBeams         = pln.propStf.numOfBeams;
+dij.numOfBeams         = numel(stf);
 dij.numOfVoxels        = prod(ct.cubeDim);
 dij.resolution         = ct.resolution;
 dij.dimensions         = ct.cubeDim;
@@ -271,7 +271,7 @@
     % compute SSDs
     stf = matRad_computeSSD(stf,ct,ctScen,param);
 
-   for i = 1:dij.numOfBeams % loop over all beams
+   for i = 1:numel(stf) % loop over all beams
 
        matRad_dispToConsole(['Beam ' num2str(i) ' of ' num2str(dij.numOfBeams) ':  \n'],param,'info');
 
@@ -295,7 +295,7 @@
        % transformation of the coordinate system need double transpose
 
        % rotation around Z axis (gantry)
-       rotMat_system_T = matRad_getRotationMatrix(pln.propStf.gantryAngles(i),pln.propStf.couchAngles(i));
+       rotMat_system_T = matRad_getRotationMatrix(stf(i).gantryAngle,stf(i).couchAngle);
 
        % Rotate coordinates (1st couch around Y axis, 2nd gantry movement)
        rot_coordsV = coordsV*rotMat_system_T;

matRad_multScen.m

@@ -56,6 +56,8 @@
         shiftCombType;          % 'individual':  no combination of shift scenarios;       number of shift scenarios is sum(multScen.numOfShiftScen) 
                                 % 'combined':    combine shift scenarios;                 number of shift scenarios is numOfShiftScen
                                 % 'permuted':    create every possible shift combination; number of shift scenarios is 8,27,64 .
+
+        isoShiftCorrelationBreak = 'fraction';    % correlation type for isoshifts ('fraction', 'beam', 'ray')
 
         % b) define range error scenarios                                                
         numOfRangeShiftScen;    % number of absolute and/or relative range scnearios. 
@@ -65,6 +67,8 @@
         rangeCombType;          % 'individual':  no combination of absolute and relative range scenarios
                                 % 'combined':    combine absolute and relative range scenarios
         rangeGenType;           % 'equidistant': equidistant range shifts, 'sampled': sample range shifts from normal distribution
+
+        rangeShiftCorrelationBreak = 'fraction';  % correlation type for rangeshifts ('fraction', 'beam', 'ray')
         scenCombType;           % 'individual':  no combination of range and setup scenarios, 
                                 % 'combined':    combine range and setup scenarios if their scenario number is consistent 
                                 % 'permuted':    create every possible combination of range and setup scenarios
@@ -89,6 +93,8 @@
         rangeAbsSD  = 1;                  % given in [mm]   
         shiftSD     = [2.25 2.25 2.25];   % given in [mm]
 
+        deepestCorrelationBreak = 'fraction';    % argmax of isoShiftCorrelationBreak, rangeShiftCorrelationBreak
+
     end
 
     % private properties which can only be changed inside matRad_multScen
@@ -104,6 +110,8 @@
     properties(Constant = true)
 
         AvailableScenCreationTYPE = {'nomScen','wcScen','impScen','rndScen'};
+        AvailableCorrelationTYPES = {'fraction', 'beam', 'ray'};
+        AvailableUncertainties = {'range', 'setup'};
     end
 
     % constant private properties which are only visible within matRad_multScen
@@ -198,6 +206,32 @@ function attachListener(obj)
           this      = calcScenProb(this);
       end
 
+      function this = matRad_recreateInstance(this, keep)
+          % function to recreate instance (usefule if sampling is used, and
+          % one wants to sample again. Possibility to keep one, either or
+          % both uncertainties fixed.
+          if ~exist('keep', 'var') || isempty(keep)
+              this = setMultScen(this);
+
+          elseif isequal(sort(keep), sort(this.AvailableUncertainties))
+              % do nothing
+
+          elseif strcmp(keep, 'range')
+               rangeShiftRel = this.relRangeShift;
+               rangeShiftAbs = this.absRangeShift;
+               this = setMultScen(this);
+
+               this.relRangeShift = rangeShiftRel;
+               this.absRangeShift = rangeShiftAbs;
+
+          elseif strcmp(keep, 'setup')
+               isoShift = this.isoShift;
+               this = setMultScen(this);
+               this.isoShift = isoShift;
+          end
+          this = calcScenProb(this);
+      end
+
       %% setters to check for valid input
       function this = set.numOfShiftScen(this,value)
 
@@ -472,7 +506,19 @@ function attachListener(obj)
                        end
                end
             end
-
+
+            % check correlationt types
+            if ~ismember(this.isoShiftCorrelationBreak, this.AvailableCorrelationTYPES) ...
+                    || ~ismember(this.rangeShiftCorrelationBreak, this.AvailableCorrelationTYPES)
+                matRad_dispToConsole('Correlationbreak not implemented as specified. \n',[],'error');
+            end
+            if strcmp(this.isoShiftCorrelationBreak, 'ray') || strcmp(this.rangeShiftCorrelationBreak, 'ray')
+                this.deepestCorrelationBreak = 'ray';
+            elseif strcmp(this.isoShiftCorrelationBreak, 'beam') || strcmp(this.rangeShiftCorrelationBreak, 'beam')
+                 this.deepestCorrelationBreak = 'beam';
+            else
+                this.deepestCorrelationBreak = 'fraction';
+            end
 
             % create linearalized mask where the i row points to the indexes of scenMask
             [x{1}, x{2}, x{3}] = ind2sub(size(this.scenMask),find(this.scenMask));
@@ -538,14 +584,3 @@ function attachListener(obj)
 
 
 end  % end of matRad_multScen class
-
-
-
-
-
-
-
-
-
-
-

matRad_sampling.m

@@ -1,4 +1,4 @@
-function [caSampRes, mSampDose, pln, resultGUInomScen]  = matRad_sampling(ct,stf,cst,pln,w,structSel,multScen,param)
+function [treatmentSimulation, scenContainer, pln, resultGUInomScen, nomScen]  = matRad_sampling(ct,stf,cst,pln,w,structSel,multScen,param)
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % matRad_randomSampling enables sampling multiple treatment scenarios
 % 
@@ -77,7 +77,7 @@
     for i=1:size(cst,1)
         for j = 1:numel(structSel)
             if strcmp(structSel{j}, cst{i,2})
-                V = [V cst{i,4}{1}];
+                V = [V {cst{i,4}{1}}];
             end
         end
     end
@@ -125,10 +125,16 @@
 resultGUInomScen.qi  = nomQi;
 resultGUInomScen.cst = cst;
 
-% only show errors and disable waitbars and figures
-param.logLevel = 4;
+nomDose                = resultGUInomScen.(pln.bioParam.quantityVis)(param.subIx);
+nomDoseLin             = single(reshape(nomDose,[],1));
+
+nomScen = MatRadScenario(nomDoseLin, param.subIx, pln.bioParam.quantityVis, 1, [0 0 0], 0, 0, size(resultGUInomScen.(pln.bioParam.quantityVis)), NaN);
+
+% only show warnings and disable waitbars and figures
+param.logLevel = 3;
 
 %% perform parallel sampling
+scenContainer = cell(pln.multScen.totNumScen,1);
 if FlagParallToolBoxLicensed
    % Create parallel pool on cluster
    p = gcp(); % If no pool, create new one.
@@ -150,7 +156,7 @@
       fprintf('matRad: Consider downloading parfor_progress function from the matlab central fileexchange to get feedback from parfor loop.\n');
       FlagParforProgressDisp = false;
    end
-
+   
    parfor i = 1:pln.multScen.totNumScen
 
           % create nominal scenario
@@ -170,14 +176,12 @@
 
           resultSamp                 = matRad_calcDoseDirect(ct,stf,plnSamp,cst,w,param);
           sampledDose                = resultSamp.(pln.bioParam.quantityVis)(param.subIx);
-          mSampDose(:,i)             = single(reshape(sampledDose,[],1));
-          caSampRes(i).bioParam      = pln.bioParam;
-          caSampRes(i).relRangeShift = plnSamp.multScen.relRangeShift;
-          caSampRes(i).absRangeShift = plnSamp.multScen.absRangeShift;
-          caSampRes(i).isoShift      = plnSamp.multScen.isoShift;
+          sampledDoseLin             = single(reshape(sampledDose,[],1));
+          scenContainer{i}           = MatRadScenario(sampledDoseLin, param.subIx, pln.bioParam.quantityVis, 1, plnSamp.multScen.isoShift, ...
+                plnSamp.multScen.relRangeShift, plnSamp.multScen.absRangeShift, size(resultSamp.(pln.bioParam.quantityVis)), 1);
 
-          caSampRes(i).dvh = matRad_calcDVH(cst,resultSamp.(pln.bioParam.quantityVis),'cum',dvhPoints);
-          caSampRes(i).qi  = matRad_calcQualityIndicators(cst,pln,resultSamp.(pln.bioParam.quantityVis),refGy,refVol,param);
+          % caSampRes(i).dvh = matRad_calcDVH(cst,resultSamp.(pln.bioParam.quantityOpt),'cum',dvhPoints);
+          % caSampRes(i).qi  = matRad_calcQualityIndicators(cst,pln,resultSamp.(pln.bioParam.quantityOpt),refGy,refVol,param);
 
           if FlagParforProgressDisp
             parfor_progress;
@@ -214,14 +218,13 @@
 
           resultSamp                 = matRad_calcDoseDirect(ct,stf,plnSamp,cst,w,param);
           sampledDose                = resultSamp.(pln.bioParam.quantityVis)(param.subIx);
-          mSampDose(:,i)             = single(reshape(sampledDose,[],1));
-          caSampRes(i).bioParam      = pln.bioParam;
-          caSampRes(i).relRangeShift = plnSamp.multScen.relRangeShift;
-          caSampRes(i).absRangeShift = plnSamp.multScen.absRangeShift;
-          caSampRes(i).isoShift      = plnSamp.multScen.isoShift;
+          sampledDoseLin             = single(reshape(sampledDose,[],1));
+          scenContainer{i}           = MatRadScenario(sampledDoseLin, param.subIx, pln.bioParam.quantityVis, 1, plnSamp.multScen.isoShift, ...
+                plnSamp.multScen.relRangeShift, plnSamp.multScen.absRangeShift, size(resultSamp.(pln.bioParam.quantityVis)), 1);
 
-          caSampRes(i).dvh = matRad_calcDVH(cst,resultSamp.(pln.bioParam.quantityVis),'cum',dvhPoints);
-          caSampRes(i).qi  = matRad_calcQualityIndicators(cst,pln,resultSamp.(pln.bioParam.quantityVis),refGy,refVol,param);
+          % caSampRes(i).dvh = matRad_calcDVH(cst,resultSamp.(pln.bioParam.quantityOpt),'cum',dvhPoints);
+          % caSampRes(i).qi  = matRad_calcQualityIndicators(cst,pln,resultSamp.(pln.bioParam.quantityOpt),refGy,refVol,param);
+
           matRad_progress(i, pln.multScen.totNumScen)
     end
 
@@ -230,4 +233,9 @@
 %% add subindices
 pln.subIx        = param.subIx;
 
+%% account for fractionation
+treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nomScen, ct, cst, pln, pln.numOfFractions);
+for i = 1:numel(scenContainer)
+    treatmentSimulation.initNewScen(scenContainer{i})
+end
 end