diff --git a/scenarios/matRad_GriddedScenariosAbstract.m b/scenarios/matRad_GriddedScenariosAbstract.m new file mode 100644 index 000000000..a52df1216 --- /dev/null +++ b/scenarios/matRad_GriddedScenariosAbstract.m @@ -0,0 +1,262 @@ +classdef (Abstract) matRad_GriddedScenariosAbstract < matRad_ScenarioModel + %UNTITLED Summary of this class goes here + % Detailed explanation goes here + + properties (AbortSet = true) + includeNominalScenario = true; + combinations = 'none'; %Can be 'none', 'shift', 'all' to ontrol creation of worst case combinations + combineRange = true; %Wether to treat absolute & relative range as one shift or as separate scenarios + end + + %Each subclass needs to define how many gridpoints it uses and if this + %can be set or not + properties (Abstract) + numOfSetupGridPoints; + numOfRangeGridPoints; + end + + properties (Constant) + validCombinationTypes = {'all','none','shift'}; + end + + methods + function this = matRad_GriddedScenariosAbstract(ct) + + if nargin == 0 + superclassArgs = {}; + else + superclassArgs = {ct}; + end + + this@matRad_ScenarioModel(superclassArgs{:}); + + %TODO: We could do this automatically in the superclass + %Octave 5 has a bug there and throws an error + %this.updateScenarios(); + end + + function set.combineRange(this,combineRange_) + valid = isscalar(combineRange_) && (isnumeric(combineRange_) || islogical(combineRange_)); + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for combineRange! Needs to be a boolean / logical value!'); + end + this.combineRange = combineRange_; + this.updateScenarios(); + end + + %% set methods + function set.includeNominalScenario(this,includeNomScen) + valid = isscalar(includeNomScen) && (isnumeric(includeNomScen) || islogical(includeNomScen)); + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for includeNominalScenario! Needs to be a boolean / logical value!'); + end + this.includeNominalScenario = includeNomScen; + this.updateScenarios(); + end + + function set.combinations(this,combinations_) + valid = any(strcmp(combinations_,this.validCombinationTypes)); + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for combinations! Needs to be one of the strings %s!',strjoin(this.validCombinationTypes,' / ')); + end + this.combinations = combinations_; + this.updateScenarios(); + end + + function scenarios = updateScenarios(this) + matRad_cfg = MatRad_Config.instance(); + + %Get the maximum, i.e., worst case shifts + wcSetupShifts = this.wcSigma * this.shiftSD; + + %% Create gridded setup shifts + %Create grid vectors for setup shifts + setupShiftGrid = zeros(this.numOfSetupGridPoints,numel(wcSetupShifts)); + if mod(this.numOfSetupGridPoints,2) == 0 && this.includeNominalScenario + matRad_cfg.dispWarning('Obtaining Setup Shifts: Including the nominal scenario with even number of grid points creates asymmetrical shifts!'); + end + + for i = 1:numel(wcSetupShifts) + setupShiftGrid(:,i) = linspace(-wcSetupShifts(i),wcSetupShifts(i),this.numOfSetupGridPoints); + if this.includeNominalScenario + + [~,ix] = min(abs(setupShiftGrid(:,i))); + setupShiftGrid(ix,i) = 0; + end + end + + %Now create vector of all shifts for different combinatorial + %settings + switch this.combinations + case 'none' + %Create independent shifts + griddedSetupShifts = []; + for i=1:size(setupShiftGrid,2) + tmpGrid = zeros(size(setupShiftGrid,1),3); + tmpGrid(:,i) = setupShiftGrid(:,i); + griddedSetupShifts = [griddedSetupShifts; tmpGrid]; + end + case {'shift','all'} + [X,Y,Z] = meshgrid(setupShiftGrid(:,1),setupShiftGrid(:,2),setupShiftGrid(:,3)); + griddedSetupShifts = [X(:), Y(:), Z(:)]; + otherwise + matRad_cfg.dispError('Invalid value for combinations! This sanity check should never be reached!'); + end + + griddedSetupShifts = matRad_ImportanceScenarios.uniqueStableRowsCompat(griddedSetupShifts); + shiftNomScenIx = find(all(griddedSetupShifts == zeros(1,3),2)); + + if ~isempty(shiftNomScenIx) || this.includeNominalScenario + if ~isempty(shiftNomScenIx) + griddedSetupShifts(shiftNomScenIx,:) = []; + end + griddedSetupShifts = [0 0 0; griddedSetupShifts]; + end + + this.totNumShiftScen = size(griddedSetupShifts,1); + + %% Create gridded range shifts + %Obtain worst case range shifts + wcRangeShifts = this.wcSigma * [this.rangeAbsSD this.rangeRelSD./100]; + + rangeShiftGrid = zeros(this.numOfRangeGridPoints,numel(wcRangeShifts)); + if mod(this.numOfRangeGridPoints,2) == 0 && this.includeNominalScenario + matRad_cfg.dispWarning('Obtaining Range Shifts: Including the nominal scenario with even number of grid points creates asymmetrical shifts!'); + end + + for i = 1:numel(wcRangeShifts) + rangeShiftGrid(:,i) = linspace(-wcRangeShifts(i),wcRangeShifts(i),this.numOfRangeGridPoints); + if this.includeNominalScenario + + [~,ix] = min(abs(rangeShiftGrid(:,i))); + rangeShiftGrid(ix,i) = 0; + end + end + + if this.combineRange + griddedRangeShifts = rangeShiftGrid; + else + [rngAbs,rngRel] = meshgrid(rangeShiftGrid(:,1),rangeShiftGrid(:,2)); + griddedRangeShifts = [rngAbs(:),rngRel(:)]; + end + + %Remove duplicate scenarios and update number of shifts + griddedRangeShifts = this.uniqueStableRowsCompat(griddedRangeShifts); + + rangeNomScenIx = find(all(griddedRangeShifts == zeros(1,2),2)); + + if ~isempty(rangeNomScenIx) || this.includeNominalScenario + if ~isempty(rangeNomScenIx) + griddedRangeShifts(rangeNomScenIx,:) = []; + end + griddedRangeShifts = [0 0; griddedRangeShifts]; + end + + this.totNumRangeScen = size(griddedRangeShifts,1); + + %Aggregate scenarios + switch this.combinations + case {'none','shift'} + scenarios = zeros(this.totNumShiftScen + this.totNumRangeScen,5); + scenarios(1:this.totNumShiftScen,1:3) = griddedSetupShifts; + scenarios(this.totNumShiftScen+1:this.totNumShiftScen+this.totNumRangeScen,4:5) = griddedRangeShifts; + + %create the linear mask of scenarios + linearMaskTmp = ones(size(scenarios,1),3); + linearMaskTmp(1:this.totNumShiftScen,2) = (1:this.totNumShiftScen)'; + linearMaskTmp(this.totNumShiftScen+1:this.totNumShiftScen+this.totNumRangeScen,3) = (1:this.totNumRangeScen)'; + + [scenarios,ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); + linearMaskTmp = linearMaskTmp(ia,:); + + case 'all' + %Prepare scenario matrix by replicating shifts + %with the number of range scenarios + scenarios = repmat(griddedSetupShifts,this.totNumRangeScen,1); + scenarios = [scenarios zeros(size(scenarios,1),2)]; + + %create the linear mask of scenarios + linearMaskTmp = ones(size(scenarios,1),3); + for r = 1:this.totNumRangeScen + offset = (r-1)*this.totNumShiftScen; + ixR = (offset + 1):(offset + this.totNumShiftScen); + scenarios(ixR,4:5) = repmat(griddedRangeShifts(r,:),this.totNumShiftScen,1); + + %Set linear mask + linearMaskTmp(ixR,2) = (1:this.totNumShiftScen)'; + linearMaskTmp(ixR,3) = r; + end + + + %create the linear mask of scenarios + [scenarios,ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); + linearMaskTmp = linearMaskTmp(ia,:); + otherwise + matRad_cfg.dispError('Invalid value for combinations! This sanity check should never be reached!'); + end + + %Handle 4D phases + phases = repmat(1:this.numOfCtScen,size(scenarios,1),1); + phases = phases(:); + scenarios = horzcat(phases, repmat(scenarios,[this.numOfCtScen 1])); + linearMaskTmp = repmat(linearMaskTmp,this.numOfCtScen,1); + linearMaskTmp(:,1) = phases; + + %Finalize meta information + this.totNumScen = size(scenarios,1); + + this.relRangeShift = scenarios(:,6); + this.absRangeShift = scenarios(:,5); + this.isoShift = scenarios(:,2:4); + + this.maxAbsRangeShift = max(this.absRangeShift); + this.maxRelRangeShift = max(this.relRangeShift); + + this.scenMask = false(this.numOfCtScen,this.totNumShiftScen,this.totNumRangeScen); + + this.scenForProb = scenarios; + this.linearMask = linearMaskTmp; + + maskIx = sub2ind(size(this.scenMask),linearMaskTmp(:,1),linearMaskTmp(:,2),linearMaskTmp(:,3)); + this.scenMask(maskIx) = true; + + %Get Scenario probability + %First, we use the Gaussian Uncertainty model for range and + %setup + Sigma = diag([this.shiftSD,this.rangeAbsSD,this.rangeRelSD./100].^2); + d = size(Sigma,1); + [cs,p] = chol(Sigma); + + tmpScenProb = (2*pi)^(-d/2) * exp(-0.5*sum((scenarios(:,2:end)/cs).^2, 2)) / prod(diag(cs)); + + %Now we combine with the 4D ct phase probabilities (multiply) + tmpPhaseProb = arrayfun(@(phase) this.phaseProbability(phase),phases); + + %Finalize probabilities + this.scenProb = tmpPhaseProb .* tmpScenProb; + this.scenWeight = this.scenProb./sum(this.scenProb); + + %TODO: Discard scenarios with probability 0? + end + end + + methods (Static) + function [uniqueStableRows,ia] = uniqueStableRowsCompat(values) + %This is a compatability wrapper to call unique without sorting + + matRad_cfg = MatRad_Config.instance(); + + if matRad_cfg.isOctave + %https://stackoverflow.com/questions/37828894/ + [~,ia,~] = unique(values,'rows','first'); + ia = sort(ia); + uniqueStableRows = values(ia,:); + else + [uniqueStableRows,ia] = unique(values,'rows','stable'); + end + end + end +end \ No newline at end of file diff --git a/scenarios/matRad_ImportanceScenarios.m b/scenarios/matRad_ImportanceScenarios.m index 3d3b286de..ea7f2718e 100644 --- a/scenarios/matRad_ImportanceScenarios.m +++ b/scenarios/matRad_ImportanceScenarios.m @@ -1,4 +1,4 @@ -classdef matRad_ImportanceScenarios < matRad_ScenarioModel +classdef matRad_ImportanceScenarios < matRad_GriddedScenariosAbstract % matRad_ImportanceScenarios % Implements gridded importance scenarios, i.e., weighted according to a % probability distribution. It is not advised to create "combined" grids @@ -28,10 +28,6 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% properties (AbortSet = true) - includeNominalScenario = true; - combinations = 'none'; %Can be 'none', 'shift', 'all' to ontrol creation of worst case combinations - combineRange = true; %Wether to treat absolute & relative range as one shift or as separate scenarios - numOfSetupGridPoints = 9; numOfRangeGridPoints = 9; end @@ -39,10 +35,6 @@ properties (SetAccess=protected) name = 'impScen'; end - - properties (Constant) - validCombinationTypes = {'all','none','shift'}; - end methods function this = matRad_ImportanceScenarios(ct) @@ -52,227 +44,32 @@ superclassArgs = {ct}; end - this@matRad_ScenarioModel(superclassArgs{:}); + this@matRad_GriddedScenariosAbstract(superclassArgs{:}); %TODO: We could do this automatically in the superclass %Octave 5 has a bug there and throws an error this.updateScenarios(); end - function set.combineRange(this,combineRange_) - valid = isscalar(combineRange_) && (isnumeric(combineRange_) || islogical(combineRange_)); - if ~valid - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for combineRange! Needs to be a boolean / logical value!'); - end - this.combineRange = combineRange_; - this.updateScenarios(); - end - - function scenarios = updateScenarios(this) - - [scenarios,linearMask,this.totNumShiftScen,this.totNumRangeScen] = this.createGriddedScenarios(this,this.includeNominalScenario,this.combinations,this.combineRange,this.numOfSetupGridPoints,this.numOfRangeGridPoints); - - %Finalize meta information - this.totNumScen = size(scenarios,1); - - - this.relRangeShift = scenarios(:,5); - this.absRangeShift = scenarios(:,4); - this.isoShift = scenarios(:,1:3); - - this.maxAbsRangeShift = max(this.absRangeShift); - this.maxRelRangeShift = max(this.relRangeShift); - - this.scenMask = false(this.numOfCtScen,this.totNumShiftScen,this.totNumRangeScen); - - this.scenForProb = scenarios; - this.linearMask = linearMask; - - maskIx = sub2ind(size(this.scenMask),linearMask(:,1),linearMask(:,2),linearMask(:,3)); - this.scenMask(maskIx) = true; - - %Get Scenario probability - Sigma = diag([this.shiftSD,this.rangeAbsSD,this.rangeRelSD./100].^2); - d = size(Sigma,1); - [cs,p] = chol(Sigma); - this.scenProb = (2*pi)^(-d/2) * exp(-0.5*sum((scenarios/cs).^2, 2)) / prod(diag(cs)); - this.scenWeight = this.scenProb./sum(this.scenProb); - end - - %% set methods - function set.includeNominalScenario(this,includeNomScen) - valid = isscalar(includeNomScen) && (isnumeric(includeNomScen) || islogical(includeNomScen)); + function set.numOfSetupGridPoints(this,numGridPoints) + valid = isscalar(numGridPoints) && numGridPoints > 0; if ~valid matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for includeNominalScenario! Needs to be a boolean / logical value!'); + matRad_cfg.dispError('Invalid number of setup grid points, needs to be a positive scalar!'); end - this.includeNominalScenario = includeNomScen; + this.numOfSetupGridPoints = inumGridPoints; this.updateScenarios(); end - function set.combinations(this,combinations_) - valid = any(strcmp(combinations_,this.validCombinationTypes)); + function set.numOfRangeGridPoints(this,numGridPoints) + valid = isscalar(numGridPoints) && numGridPoints > 0; if ~valid matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for combinations! Needs to be one of the strings %s!',strjoin(this.validCombinationTypes,' / ')); + matRad_cfg.dispError('Invalid number of range grid points, needs to be a positive scalar!'); end - this.combinations = combinations_; + this.numOfRAngeGridPoints = inumGridPoints; this.updateScenarios(); end end - - methods (Static) - function [scenarios,linearMask,totNumShiftScen,totNumRangeScen] = createGriddedScenarios(scenarioModel,includeNominalScenario,combinations,combineRange,numOfSetupGridPoints,numOfRangeGridPoints) - matRad_cfg = MatRad_Config.instance(); - - %Get the maximum, i.e., worst case shifts - wcSetupShifts = scenarioModel.wcSigma * scenarioModel.shiftSD; - - %% Create gridded setup shifts - %Create grid vectors for setup shifts - setupShiftGrid = zeros(numOfSetupGridPoints,numel(wcSetupShifts)); - if mod(numOfSetupGridPoints,2) == 0 && includeNominalScenario - matRad_cfg.dispWarning('Obtaining Setup Shifts: Including the nominal scenario with even number of grid points creates asymmetrical shifts!'); - end - - for i = 1:numel(wcSetupShifts) - setupShiftGrid(:,i) = linspace(-wcSetupShifts(i),wcSetupShifts(i),numOfSetupGridPoints); - if includeNominalScenario - - [~,ix] = min(abs(setupShiftGrid(:,i))); - setupShiftGrid(ix,i) = 0; - end - end - - %Now create vector of all shifts for different combinatorial - %settings - switch combinations - case 'none' - %Create independent shifts - griddedSetupShifts = []; - for i=1:size(setupShiftGrid,2) - tmpGrid = zeros(size(setupShiftGrid,1),3); - tmpGrid(:,i) = setupShiftGrid(:,i); - griddedSetupShifts = [griddedSetupShifts; tmpGrid]; - end - case {'shift','all'} - [X,Y,Z] = meshgrid(setupShiftGrid(:,1),setupShiftGrid(:,2),setupShiftGrid(:,3)); - griddedSetupShifts = [X(:), Y(:), Z(:)]; - otherwise - matRad_cfg.dispError('Invalid value for combinations! This sanity check should never be reached!'); - end - - griddedSetupShifts = matRad_ImportanceScenarios.uniqueStableRowsCompat(griddedSetupShifts); - shiftNomScenIx = find(all(griddedSetupShifts == zeros(1,3),2)); - - if ~isempty(shiftNomScenIx) || includeNominalScenario - if ~isempty(shiftNomScenIx) - griddedSetupShifts(shiftNomScenIx,:) = []; - end - griddedSetupShifts = [0 0 0; griddedSetupShifts]; - end - - totNumShiftScen = size(griddedSetupShifts,1); - - %% Create gridded range shifts - %Obtain worst case range shifts - wcRangeShifts = scenarioModel.wcSigma * [scenarioModel.rangeAbsSD scenarioModel.rangeRelSD./100]; - - rangeShiftGrid = zeros(numOfRangeGridPoints,numel(wcRangeShifts)); - if mod(numOfRangeGridPoints,2) == 0 && includeNominalScenario - matRad_cfg.dispWarning('Obtaining Range Shifts: Including the nominal scenario with even number of grid points creates asymmetrical shifts!'); - end - - for i = 1:numel(wcRangeShifts) - rangeShiftGrid(:,i) = linspace(-wcRangeShifts(i),wcRangeShifts(i),numOfRangeGridPoints); - if includeNominalScenario - - [~,ix] = min(abs(rangeShiftGrid(:,i))); - rangeShiftGrid(ix,i) = 0; - end - end - - if combineRange - griddedRangeShifts = rangeShiftGrid; - else - [rngAbs,rngRel] = meshgrid(rangeShiftGrid(:,1),rangeShiftGrid(:,2)); - griddedRangeShifts = [rngAbs(:),rngRel(:)]; - end - - %Remove duplicate scenarios and update number of shifts - griddedRangeShifts = matRad_ImportanceScenarios.uniqueStableRowsCompat(griddedRangeShifts); - - rangeNomScenIx = find(all(griddedRangeShifts == zeros(1,2),2)); - - if ~isempty(rangeNomScenIx) || includeNominalScenario - if ~isempty(rangeNomScenIx) - griddedRangeShifts(rangeNomScenIx,:) = []; - end - griddedRangeShifts = [0 0; griddedRangeShifts]; - end - - totNumRangeScen = size(griddedRangeShifts,1); - - - %Aggregate scenarios - switch combinations - case {'none','shift'} - scenarios = zeros(totNumShiftScen + totNumRangeScen,5); - scenarios(1:totNumShiftScen,1:3) = griddedSetupShifts; - scenarios(totNumShiftScen+1:totNumShiftScen+totNumRangeScen,4:5) = griddedRangeShifts; - - - - %create the linear mask of scenarios - linearMask = ones(size(scenarios,1),3); - linearMask(1:totNumShiftScen,2) = (1:totNumShiftScen)'; - linearMask(totNumShiftScen+1:totNumShiftScen+totNumRangeScen,3) = (1:totNumRangeScen)'; - - [scenarios,ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); - linearMask = linearMask(ia,:); - - case 'all' - %Prepare scenario matrix by replicating shifts - %with the number of range scenarios - scenarios = repmat(griddedSetupShifts,totNumRangeScen,1); - scenarios = [scenarios zeros(size(scenarios,1),2)]; - - %create the linear mask of scenarios - linearMask = ones(size(scenarios,1),3); - for r = 1:totNumRangeScen - offset = (r-1)*totNumShiftScen; - ixR = (offset + 1):(offset + totNumShiftScen); - scenarios(ixR,4:5) = repmat(griddedRangeShifts(r,:),totNumShiftScen,1); - - %Set linear mask - linearMask(ixR,2) = (1:totNumShiftScen)'; - linearMask(ixR,3) = r; - end - - - %create the linear mask of scenarios - [scenarios,ia] = matRad_ImportanceScenarios.uniqueStableRowsCompat(scenarios); - linearMask = linearMask(ia,:); - otherwise - matRad_cfg.dispError('Invalid value for combinations! This sanity check should never be reached!'); - end - end - - function [uniqueStableRows,ia] = uniqueStableRowsCompat(values) - %This is a compatability wrapper to call unique without sorting - - matRad_cfg = MatRad_Config.instance(); - - if matRad_cfg.isOctave - %https://stackoverflow.com/questions/37828894/ - [~,ia,~] = unique(values,'rows','first'); - ia = sort(ia); - uniqueStableRows = values(ia,:); - else - [uniqueStableRows,ia] = unique(values,'rows','stable'); - end - end - end end diff --git a/scenarios/matRad_NominalScenario.m b/scenarios/matRad_NominalScenario.m index 99a08d969..b97666ec5 100644 --- a/scenarios/matRad_NominalScenario.m +++ b/scenarios/matRad_NominalScenario.m @@ -46,11 +46,12 @@ %TODO: In the context of an uncertainty model, we should %consider assigning probability according to the model, and %just leaving the weight 1 - this.scenForProb = [0 0 0 0 0]; - this.scenProb = 1; + this.scenForProb = zeros(this.numOfCtScen,5); + this.scenProb = this.phaseProbability; %Scenario weight - this.scenWeight = 1; + this.scenWeight = ones(this.numOfCtScen,1)./this.numOfCtScen; + this.scenWeight = this.phaseProbability; %set variables this.totNumShiftScen = 1; @@ -58,9 +59,9 @@ this.totNumScen = this.numOfCtScen; %Individual shifts - this.relRangeShift = 0; - this.absRangeShift = 0; - this.isoShift = [0 0 0]; + this.relRangeShift = zeros(this.numOfCtScen,1); + this.absRangeShift = zeros(this.numOfCtScen,1); + this.isoShift = zeros(this.numOfCtScen,3); this.maxAbsRangeShift = max(this.absRangeShift); this.maxRelRangeShift = max(this.absRangeShift); @@ -77,11 +78,18 @@ Sigma = diag([this.shiftSD,this.rangeAbsSD,this.rangeRelSD./100].^2); d = size(Sigma,1); [cs,p] = chol(Sigma); - this.scenProb = (2*pi)^(-d/2) * exp(-0.5*sum((this.scenForProb/cs).^2, 2)) / prod(diag(cs)); + tmpScenProb = (2*pi)^(-d/2) * exp(-0.5*sum((this.scenForProb/cs).^2, 2)) / prod(diag(cs)); + + %Multiply with 4D phase probability + this.scenProb = this.phaseProbability .* tmpScenProb; + + %Get relative (normalized) weight of the scenario this.scenWeight = this.scenProb./sum(this.scenProb); %Return variable - scenarios = [0 0 0 0 0]; + scenarios = (1:this.numOfCtScen)'; + scenarios = [scenarios zeros(this.numOfCtScen,5)]; + this.scenForProb = scenarios; if totNumScen ~= this.totNumScen matRad_cfg = MatRad_Config.instance(); diff --git a/scenarios/matRad_ScenarioModel.m b/scenarios/matRad_ScenarioModel.m index f3e642206..8f5c0c0e7 100644 --- a/scenarios/matRad_ScenarioModel.m +++ b/scenarios/matRad_ScenarioModel.m @@ -1,4 +1,4 @@ -classdef matRad_ScenarioModel < handle +classdef (Abstract) matRad_ScenarioModel < handle % matRad_ScenarioModel % This is an abstract interface class to define Scenario Models for use in % robust treatment planning and uncertainty analysis. @@ -30,7 +30,9 @@ rangeRelSD = 3.5; % given in % rangeAbsSD = 1; % given in [mm] shiftSD = [2.25 2.25 2.25]; % given in [mm] - wcSigma = 1; % Multiplier to compute the worst case / maximum shifts + wcSigma = 1; % Multiplier to compute the worst case / maximum shifts + + phaseProbability = 1; % Probability of being in a phase of a 4D CT (vector with probability value for each phase) end properties (Abstract,SetAccess=protected) @@ -44,7 +46,6 @@ properties (SetAccess = protected) numOfCtScen; % total number of CT scenarios - % these parameters will be filled according to the choosen scenario type isoShift; relRangeShift; @@ -71,6 +72,8 @@ else this.numOfCtScen = ct.numOfCtScen; end + + this.phaseProbability = ones(this.numOfCtScen,1)./this.numOfCtScen; %Equal probability to be in each phase of the 4D ct %TODO: We could do this here automatically in the constructur, but %Octave 5 has a bug here and throws an error @@ -128,7 +131,16 @@ function listAllScenarios(this) this.updateScenarios(); end - + function set.phaseProbability(this,phaseProbability) + valid = isnumeric(phaseProbability) && iscolumn(phaseProbability) && all(phaseProbability <= 1) && all(phaseProbability >= 0); + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for phaseProbability! Needs to be a valid column vector of probabilities [0,1]!'); + end + this.phaseProbability = phaseProbability; + this.updateScenarios(); + end + function scenarios = updateScenarios(this) %This function will always update the scenarios given the diff --git a/scenarios/matRad_WorstCaseScenarios.m b/scenarios/matRad_WorstCaseScenarios.m index 8f2103161..411379c9b 100644 --- a/scenarios/matRad_WorstCaseScenarios.m +++ b/scenarios/matRad_WorstCaseScenarios.m @@ -1,4 +1,4 @@ -classdef matRad_WorstCaseScenarios < matRad_ScenarioModel +classdef matRad_WorstCaseScenarios < matRad_GriddedScenariosAbstract % matRad_WorstCaseScenarios % Implements single worst-case shifts per dimension.% % @@ -22,18 +22,14 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - properties - includeNominalScenario = true; - combinations = 'none'; %Can be 'none', 'shift', 'all' to ontrol creation of worst case combinations - combineRange = true; %Wether to treat absolute & relative range as one shift or as separate scenarios - end properties (SetAccess=protected) name = 'wcScen'; end - properties (Constant) - validCombinationTypes = {'all','none','shift'}; + properties (Hidden) + numOfSetupGridPoints = 2; + numOfRangeGridPoints = 2; end methods @@ -44,80 +40,27 @@ superclassArgs = {ct}; end - this@matRad_ScenarioModel(superclassArgs{:}); + this@matRad_GriddedScenariosAbstract(superclassArgs{:}); %TODO: We could do this automatically in the superclass %Octave 5 has a bug there and throws an error this.updateScenarios(); end - - function set.includeNominalScenario(this,includeNomScen) - valid = isscalar(includeNomScen) && (isnumeric(includeNomScen) || islogical(includeNomScen)); - if ~valid - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for includeNominalScenario! Needs to be a boolean / logical value!'); - end - this.includeNominalScenario = includeNomScen; - this.updateScenarios(); - end - - function set.combineRange(this,combineRange_) - valid = isscalar(combineRange_) && (isnumeric(combineRange_) || islogical(combineRange_)); - if ~valid - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for combineRange! Needs to be a boolean / logical value!'); - end - this.combineRange = combineRange_; - this.updateScenarios(); - end - - function set.combinations(this,combinations_) - valid = any(strcmp(combinations_,this.validCombinationTypes)); - if ~valid - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for combinations! Needs to be one of the strings %s!',strjoin(this.validCombinationTypes,' / ')); - end - this.combinations = combinations_; - this.updateScenarios(); - end function scenarios = updateScenarios(this) if this.includeNominalScenario - nPoints = 3; + this.numOfSetupGridPoints = 3; + this.numOfRangeGridPoints = 3; else - nPoints = 2; + this.numOfSetupGridPoints = 2; + this.numOfRangeGridPoints = 2; end %Use the static gridded shift function from %ImportanceScenarios. We set inclusion of nominal scenarios to %false and handle it automatically via the grid point number - [scenarios,linearMask,this.totNumShiftScen,this.totNumRangeScen] = matRad_ImportanceScenarios.createGriddedScenarios(this,false,this.combinations,this.combineRange,nPoints,nPoints); - - %Finalize meta information - this.totNumScen = size(scenarios,1); - - this.relRangeShift = scenarios(:,5); - this.absRangeShift = scenarios(:,4); - this.isoShift = scenarios(:,1:3); - - this.maxAbsRangeShift = max(this.absRangeShift); - this.maxRelRangeShift = max(this.relRangeShift); - - this.scenMask = false(this.numOfCtScen,this.totNumShiftScen,this.totNumRangeScen); - - this.scenForProb = scenarios; - this.linearMask = linearMask; - - maskIx = sub2ind(size(this.scenMask),linearMask(:,1),linearMask(:,2),linearMask(:,3)); - this.scenMask(maskIx) = true; - - %Get Scenario probability - Sigma = diag([this.shiftSD,this.rangeAbsSD,this.rangeRelSD./100].^2); - d = size(Sigma,1); - [cs,p] = chol(Sigma); - this.scenProb = (2*pi)^(-d/2) * exp(-0.5*sum((scenarios/cs).^2, 2)) / prod(diag(cs)); - this.scenWeight = this.scenProb./sum(this.scenProb); + scenarios = this.updateScenarios@matRad_GriddedScenariosAbstract(); end end