From a3644fc765fb13d724598c7ab73e7f79961c39c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Fri, 12 Jan 2018 17:14:04 +0100 Subject: [PATCH 01/13] Init commit to perform class based simulation analysis. --- tools/samplingAnalysis/MatRadScenario.m | 138 ++++++++++++++++++++++ tools/samplingAnalysis/MatRadSimulation.m | 107 +++++++++++++++++ 2 files changed, 245 insertions(+) create mode 100644 tools/samplingAnalysis/MatRadScenario.m create mode 100644 tools/samplingAnalysis/MatRadSimulation.m diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m new file mode 100644 index 000000000..156655373 --- /dev/null +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -0,0 +1,138 @@ +classdef MatRadScenario < handle + + properties (Access = private) + dvhQiReady = false + % during calcQiDVH + cst + pln + end + + properties (SetAccess = private) + shift; + relRangeShift; + absRangeShift; + ctShiftIdentifier; + weight; + radiationQuantity; + + % computed properties which need to be recalculated when dose + % changes. is not in Dependent properties because calculation is + % expensive + dvh; + qi; + end + + properties (SetObservable, AbortSet, SetAccess = private) + % container of scenarios + dose; + end + + properties (Dependent=true, SetAccess = private) + end + + methods + function obj = MatRadScenario(dose, weight, radiationQuantity, ctShiftIdentifier, shift, ... + relRangeShift, absRangeShift) + + obj.dose = dose; + obj.weight = weight; + obj.radiationQuantity = radiationQuantity; + + obj.ctShiftIdentifier = ctShiftIdentifier; + obj.shift = shift; + obj.relRangeShift = relRangeShift; + obj.absRangeShift = absRangeShift; + + addlistener(obj,'dose','PostSet',@obj.handleChangeOfDose); + end % eof constructor + + function obj = calcQiDVH(obj, cst, pln, dvhType, doseGrid, refGy, refVol) + obj.cst = cst; + obj.pln = pln; + obj.dvh = matRad_calcDVH(cst, obj.dose, dvhType, doseGrid); + obj.qi = matRad_calcQualityIndicators(cst, pln, obj.dose, refGy, refVol); + obj.dvhQiReady = true; + end + + function plotQiDVH(obj) + if obj.dvhQiReady + subplot(2,1,1) + matRad_showDVH(obj.dvh,obj.cst,obj.pln); + subplot(2,1,2) + matRad_showQualityIndicators(obj.qi); + else + error('You need to compute first.'); + end + end + + function plotDoseSlice(obj, ax, ct, plane, slice) + matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice) + end + + function singleStructDVH = getSingleStructDVH(obj, voi) + ix = obj.getIndexOfStruct(obj.dvh, voi, 'name'); + singleStructDVH = obj.dvh(ix); + end + + function singleStructQi = getSingleStructQi(obj, voi) + ix = obj.getIndexOfStruct(obj.qi, voi, 'name'); + singleStructQi = obj.qi(ix); + end + + function dvh = get.dvh(obj) + if isempty(obj.dvh) + error('DVH is not yet calculated.'); + else + dvh = obj.dvh; + end + end + + function qi = get.qi(obj) + if isempty(obj.qi) + error('QI are not yet calculated.'); + else + qi = obj.qi; + end + end + + + end % eof methods + + methods (Access = private) + function obj = handleChangeOfDose(obj,src,data) + obj.dvhQiReady = false; + obj.resetComputedProperties; + fprintf('Changed.'); + end + + function obj = resetComputedProperties(obj) + obj.dvhQiReady = false; + obj.dvh = []; + obj.qi = []; + obj.cst = []; + obj.pln = []; + end + + end + + methods (Static, Access = private) + + function ix = getIndexOfStruct(voiInStruct, voi, fieldname) + if isnumeric(voi) + ix = voi; + elseif ischar(voi) + for i = 1:numel(voiInStruct) + if voiInStruct(i).(fieldname) == voi + ix = i; + break; + end + end + else + ix = NaN; + error('VOI needs to be either string or numeric'); + end + end + + end + +end % eof classdef diff --git a/tools/samplingAnalysis/MatRadSimulation.m b/tools/samplingAnalysis/MatRadSimulation.m new file mode 100644 index 000000000..4b432df46 --- /dev/null +++ b/tools/samplingAnalysis/MatRadSimulation.m @@ -0,0 +1,107 @@ +classdef MatRadSimulation < handle + + properties (Access = private) + % computed properties which need to be recalculated when dose + % changes. is not in Dependent properties because calculation is + % expensive + doseContainer + end + + properties (SetAccess = private) + % should only be changed with constructor + radiationQuantity + subIx + doseCubeDimensions + nominalScenario + cst + + statisticsComputed = false + end + + properties (SetObservable, AbortSet, Access = private) + % container of scenarios + scenContainer + end + + properties (Dependent = true) + numOfScen + dvhContainer + qiContainer + end + + properties (Dependent = true, Access = private) + nextScenIdentifier + end + + methods + function obj = MatRadSimulation(radiationQuantity, nominalScenario, cst, subIx, expectedNumOfScen) + % optionals + if ~exist('expectedNumOfScen', 'var') || isempty(expectedNumOfScen) + expectedNumOfScen = 0; + end + + obj.radiationQuantity = radiationQuantity; + obj.nominalScenario = nominalScenario; + obj.nominalScenario = cst; + obj.subIx = subIx; + + obj.scenContainer = cell(1, expectedNumOfScen); + addlistener(obj,'scenContainer','PostSet',@obj.handleChangeOfScen); + end % eof constructor + + function obj = initNewScen(obj, scenario) + if obj.numOfScen == 0 + obj.deriveConstantsFromFirstScenario(scenario); + end + if ~obj.isValidScen(scenario) + error('Scenario not valid.'); + end + + obj.scenContainer{obj.nextScenIdentifier} = scenario; + end % eof initNewScen + + function obj = runAnalysis(obj) + obj.fillDoseContainer(); + end + + function numOfScen = get.numOfScen(obj) + numOfScen = sum(~cellfun(@isempty,obj.scenContainer)); + end % eof get.NumOfScen + + function numOfScen = get.nextScenIdentifier(obj) + numOfScen = obj.numOfScen + 1; + end % eof get.NumOfScen + end + + methods (Static) + + end + + methods (Access = private) + function validScen = isValidScen(obj, scenario) + if scenario.radiationQuantity ~= obj.radiationQuantity + error('Scenarios can only be added if they are of the same quantity.'); + else + validScen = true; + end + + end % eof isValidScen + + function obj = deriveConstantsFromFirstScenario(obj, scenario) + obj.doseCubeDimensions = size(scenario.dose); + end % eof deriveConstantsFrom + + function obj = handleChangeOfScen(obj,~,~) + obj.statisticsComputed = false; + fprintf('Changed.'); + end + + function obj = fillDoseContainer(obj) + for i=1:obj.numOfScen + obj.doseContainer(:,i) = obj.scenContainer{i}.dose(:); + end + end + + end + +end From b510b3120a95fc2da5f16eedfdda3969c1e12e12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Fri, 12 Jan 2018 17:43:00 +0100 Subject: [PATCH 02/13] Add support for 1D dose --- tools/samplingAnalysis/MatRadScenario.m | 37 +++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 156655373..ac9b037a7 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -2,12 +2,15 @@ properties (Access = private) dvhQiReady = false + doseCubeDim % during calcQiDVH cst pln end properties (SetAccess = private) + + subIx; shift; relRangeShift; absRangeShift; @@ -25,16 +28,26 @@ properties (SetObservable, AbortSet, SetAccess = private) % container of scenarios dose; + doseLin; end properties (Dependent=true, SetAccess = private) end methods - function obj = MatRadScenario(dose, weight, radiationQuantity, ctShiftIdentifier, shift, ... - relRangeShift, absRangeShift) + function obj = MatRadScenario(dose, subIx, weight, radiationQuantity, ctShiftIdentifier, shift, ... + relRangeShift, absRangeShift, doseCubeDim) - obj.dose = dose; + if ndims(dose) == 3 + obj.dose = dose; + obj.subIx = []; + obj.doseCubeDim = size(obj.dose); + else + obj.doseCubeDim = doseCubeDim; + obj.subIx = subIx; + obj.doseLin = dose; + obj.dose = []; + end obj.weight = weight; obj.radiationQuantity = radiationQuantity; @@ -44,6 +57,7 @@ obj.absRangeShift = absRangeShift; addlistener(obj,'dose','PostSet',@obj.handleChangeOfDose); + addlistener(obj,'doseLin','PostSet',@obj.handleChangeOfDose); end % eof constructor function obj = calcQiDVH(obj, cst, pln, dvhType, doseGrid, refGy, refVol) @@ -79,6 +93,18 @@ function plotDoseSlice(obj, ax, ct, plane, slice) singleStructQi = obj.qi(ix); end + function dose = get.dose(obj) + if isempty(obj.dose) + dose = obj.get3dDoseFromLin(obj.doseLin, obj.subIx, obj.doseCubeDim); + else + dose = obj.dose; + end + end + + function set.dose(obj, v) + obj.dose = v; + end + function dvh = get.dvh(obj) if isempty(obj.dvh) error('DVH is not yet calculated.'); @@ -133,6 +159,11 @@ function plotDoseSlice(obj, ax, ct, plane, slice) end end + function dose3D = get3dDoseFromLin(doseLin, subIx, doseCubeDim) + dose3D = zeros(doseCubeDim(1), doseCubeDim(2), doseCubeDim(3)); + dose3D(subIx) = doseLin; + end + end end % eof classdef From 0f3d6c8aa4e2bef8f735dd01a155950556f35876 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Rapahel=20Mu=CC=88ller?= Date: Sun, 14 Jan 2018 21:05:15 +0100 Subject: [PATCH 03/13] Add dvh qi calculation. Add structurewise dvh qi. --- tools/samplingAnalysis/MatRadScenario.m | 39 ++++--- tools/samplingAnalysis/MatRadSimulation.m | 120 ++++++++++++++++++++-- 2 files changed, 137 insertions(+), 22 deletions(-) diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index ac9b037a7..0d826e60e 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -2,21 +2,20 @@ properties (Access = private) dvhQiReady = false - doseCubeDim % during calcQiDVH cst pln end properties (SetAccess = private) - + weight; subIx; shift; relRangeShift; absRangeShift; ctShiftIdentifier; - weight; radiationQuantity; + doseCubeDim; % computed properties which need to be recalculated when dose % changes. is not in Dependent properties because calculation is @@ -28,27 +27,29 @@ properties (SetObservable, AbortSet, SetAccess = private) % container of scenarios dose; - doseLin; end properties (Dependent=true, SetAccess = private) end methods - function obj = MatRadScenario(dose, subIx, weight, radiationQuantity, ctShiftIdentifier, shift, ... - relRangeShift, absRangeShift, doseCubeDim) + function obj = MatRadScenario(dose, subIx, radiationQuantity, ctShiftIdentifier, shift, ... + relRangeShift, absRangeShift, doseCubeDim, weight) + obj.dose = dose; if ndims(dose) == 3 - obj.dose = dose; - obj.subIx = []; + obj.subIx = (1:numel(obj.dose))'; obj.doseCubeDim = size(obj.dose); else obj.doseCubeDim = doseCubeDim; obj.subIx = subIx; - obj.doseLin = dose; - obj.dose = []; end - obj.weight = weight; + + if exist('weight','var') && ~isempty(weight) + obj.weight = weight; + else + obj.weight = NaN; + end obj.radiationQuantity = radiationQuantity; obj.ctShiftIdentifier = ctShiftIdentifier; @@ -57,7 +58,7 @@ obj.absRangeShift = absRangeShift; addlistener(obj,'dose','PostSet',@obj.handleChangeOfDose); - addlistener(obj,'doseLin','PostSet',@obj.handleChangeOfDose); + % addlistener(obj,'doseLin','PostSet',@obj.handleChangeOfDose); end % eof constructor function obj = calcQiDVH(obj, cst, pln, dvhType, doseGrid, refGy, refVol) @@ -95,7 +96,9 @@ function plotDoseSlice(obj, ax, ct, plane, slice) function dose = get.dose(obj) if isempty(obj.dose) - dose = obj.get3dDoseFromLin(obj.doseLin, obj.subIx, obj.doseCubeDim); + dose = []; + elseif ndims(obj.dose) ~= 3 + dose = obj.get3dDoseFromLin(obj.dose, obj.subIx, obj.doseCubeDim); else dose = obj.dose; end @@ -105,6 +108,14 @@ function plotDoseSlice(obj, ax, ct, plane, slice) obj.dose = v; end +% function doseLin = get.doseLin(obj) +% if ~isempty(obj.doseLin) +% doseLin = obj.doseLin; +% else +% doseLin = obj.dose(:); +% end +% end + function dvh = get.dvh(obj) if isempty(obj.dvh) error('DVH is not yet calculated.'); @@ -148,7 +159,7 @@ function plotDoseSlice(obj, ax, ct, plane, slice) ix = voi; elseif ischar(voi) for i = 1:numel(voiInStruct) - if voiInStruct(i).(fieldname) == voi + if strcmp(voiInStruct(i).(fieldname), voi) ix = i; break; end diff --git a/tools/samplingAnalysis/MatRadSimulation.m b/tools/samplingAnalysis/MatRadSimulation.m index 4b432df46..4d3ac175f 100644 --- a/tools/samplingAnalysis/MatRadSimulation.m +++ b/tools/samplingAnalysis/MatRadSimulation.m @@ -10,12 +10,27 @@ properties (SetAccess = private) % should only be changed with constructor radiationQuantity - subIx doseCubeDimensions nominalScenario + ct cst + pln + + % doseStat + meanCube + stdCube + gammaCube + gammaDoseAgreement + gammaDistAgreement + percentiles statisticsComputed = false + + % computed properties which need to be recalculated when dose + % changes. is not in Dependent properties because calculation is + % expensive + dvhContainer + qiContainer end properties (SetObservable, AbortSet, Access = private) @@ -25,16 +40,18 @@ properties (Dependent = true) numOfScen - dvhContainer - qiContainer + weights end properties (Dependent = true, Access = private) nextScenIdentifier + doseCubeDim + weights_unnormalised + subIx end methods - function obj = MatRadSimulation(radiationQuantity, nominalScenario, cst, subIx, expectedNumOfScen) + function obj = MatRadSimulation(radiationQuantity, nominalScenario, ct, cst, pln, expectedNumOfScen) % optionals if ~exist('expectedNumOfScen', 'var') || isempty(expectedNumOfScen) expectedNumOfScen = 0; @@ -42,14 +59,16 @@ obj.radiationQuantity = radiationQuantity; obj.nominalScenario = nominalScenario; - obj.nominalScenario = cst; - obj.subIx = subIx; + obj.ct = ct; + obj.cst = cst; + obj.pln = pln; + obj.scenContainer = NaN * ones(expectedNumOfScen, 1); obj.scenContainer = cell(1, expectedNumOfScen); addlistener(obj,'scenContainer','PostSet',@obj.handleChangeOfScen); end % eof constructor - function obj = initNewScen(obj, scenario) + function obj = initNewScen(obj, scenario, w) if obj.numOfScen == 0 obj.deriveConstantsFromFirstScenario(scenario); end @@ -57,17 +76,46 @@ error('Scenario not valid.'); end + % obj.weights_unnormalised(obj.nextScenIdentifier,1) = w; obj.scenContainer{obj.nextScenIdentifier} = scenario; end % eof initNewScen - function obj = runAnalysis(obj) + function obj = runAnalysis(obj, gammaCriteria, percentiles) + obj.gammaDoseAgreement = gammaCriteria(1); + obj.gammaDistAgreement = gammaCriteria(2); + obj.percentiles = percentiles; + obj.fillDoseContainer(); + obj.computeMeanStdCube(); + obj.computeGammaCube(); + obj.computeAllDvhQi(); + obj.computeStructureWiseDvhQi(); + end + function weights_unnormalised = get.weights_unnormalised(obj) + weights_unnormalised = NaN * ones(obj.numOfScen, 1); + for i = 1:obj.numOfScen + weights_unnormalised(i) = obj.scenContainer{i}.weight; + end + end % eof get.NumOfScen + + function weights = get.weights(obj) + weights = obj.weights_unnormalised ./ sum(obj.weights_unnormalised); + end % eof get.NumOfScen + function numOfScen = get.numOfScen(obj) numOfScen = sum(~cellfun(@isempty,obj.scenContainer)); end % eof get.NumOfScen + + function doseCubeDim = get.doseCubeDim(obj) + doseCubeDim = obj.nominalScenario.doseCubeDim; + end % eof get.doseCubeDim + function subIx = get.subIx(obj) + subIx = obj.nominalScenario.subIx; + end % eof get.subIx + function numOfScen = get.nextScenIdentifier(obj) numOfScen = obj.numOfScen + 1; end % eof get.NumOfScen @@ -78,6 +126,42 @@ end methods (Access = private) + function computeAllDvhQi(obj, refVol, refGy) + if ~exist('refVol', 'var') || isempty(refVol) || ~exist('refGy', 'var') || isempty(refVol) + refVol = [2 5 50 95 98]; + refGy = linspace(0,max(obj.nominalScenario.dose(:)),6); + end + % compute nomQIDVH + obj.nominalScenario.calcQiDVH(obj.cst, obj.pln, 'cum', [], refGy, refVol); + doseGrid = obj.nominalScenario.dvh(1).doseGrid; + for i = 1:obj.numOfScen + obj.scenContainer{i}.calcQiDVH(obj.cst, obj.pln, 'cum', doseGrid, refGy, refVol); + end + end + + function computeStructureWiseDvhQi(obj) + dvhContainer = struct(); + qiContainer = struct(); + % reassing dvh to stats structure + for i = 1:size(obj.cst,1) + dvhContainer(i).name = obj.cst{i,2}; + qiContainer(i).name = obj.cst{i,2}; + dvhContainer(i).weights = obj.weights; + qiContainer(i).weights = obj.weights; + dvh = obj.scenContainer{1}.getSingleStructDVH(dvhContainer(i).name); + dvhContainer(i).doseGrid = NaN * ones(obj.numOfScen, numel(dvh.doseGrid)); + dvhContainer(i).volumePoints = NaN * ones(obj.numOfScen, numel(dvh.volumePoints)); + for j = 1:obj.numOfScen + dvh = obj.scenContainer{j}.getSingleStructDVH(dvhContainer(i).name); + dvhContainer(i).doseGrid(j,:) = dvh.doseGrid; + dvhContainer(i).volumePoints(j,:) = dvh.volumePoints; + qiContainer(i).qi(j) = obj.scenContainer{j}.getSingleStructQi(qiContainer(i).name); + end + end + obj.dvhContainer = dvhContainer; + obj.qiContainer = qiContainer; + end + function validScen = isValidScen(obj, scenario) if scenario.radiationQuantity ~= obj.radiationQuantity error('Scenarios can only be added if they are of the same quantity.'); @@ -93,6 +177,10 @@ function obj = handleChangeOfScen(obj,~,~) obj.statisticsComputed = false; + obj.dvhContainer = struct(); + obj.qiContainer = struct(); + obj.doseContainer = []; + fprintf('Changed.'); end @@ -101,6 +189,22 @@ obj.doseContainer(:,i) = obj.scenContainer{i}.dose(:); end end + + function computeMeanStdCube(obj) + obj.meanCube = zeros(obj.doseCubeDim); + obj.stdCube = zeros(obj.doseCubeDim); + + ix = obj.subIx; + + obj.meanCube(ix) = (sum(obj.doseContainer * diag(obj.weights),2)); + obj.stdCube(ix) = std(obj.doseContainer, obj.weights,2); + end + + function computeGammaCube(obj) + obj.gammaCube = zeros(obj.doseCubeDim); + + obj.gammaCube = matRad_gammaIndex(obj.nominalScenario.dose,obj.meanCube,[obj.ct.resolution.x obj.ct.resolution.y obj.ct.resolution.z],[obj.gammaDoseAgreement obj.gammaDistAgreement]); + end end From 9de2eb9f4ebde342c9c842a832ff11bfad359d64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 15 Jan 2018 17:08:36 +0100 Subject: [PATCH 04/13] Draft to account for fractionation. --- MatRadMultScen.m | 135 +++++++++++++ matRad_sampling.m | 41 ++-- tools/samplingAnalysis/MatRadScenario.m | 28 +-- tools/samplingAnalysis/MatRadSimulation.m | 223 ++++++++++++++++++---- tools/samplingAnalysis/matRad_calcStudy.m | 62 +++--- 5 files changed, 399 insertions(+), 90 deletions(-) create mode 100644 MatRadMultScen.m diff --git a/MatRadMultScen.m b/MatRadMultScen.m new file mode 100644 index 000000000..f88698dca --- /dev/null +++ b/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 diff --git a/matRad_sampling.m b/matRad_sampling.m index f2fb2a656..7e2eda653 100644 --- a/matRad_sampling.m +++ b/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 % @@ -122,6 +122,11 @@ resultGUInomScen.qi = nomQi; resultGUInomScen.cst = cst; +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; @@ -148,6 +153,7 @@ FlagParforProgressDisp = false; end + scenContainer = cell(pln.multScen.totNumScen,1); parfor i = 1:pln.multScen.totNumScen % create nominal scenario @@ -166,15 +172,13 @@ plnSamp.multScen.scenProb = 1; resultSamp = matRad_calcDoseDirect(ct,stf,plnSamp,cst,w,param); - sampledDose = resultSamp.(pln.bioParam.quantityOpt)(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; + sampledDose = resultSamp.(pln.bioParam.quantityVis)(param.subIx); + 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.quantityOpt),'cum',dvhPoints); - caSampRes(i).qi = matRad_calcQualityIndicators(cst,pln,resultSamp.(pln.bioParam.quantityOpt),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; @@ -210,15 +214,14 @@ plnSamp.multScen.scenProb = 1; resultSamp = matRad_calcDoseDirect(ct,stf,plnSamp,cst,w,param); - sampledDose = resultSamp.(pln.bioParam.quantityOpt)(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; + sampledDose = resultSamp.(pln.bioParam.quantityVis)(param.subIx); + 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.quantityOpt),'cum',dvhPoints); - caSampRes(i).qi = matRad_calcQualityIndicators(cst,pln,resultSamp.(pln.bioParam.quantityOpt),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 @@ -227,4 +230,8 @@ %% add subindices pln.subIx = param.subIx; +%% account for fractionation +treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nomScen, ct, cst, pln, pln.numOfFractions); +treatmentSimulation.initFractionatedTreatments(scenContainer, multScen.fraction_index, 'random') + end diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 0d826e60e..0c625c8c2 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -30,6 +30,7 @@ end properties (Dependent=true, SetAccess = private) + doseLin; end methods @@ -38,7 +39,7 @@ obj.dose = dose; if ndims(dose) == 3 - obj.subIx = (1:numel(obj.dose))'; + obj.subIx = subIx; obj.doseCubeDim = size(obj.dose); else obj.doseCubeDim = doseCubeDim; @@ -58,9 +59,9 @@ obj.absRangeShift = absRangeShift; addlistener(obj,'dose','PostSet',@obj.handleChangeOfDose); - % addlistener(obj,'doseLin','PostSet',@obj.handleChangeOfDose); end % eof constructor + function obj = calcQiDVH(obj, cst, pln, dvhType, doseGrid, refGy, refVol) obj.cst = cst; obj.pln = pln; @@ -108,17 +109,20 @@ function plotDoseSlice(obj, ax, ct, plane, slice) obj.dose = v; end -% function doseLin = get.doseLin(obj) -% if ~isempty(obj.doseLin) -% doseLin = obj.doseLin; -% else -% doseLin = obj.dose(:); -% end -% end + function doseLin = get.doseLin(obj) + % is a derived quantity to ensure subIx like linear output of + % dose. + if isempty(obj.dose) + doseLin = []; + else + doseLin = obj.dose(obj.subIx); + end + end function dvh = get.dvh(obj) if isempty(obj.dvh) - error('DVH is not yet calculated.'); + dvh = NaN; + warning('DVH is not yet calculated.'); else dvh = obj.dvh; end @@ -126,7 +130,8 @@ function plotDoseSlice(obj, ax, ct, plane, slice) function qi = get.qi(obj) if isempty(obj.qi) - error('QI are not yet calculated.'); + qi = NaN; + warning('QI are not yet calculated.'); else qi = obj.qi; end @@ -139,7 +144,6 @@ function plotDoseSlice(obj, ax, ct, plane, slice) function obj = handleChangeOfDose(obj,src,data) obj.dvhQiReady = false; obj.resetComputedProperties; - fprintf('Changed.'); end function obj = resetComputedProperties(obj) diff --git a/tools/samplingAnalysis/MatRadSimulation.m b/tools/samplingAnalysis/MatRadSimulation.m index 4d3ac175f..729702458 100644 --- a/tools/samplingAnalysis/MatRadSimulation.m +++ b/tools/samplingAnalysis/MatRadSimulation.m @@ -17,9 +17,10 @@ pln % doseStat - meanCube - stdCube - gammaCube + meanDoseScenario + stdDoseScenario + gammaDoseScenario + gammaDoseAgreement gammaDistAgreement percentiles @@ -29,11 +30,15 @@ % computed properties which need to be recalculated when dose % changes. is not in Dependent properties because calculation is % expensive + dvhDoseGrid dvhContainer qiContainer + + dvhStatistics + qiStatistics end - properties (SetObservable, AbortSet, Access = private) + properties (SetObservable, AbortSet, SetAccess = private) % container of scenarios scenContainer end @@ -41,6 +46,7 @@ properties (Dependent = true) numOfScen weights + metric end properties (Dependent = true, Access = private) @@ -64,11 +70,11 @@ obj.pln = pln; obj.scenContainer = NaN * ones(expectedNumOfScen, 1); - obj.scenContainer = cell(1, expectedNumOfScen); + obj.scenContainer = cell(expectedNumOfScen, 1); addlistener(obj,'scenContainer','PostSet',@obj.handleChangeOfScen); end % eof constructor - function obj = initNewScen(obj, scenario, w) + function obj = initNewScen(obj, scenario) if obj.numOfScen == 0 obj.deriveConstantsFromFirstScenario(scenario); end @@ -80,6 +86,44 @@ obj.scenContainer{obj.nextScenIdentifier} = scenario; end % eof initNewScen + function obj = initFractionatedTreatments(obj, scenarios, treatmentNumber, samplingMethod) + if ~strcmp(samplingMethod, 'random') + error('Only random sampling possible at the moment.') + end + if obj.numOfScen == 0 + obj.deriveConstantsFromFirstScenario(scenarios{1}); + end + if size(treatmentNumber,1) ~= 1 % must be row vector to be iterable + treatmentNumber = treatmentNumber'; + end + treatmentNumber_unique = unique(treatmentNumber); + + % go through complete treatment scenarios + for treatment_currIx = treatmentNumber_unique + % select all scenarios which belong to current treatment_ix + scenForCurrentTreatment = find(treatmentNumber == treatment_currIx); + cumDose = zeros('like',scenarios{1}.doseLin); + ctIndex = NaN * ones(numel(scenForCurrentTreatment),1); + shifts = NaN * ones(numel(scenForCurrentTreatment),3); + relRangeShifts = NaN * ones(numel(scenForCurrentTreatment),1); + absRangeShifts = NaN * ones(numel(scenForCurrentTreatment),1); + runIx = 1; + for s = scenForCurrentTreatment + if ~obj.isValidScen(scenarios{s}) + error('Scenario not valid.'); + end + cumDose = cumDose + scenarios{s}.doseLin; + ctIndex(runIx) = scenarios{s}.ctShiftIdentifier; + shifts(runIx, :) = scenarios{s}.shift; + relRangeShifts(runIx) = scenarios{s}.relRangeShift; + absRangeShifts(runIx) = scenarios{s}.absRangeShift; + runIx = runIx + 1; + end + obj.scenContainer{obj.nextScenIdentifier} = MatRadScenario(cumDose, scenarios{1}.subIx, 'phsicalDose', ctIndex, shifts, ... + relRangeShifts, absRangeShifts, scenarios{1}.doseCubeDim, 1); + end + end + function obj = runAnalysis(obj, gammaCriteria, percentiles) obj.gammaDoseAgreement = gammaCriteria(1); obj.gammaDistAgreement = gammaCriteria(2); @@ -90,7 +134,7 @@ obj.computeGammaCube(); obj.computeAllDvhQi(); obj.computeStructureWiseDvhQi(); - + obj.computeStructureWiseStatistics(); end function weights_unnormalised = get.weights_unnormalised(obj) @@ -119,10 +163,36 @@ function numOfScen = get.nextScenIdentifier(obj) numOfScen = obj.numOfScen + 1; end % eof get.NumOfScen + + function metric = get.metric(obj) + percentileNames = cell(numel(obj.percentiles),1); + % create fieldnames + for i = 1:numel(obj.percentiles) + percentileNames{i} = ['P',num2str(obj.percentiles(i)*100)]; + end + % create table rownames + metric = vertcat({'mean';'min';'max';'std'},percentileNames{:}); + end end methods (Static) + function S = wMean(X,w) + if exist('w','var') || ~isempty(w) + if isvector(X) && isvector(w) + S = reshape(w,1,[]) * reshape(X,[],1) / sum(w); + else + % row-wise + S = reshape(w,1,[]) * X ./ sum(w); + end + + else + S = mean(X); + end + end % eof wMean + + + end methods (Access = private) @@ -137,34 +207,65 @@ function computeAllDvhQi(obj, refVol, refGy) for i = 1:obj.numOfScen obj.scenContainer{i}.calcQiDVH(obj.cst, obj.pln, 'cum', doseGrid, refGy, refVol); end + obj.dvhDoseGrid = doseGrid; end function computeStructureWiseDvhQi(obj) - dvhContainer = struct(); - qiContainer = struct(); + obj.dvhContainer = struct(); + obj.qiContainer = struct(); % reassing dvh to stats structure for i = 1:size(obj.cst,1) - dvhContainer(i).name = obj.cst{i,2}; - qiContainer(i).name = obj.cst{i,2}; - dvhContainer(i).weights = obj.weights; - qiContainer(i).weights = obj.weights; - dvh = obj.scenContainer{1}.getSingleStructDVH(dvhContainer(i).name); - dvhContainer(i).doseGrid = NaN * ones(obj.numOfScen, numel(dvh.doseGrid)); - dvhContainer(i).volumePoints = NaN * ones(obj.numOfScen, numel(dvh.volumePoints)); + obj.dvhContainer(i).name = obj.cst{i,2}; + obj.qiContainer(i).name = obj.cst{i,2}; + %obj.dvhContainer(i).weights = obj.weights; + %obj.qiContainer(i).weights = obj.weights; + dvh = obj.scenContainer{1}.getSingleStructDVH(obj.dvhContainer(i).name); + obj.dvhContainer(i).doseGrid = NaN * ones(1, numel(dvh.doseGrid)); + obj.dvhContainer(i).volumePoints = NaN * ones(obj.numOfScen, numel(dvh.volumePoints)); for j = 1:obj.numOfScen - dvh = obj.scenContainer{j}.getSingleStructDVH(dvhContainer(i).name); - dvhContainer(i).doseGrid(j,:) = dvh.doseGrid; - dvhContainer(i).volumePoints(j,:) = dvh.volumePoints; - qiContainer(i).qi(j) = obj.scenContainer{j}.getSingleStructQi(qiContainer(i).name); + dvh = obj.scenContainer{j}.getSingleStructDVH(obj.dvhContainer(i).name); + if j == 1 + obj.dvhContainer(i).doseGrid(1,:) = dvh.doseGrid; + else + if obj.dvhContainer(i).doseGrid(1,:) ~= dvh.doseGrid + error('Dose grids are not equal.') + end + end + obj.dvhContainer(i).volumePoints(j,:) = dvh.volumePoints; + obj.qiContainer(i).qi(j) = obj.scenContainer{j}.getSingleStructQi(obj.qiContainer(i).name); end end - obj.dvhContainer = dvhContainer; - obj.qiContainer = qiContainer; + end + + function computeStructureWiseStatistics(obj) + % create statstics where structure based results (QI and DVH) are available + for i = 1:numel(obj.dvhContainer) + obj.dvhStatistics(i).VOIname = obj.dvhContainer(i).name; + obj.dvhStatistics(i).stat = obj.calcDVHStat(obj.dvhContainer(i).volumePoints, obj.percentiles, obj.weights); + + obj.qiStatistics(i).name = obj.qiContainer(i).name; + obj.qiStatistics(i).stat = obj.calcQiStat(obj.qiContainer(i).qi, obj.percentiles, obj.weights); + end end function validScen = isValidScen(obj, scenario) - if scenario.radiationQuantity ~= obj.radiationQuantity + if ~strcmp(scenario.radiationQuantity,obj.radiationQuantity) error('Scenarios can only be added if they are of the same quantity.'); + validScen = false; + else + validScen = true; + end + + if ~(scenario.subIx == obj.subIx) + error('Scenarios can only be added if they feature the same sub indices'); + validScen = false; + else + validScen = true; + end + + if ~(scenario.doseCubeDim == obj.doseCubeDimensions) + error('Scenarios can only be added if their cube dimension agree.'); + validScen = false; else validScen = true; end @@ -172,13 +273,15 @@ function computeStructureWiseDvhQi(obj) end % eof isValidScen function obj = deriveConstantsFromFirstScenario(obj, scenario) - obj.doseCubeDimensions = size(scenario.dose); + obj.doseCubeDimensions = size(scenario.dose); end % eof deriveConstantsFrom function obj = handleChangeOfScen(obj,~,~) obj.statisticsComputed = false; obj.dvhContainer = struct(); + obj.dvhStatistics = struct(); obj.qiContainer = struct(); + obj.qiStatistics = struct(); obj.doseContainer = []; fprintf('Changed.'); @@ -186,26 +289,80 @@ function computeStructureWiseDvhQi(obj) function obj = fillDoseContainer(obj) for i=1:obj.numOfScen - obj.doseContainer(:,i) = obj.scenContainer{i}.dose(:); + obj.doseContainer(:,i) = obj.scenContainer{i}.doseLin; end end function computeMeanStdCube(obj) - obj.meanCube = zeros(obj.doseCubeDim); - obj.stdCube = zeros(obj.doseCubeDim); + meanCube = zeros(obj.doseCubeDim); + stdCube = zeros(obj.doseCubeDim); ix = obj.subIx; - obj.meanCube(ix) = (sum(obj.doseContainer * diag(obj.weights),2)); - obj.stdCube(ix) = std(obj.doseContainer, obj.weights,2); + meanCube(ix) = (sum(obj.doseContainer * diag(obj.weights),2)); + stdCube(ix) = std(obj.doseContainer, obj.weights,2); + + obj.meanDoseScenario = MatRadScenario(meanCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN) + obj.stdDoseScenario = MatRadScenario(stdCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN) end function computeGammaCube(obj) - obj.gammaCube = zeros(obj.doseCubeDim); - - obj.gammaCube = matRad_gammaIndex(obj.nominalScenario.dose,obj.meanCube,[obj.ct.resolution.x obj.ct.resolution.y obj.ct.resolution.z],[obj.gammaDoseAgreement obj.gammaDistAgreement]); + gammaCube = zeros(obj.doseCubeDim); + gammaCube = matRad_gammaIndex(obj.nominalScenario.dose,obj.meanDoseScenario.dose,[obj.ct.resolution.x obj.ct.resolution.y obj.ct.resolution.z],[obj.gammaDoseAgreement obj.gammaDistAgreement]); + + obj.gammaDoseScenario = MatRadScenario(gammaCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); end - + + function dvhStat = calcDVHStat(obj, volumePoints, percentiles,w) + dvhMat = volumePoints; % rows are different scenarios, columns positions + % for statistical reasons, treat NaN as 0 + dvhMat(isnan(dvhMat)) = 0; + + dvhStat.mean.volumePoints = obj.wMean(dvhMat,w); + dvhStat.min.volumePoints = min(dvhMat); + dvhStat.max.volumePoints = max(dvhMat); + dvhStat.std.volumePoints = std(dvhMat,w); + + dvhStat.percDVH = NaN * ones(numel(percentiles),size(volumePoints, 2)); + + for j = 1:size(dvhMat,2) + wQ = matRad_weightedQuantile(dvhMat(:,j), percentiles, w', false, 'none'); + dvhStat.percDVH(:,j) = wQ; + end + + end % eof calcDVHStat + + function qiStat = calcQiStat(obj, qi,percentiles,w) + fields = fieldnames(qi); + % remove name field + if sum(strcmp('name', fields)) >= 1 + qi = rmfield(qi, 'name'); + end + fields = fieldnames(qi); + qiStruct = qi; + + % create helper matlab structure which will be converted to table + qiStatH = struct(); + for j = 1:numel(fields) + if numel([qiStruct(:).(fields{j})]) == numel(w) + qiStatH(1).(fields{j}) = obj.wMean([qiStruct(:).(fields{j})],w); + qiStatH(2).(fields{j}) = min([qiStruct(:).(fields{j})]); + qiStatH(3).(fields{j}) = max([qiStruct(:).(fields{j})]); + qiStatH(4).(fields{j}) = std([qiStruct(:).(fields{j})],w); + wQ = matRad_weightedQuantile([qiStruct(:).(fields{j})], percentiles, w', false, 'none'); + for k = 1:numel(wQ) + sIx = k + 4; + qiStatH(sIx).(fields{j}) = wQ(k); + end + else + for k = 1:(4 + numel(percentiles)) + qiStatH(k).(fields{j}) = []; + end + end + end + qiStat = struct2table(qiStatH); + qiStat.Properties.RowNames = obj.metric; + end % eof calcQiStat end end diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index 3fa0a913c..3f1cbb925 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -44,9 +44,15 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) if ~isfield(param,'outputPath') param.outputPath = mfilename('fullpath'); end +if ~isfield(param, 'percentiles') + param.percentiles = [0.01 0.05 0.125 0.25 0.5 0.75 0.875 0.95 0.99]; +end +if ~isfield(param, 'criteria') + param.gammaCriteria = [2 2]; +end % require minimum number of scenarios to ensure proper statistics -if multScen.numOfRangeShiftScen + sum(multScen.numOfShiftScen) < 20 +if multScen.totNumScen < 20 % multScen.numOfRangeShiftScen + sum(multScen.numOfShiftScen) < 20 matRad_dispToConsole('Detected a low number of scenarios. Proceeding is not recommended.',param,'warning'); param.sufficientStatistics = false; pause(1); @@ -95,7 +101,7 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) %% perform calculation and save tic -[caSampRes, mSampDose, pln, resultGUInomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,multScen,param); +[treatmentSimulation, scenContainer, pln, resultGUInomScen, nomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,multScen,param); param.computationTime = toc; param.reportPath = fullfile('report','data'); @@ -104,30 +110,30 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) %% perform analysis % start here loading resultSampling.mat if something went wrong during analysis or report generation -[structureStat, doseStat, param] = matRad_samplingAnalysis(ct,cst,pln,caSampRes,mSampDose,resultGUInomScen,param); - -%% generate report -listOfQI = {'mean', 'std', 'max', 'min', 'D_2', 'D_5', 'D_50', 'D_95', 'D_98'}; - -cd(param.outputPath) -mkdir(fullfile('report','data')); -mkdir(fullfile('report','data','frames')); -mkdir(fullfile('report','data','figures')); -copyfile(fullfile(matRadPath,'tools','samplingAnalysis','main_template.tex'),fullfile('report','main.tex')); - -% generate actual latex report -matRad_latexReport(ct, cst, pln, resultGUInomScen, structureStat, doseStat, mSampDose, listOfQI, param); - -cd('report'); -if ispc - executeLatex = 'xelatex --shell-escape --interaction=nonstopmode main.tex'; -elseif isunix - executeLatex = '/Library/TeX/texbin/xelatex --shell-escape --interaction=nonstopmode main.tex'; -end +treatmentSimulation.runAnalysis(param.gammaCriteria, param.percentiles); -response = system(executeLatex); -if response == 127 % means not found - warning('Could not find tex distribution. Please compile manually.'); -else - system(executeLatex); -end +% %% generate report +% listOfQI = {'mean', 'std', 'max', 'min', 'D_2', 'D_5', 'D_50', 'D_95', 'D_98'}; +% +% cd(param.outputPath) +% mkdir(fullfile('report','data')); +% mkdir(fullfile('report','data','frames')); +% mkdir(fullfile('report','data','figures')); +% copyfile(fullfile(matRadPath,'tools','samplingAnalysis','main_template.tex'),fullfile('report','main.tex')); +% +% % generate actual latex report +% matRad_latexReport(ct, cst, pln, resultGUInomScen, structureStat, doseStat, mSampDose, listOfQI, param); +% +% cd('report'); +% if ispc +% executeLatex = 'xelatex --shell-escape --interaction=nonstopmode main.tex'; +% elseif isunix +% executeLatex = '/Library/TeX/texbin/xelatex --shell-escape --interaction=nonstopmode main.tex'; +% end +% +% response = system(executeLatex); +% if response == 127 % means not found +% warning('Could not find tex distribution. Please compile manually.'); +% else +% system(executeLatex); +% end From 569148f4cf0e4c269ac8214066e785eac1d6ed0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 12 Feb 2018 13:20:04 +0100 Subject: [PATCH 05/13] Include plotDVH band. --- matRad_multScen.m | 11 ------- matRad_sampling.m | 25 ++++++++------- plotting/matRad_plotDVHBand.m | 26 ++++++++------- tools/samplingAnalysis/MatRadScenario.m | 12 +++++-- tools/samplingAnalysis/MatRadSimulation.m | 39 ++++++++++++++++------- 5 files changed, 65 insertions(+), 48 deletions(-) diff --git a/matRad_multScen.m b/matRad_multScen.m index 4e370c7bb..6b37b0fb0 100644 --- a/matRad_multScen.m +++ b/matRad_multScen.m @@ -538,14 +538,3 @@ function attachListener(obj) end % end of matRad_multScen class - - - - - - - - - - - diff --git a/matRad_sampling.m b/matRad_sampling.m index 7e2eda653..09fc4fc13 100644 --- a/matRad_sampling.m +++ b/matRad_sampling.m @@ -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 @@ -99,13 +99,13 @@ matRad_dispToConsole(['matRad: Realizations variable will need: ' num2str(StorageInfo.bytes/1e9) ' GB \n'],param,'info'); % check if parallel toolbox is installed and license can be checked out -try - ver('distcomp') - FlagParallToolBoxLicensed = license('test','Distrib_Computing_Toolbox'); -catch - FlagParallToolBoxLicensed = false; -end - +% try +% ver('distcomp') +% FlagParallToolBoxLicensed = license('test','Distrib_Computing_Toolbox'); +% catch +% FlagParallToolBoxLicensed = false; +% end +FlagParallToolBoxLicensed = false; %% calculate nominal scenario nomScenTimer = tic; resultGUInomScen = matRad_calcDoseDirect(ct,stf,plnNominal,cst,w,param); @@ -131,6 +131,7 @@ 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. @@ -152,8 +153,7 @@ fprintf('matRad: Consider downloading parfor_progress function from the matlab central fileexchange to get feedback from parfor loop.\n'); FlagParforProgressDisp = false; end - - scenContainer = cell(pln.multScen.totNumScen,1); + parfor i = 1:pln.multScen.totNumScen % create nominal scenario @@ -232,6 +232,7 @@ %% account for fractionation treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nomScen, ct, cst, pln, pln.numOfFractions); -treatmentSimulation.initFractionatedTreatments(scenContainer, multScen.fraction_index, 'random') - +for i = 1:numel(scenContainer) + treatmentSimulation.initNewScen(scenContainer{i}) +end end diff --git a/plotting/matRad_plotDVHBand.m b/plotting/matRad_plotDVHBand.m index a6b92bbba..fb397a7fe 100644 --- a/plotting/matRad_plotDVHBand.m +++ b/plotting/matRad_plotDVHBand.m @@ -1,4 +1,4 @@ -function matRad_plotDVHBand(nominalDVH, structureStat, doseLabel) +function matRad_plotDVHBand(doseGrid, nominalDVH, structureStat, percentiles, doseLabel, figureName) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % matRad_plotDVHBand to plot dose volume bands % @@ -30,20 +30,22 @@ function matRad_plotDVHBand(nominalDVH, structureStat, doseLabel) if ~exist('doseLabel', 'var') || isempty(doseLabel) doseLabel = 'a.u.'; end +if ~exist('figureName', 'var') || isempty(figureName) + figure; +else + figure('Name', figureName); +end -figure; -numOfConf = floor(size(structureStat.dvhStat.percDVH,1) / 2); -% DVH -doseGrid = structureStat.dvhStat.mean.doseGrid; +numOfConf = floor(numel(percentiles) / 2); % plot nominal plan [y, argmin] = cutAtArgmin(nominalDVH.volumePoints); -x = nominalDVH.doseGrid(1:argmin); +x = nominalDVH.doseGrid(1:argmin); h(1) = plot(x,y,'LineWidth',2, 'Color', 'k', 'DisplayName', 'nominal'); hold on; % plot mean -[y, argmin] = cutAtArgmin(structureStat.dvhStat.mean.volumePoints); -x = structureStat.dvhStat.mean.doseGrid(1:argmin); +[y, argmin] = cutAtArgmin(structureStat.mean.volumePoints); +x = doseGrid(1:argmin); h(2) = plot(x,y,'--','LineWidth',2, 'Color', 'k', 'DisplayName', '\mu'); % plot dvh confidence bands % colors @@ -54,10 +56,10 @@ function matRad_plotDVHBand(nominalDVH, structureStat, doseLabel) for j = 1:numOfConf hIx = hIx + 1; lIx = j; - hIx = size(structureStat.dvhStat.percDVH,1) - (j-1); - lowerLimit = structureStat.dvhStat.percDVH(lIx,:); - upperLimit = structureStat.dvhStat.percDVH(hIx,:); - confIn = structureStat.percentiles(hIx) - structureStat.percentiles(lIx); + hIx = size(structureStat.percDVH,1) - (j-1); + lowerLimit = structureStat.percDVH(lIx,:); + upperLimit = structureStat.percDVH(hIx,:); + confIn = percentiles(hIx) - percentiles(lIx); confName = ['C', num2str(round(confIn * 100,0))]; h(hIx) = matRad_shadowPlot(doseGrid, lowerLimit, upperLimit, colors(j,:), confName, alphaTrans); end diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 0c625c8c2..08bec81c9 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -81,8 +81,8 @@ function plotQiDVH(obj) end end - function plotDoseSlice(obj, ax, ct, plane, slice) - matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice) + function plotDoseSlice(obj, ax, ct, plane, slice,legendOn) + matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice,[],[],jet,[],[],[],[],[],legendOn); end function singleStructDVH = getSingleStructDVH(obj, voi) @@ -181,4 +181,12 @@ function plotDoseSlice(obj, ax, ct, plane, slice) end + methods (Static) + function cubes = computeDerivedCubes(physicalDose, alpha, beta, bioParam) + RBExD = physicalDose * dij.RBE; + %%% WRITE HERE %%% + + end + end + end % eof classdef diff --git a/tools/samplingAnalysis/MatRadSimulation.m b/tools/samplingAnalysis/MatRadSimulation.m index 729702458..3e9bca1dd 100644 --- a/tools/samplingAnalysis/MatRadSimulation.m +++ b/tools/samplingAnalysis/MatRadSimulation.m @@ -119,7 +119,7 @@ absRangeShifts(runIx) = scenarios{s}.absRangeShift; runIx = runIx + 1; end - obj.scenContainer{obj.nextScenIdentifier} = MatRadScenario(cumDose, scenarios{1}.subIx, 'phsicalDose', ctIndex, shifts, ... + obj.scenContainer{obj.nextScenIdentifier} = MatRadScenario(cumDose, scenarios{1}.subIx, obj.radiationQuantity, ctIndex, shifts, ... relRangeShifts, absRangeShifts, scenarios{1}.doseCubeDim, 1); end end @@ -135,6 +135,7 @@ obj.computeAllDvhQi(); obj.computeStructureWiseDvhQi(); obj.computeStructureWiseStatistics(); + obj.statisticsComputed = true; end function weights_unnormalised = get.weights_unnormalised(obj) @@ -173,6 +174,26 @@ % create table rownames metric = vertcat({'mean';'min';'max';'std'},percentileNames{:}); end + + % plot + function plotDVHband(obj, listOfStruct) + if ~obj.statisticsComputed + warning('Please run statistical analysis first.'); + else + + for i = 1:numel(listOfStruct) + for j = 1:numel(obj.dvhStatistics) + if strcmp(obj.nominalScenario.dvh(j).name, listOfStruct{i}) + nomDVH = obj.nominalScenario.dvh(j); + end + if strcmp(obj.dvhStatistics(j).name, listOfStruct{i}) + structStat = obj.dvhStatistics(j).stat; + end + end + matRad_plotDVHBand(obj.dvhDoseGrid, nomDVH, structStat, obj.percentiles, obj.radiationQuantity, listOfStruct{i}) + end + end + end end methods (Static) @@ -189,10 +210,6 @@ S = mean(X); end end % eof wMean - - - - end methods (Access = private) @@ -240,7 +257,7 @@ function computeStructureWiseDvhQi(obj) function computeStructureWiseStatistics(obj) % create statstics where structure based results (QI and DVH) are available for i = 1:numel(obj.dvhContainer) - obj.dvhStatistics(i).VOIname = obj.dvhContainer(i).name; + obj.dvhStatistics(i).name = obj.dvhContainer(i).name; obj.dvhStatistics(i).stat = obj.calcDVHStat(obj.dvhContainer(i).volumePoints, obj.percentiles, obj.weights); obj.qiStatistics(i).name = obj.qiContainer(i).name; @@ -250,22 +267,22 @@ function computeStructureWiseStatistics(obj) function validScen = isValidScen(obj, scenario) if ~strcmp(scenario.radiationQuantity,obj.radiationQuantity) - error('Scenarios can only be added if they are of the same quantity.'); validScen = false; + error('Scenarios can only be added if they are of the same quantity.'); else validScen = true; end if ~(scenario.subIx == obj.subIx) - error('Scenarios can only be added if they feature the same sub indices'); validScen = false; + error('Scenarios can only be added if they feature the same sub indices'); else validScen = true; end if ~(scenario.doseCubeDim == obj.doseCubeDimensions) - error('Scenarios can only be added if their cube dimension agree.'); validScen = false; + error('Scenarios can only be added if their cube dimension agree.'); else validScen = true; end @@ -302,8 +319,8 @@ function computeMeanStdCube(obj) meanCube(ix) = (sum(obj.doseContainer * diag(obj.weights),2)); stdCube(ix) = std(obj.doseContainer, obj.weights,2); - obj.meanDoseScenario = MatRadScenario(meanCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN) - obj.stdDoseScenario = MatRadScenario(stdCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN) + obj.meanDoseScenario = MatRadScenario(meanCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); + obj.stdDoseScenario = MatRadScenario(stdCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); end function computeGammaCube(obj) From b19fd5947f7da6e87072d02fd916f817c3fb1a11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Thu, 1 Mar 2018 16:41:30 +0100 Subject: [PATCH 06/13] Prototype of correlation break. NomScen breaks currently. --- matRad.m | 6 +-- matRad_multScen.m | 48 +++++++++++++++++++- tools/samplingAnalysis/MatRadScenario.m | 3 ++ tools/samplingAnalysis/matRad_calcStudy.m | 39 +++++++++++++++- tools/samplingAnalysis/setupStudy_template.m | 24 ++++++---- 5 files changed, 105 insertions(+), 15 deletions(-) diff --git a/matRad.m b/matRad.m index 3f37b7d1c..fbb4f8848 100644 --- a/matRad.m +++ b/matRad.m @@ -22,10 +22,10 @@ % load patient data, i.e. ct, voi, cst %load HEAD_AND_NECK -load TG119.mat +% load TG119.mat %load PROSTATE.mat %load LIVER.mat -%load BOXPHANTOM.mat +load BOXPHANTOM_TINY.mat % meta information for treatment plan pln.bixelWidth = 5; % [mm] / also corresponds to lateral spot spacing for particles @@ -35,7 +35,7 @@ pln.numOfVoxels = prod(ct.cubeDim); pln.isoCenter = ones(pln.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0); pln.voxelDimensions = ct.cubeDim; -pln.radiationMode = 'photons'; % either photons / protons / carbon +pln.radiationMode = 'protons'; % either photons / protons / carbon pln.bioOptimization = 'none_physicalDose'; % none_physicalDose: physical optimization; constRBE_RBExD; constant RBE of 1.1; % MCN_effect; McNamara-variable RBE model for protons (effect based) MCN_RBExD; McNamara-variable RBE model for protons (RBExD) based % WED_effect; Wedenberg-variable RBE model for protons (effect based) WED_RBExD; Wedenberg-variable RBE model for protons (RBExD) based diff --git a/matRad_multScen.m b/matRad_multScen.m index 6b37b0fb0..8f28da4b1 100644 --- a/matRad_multScen.m +++ b/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)); diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 08bec81c9..43ce49663 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -61,6 +61,9 @@ addlistener(obj,'dose','PostSet',@obj.handleChangeOfDose); end % eof constructor + function obj = cumulateDose(obj, dose) + obj.dose = obj.dose + dose; + end function obj = calcQiDVH(obj, cst, pln, dvhType, doseGrid, refGy, refVol) obj.cst = cst; diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index 3f1cbb925..ef45a9866 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -32,6 +32,19 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [treatmentSimulation, scenContainer] = mergeSplittedScen(scenContainerParts) + scenContainer = scenContainerParts{1}; + for j = 1:numel(scenContainer) + for k = 2:numel(scenContainerParts) + scenContainer{j}.cumulateDose(scenContainerParts{k}{j}); + end + end + treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nomScen, ct, cst, pln, pln.numOfFractions); + for j = 1:numel(scenContainer) + treatmentSimulation.initNewScen(scenContainer{j}) + end +end + if exist('param','var') if ~isfield(param,'logLevel') param.logLevel = 4; @@ -100,8 +113,32 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) pln.sampling = true; %% perform calculation and save +multScen = multScen.matRad_createValidInstance(); + tic -[treatmentSimulation, scenContainer, pln, resultGUInomScen, nomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,multScen,param); +if strcmp(multScen.deepestCorrelationBreak, 'fraction') + [treatmentSimulation, scenContainer, pln, resultGUInomScen, nomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,multScen,param); +elseif strcmp(multScen.deepestCorrelationBreak, 'beam') + scenContainerParts = cell(numel(stf),1); + for i = 1:numel(stf) + stfPart = stf(i); + plnPart = pln; + plnPart.isoCenter = pln.isoCenter(i,:); + if strcmp(multScen.rangeShiftCorrelationBreak, 'beam') && strcmp(multScen.isoShiftCorrelationBreak, 'beam') + multScen = multScen.matRad_recreateInstance; + elseif strcmp(multScen.isoShiftCorrelationBreak, 'beam') + multScen = multScen.matRad_recreateInstance('range'); + elseif strcmp(multScen.rangeShiftCorrelationBreak, 'beam') + multScen = multScen.matRad_recreateInstance('setup'); + else + erorr('Incosistency'); + end + + [~, scenContainerPart, pln, ~, ~] = matRad_sampling(ct,stfPart,cst,plnPart,resultGUI.w,structSel,multScen,param); + scenContainerParts{i} = scenContainerPart; + end + +end param.computationTime = toc; param.reportPath = fullfile('report','data'); diff --git a/tools/samplingAnalysis/setupStudy_template.m b/tools/samplingAnalysis/setupStudy_template.m index 4f9697206..50ed81305 100644 --- a/tools/samplingAnalysis/setupStudy_template.m +++ b/tools/samplingAnalysis/setupStudy_template.m @@ -4,18 +4,22 @@ multScen = matRad_multScen([],'impScen'); % a) define shift scenarios -multScen.numOfShiftScen = [0 0 0]; % number of shifts in x y and z direction -multScen.shiftSize = [3 3 3]; % maximum shift [mm] % (e.g. prostate cases 5mm otherwise 3mm) -multScen.shiftGenType = 'equidistant'; % equidistant: equidistant shifts, sampled: sample shifts from normal distribution -multScen.shiftCombType = 'individual'; % individual: no combination of shift scenarios; number of shift scenarios is sum(multScen.numOfShiftScen) +multScen.numOfShiftScen = [2 2 2]; % number of shifts in x y and z direction +multScen.shiftSize = [4.5 4.5 4.5]; % maximum shift [mm] % (e.g. prostate cases 5mm otherwise 3mm) +multScen.shiftGenType = 'sampled'; % equidistant: equidistant shifts, sampled: sample shifts from normal distribution +multScen.shiftCombType = 'permuted'; % individual: no combination of shift scenarios; number of shift scenarios is sum(multScen.numOfShiftScen) % permuted: create every possible shift combination; number of shift scenarios is 8,27,64 ... +multScen.isoShiftCorrelationBreak = 'fraction'; + % b) define range error scenarios multScen.numOfRangeShiftScen = 30; % number of absolute and/or relative range scnearios. % if absolute and relative range scenarios are defined then multScen.rangeCombType defines the resulting number of range scenarios multScen.maxAbsRangeShift = 2; % maximum absolute over and undershoot in mm multScen.maxRelRangeShift = 4.5; % maximum relative over and undershoot in % multScen.rangeCombType = 'combined'; % individual: no combination of absolute and relative range scenarios; combined: combine absolute and relative range scenarios -multScen.rangeGenType = 'equidistant'; % equidistant: equidistant range shifts, sampled: sample range shifts from normal distribution +multScen.rangeGenType = 'sampled'; % equidistant: equidistant range shifts, sampled: sample range shifts from normal distribution +mutlScen.rangeShiftCorrelationBreak = 'beam'; + multScen.scenCombType = 'individual'; % 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 @@ -23,11 +27,11 @@ % define standard deviation of normal distribution - important for probabilistic treatment planning -multScen.rangeRelSD = 1.8; % given in [%] -multScen.rangeAbsSD = 0.8; % given in [mm] -multScen.shiftSD = [2 2 2]; % given in [mm] +multScen.rangeRelSD = 0; % given in [%] +multScen.rangeAbsSD = 1.8; % given in [mm] +multScen.shiftSD = [1.8 1.8 1.8]; % given in [mm] + -multScen = multScen.matRad_createValidInstance(); %% path for output pdf and mat param.outputPath = pwd; @@ -42,6 +46,6 @@ % param.criteria = [3 3]; %%% [X % Y mm] %% start calculation -matRad_calcStudy(examineStructures, multScen, param); +matRad_calcStudy(examineStructures, multScen, [], param); % exit; From a211ff15a7e346b6d9b65a7619496cd92a0d63b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 5 Mar 2018 13:56:08 +0100 Subject: [PATCH 07/13] Drafted beam wise shifts. --- matRad_calcDoseDirect.m | 2 +- matRad_sampling.m | 14 ++++---- tools/samplingAnalysis/matRad_calcStudy.m | 43 ++++++++++++++++++++--- 3 files changed, 46 insertions(+), 13 deletions(-) diff --git a/matRad_calcDoseDirect.m b/matRad_calcDoseDirect.m index d465dc855..ae9c6618a 100644 --- a/matRad_calcDoseDirect.m +++ b/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 diff --git a/matRad_sampling.m b/matRad_sampling.m index 09fc4fc13..074f385a0 100644 --- a/matRad_sampling.m +++ b/matRad_sampling.m @@ -99,13 +99,13 @@ matRad_dispToConsole(['matRad: Realizations variable will need: ' num2str(StorageInfo.bytes/1e9) ' GB \n'],param,'info'); % check if parallel toolbox is installed and license can be checked out -% try -% ver('distcomp') -% FlagParallToolBoxLicensed = license('test','Distrib_Computing_Toolbox'); -% catch -% FlagParallToolBoxLicensed = false; -% end -FlagParallToolBoxLicensed = false; +try + ver('distcomp') + FlagParallToolBoxLicensed = license('test','Distrib_Computing_Toolbox'); +catch + FlagParallToolBoxLicensed = false; +end + %% calculate nominal scenario nomScenTimer = tic; resultGUInomScen = matRad_calcDoseDirect(ct,stf,plnNominal,cst,w,param); diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index ef45a9866..0ed6de6c3 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -32,19 +32,39 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function [treatmentSimulation, scenContainer] = mergeSplittedScen(scenContainerParts) +function [treatmentSimulation, scenContainer, nominalScenario] = mergeSplittedScen(scenContainerParts, nomScenParts) + % merge nominal scenario + nominalScenario = nomScenParts{1}; + for k = 2:numel(nomScenParts) + nominalScenario.cumulateDose(nomScenParts{k}.dose); + end scenContainer = scenContainerParts{1}; for j = 1:numel(scenContainer) for k = 2:numel(scenContainerParts) - scenContainer{j}.cumulateDose(scenContainerParts{k}{j}); + scenContainer{j}.cumulateDose(scenContainerParts{k}{j}.dose); end end - treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nomScen, ct, cst, pln, pln.numOfFractions); + treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nominalScenario, ct, cst, pln, pln.numOfFractions); for j = 1:numel(scenContainer) treatmentSimulation.initNewScen(scenContainer{j}) end end +function stf = copyWeightsToStf(stf, w) + if sum([stf.totalNumOfBixels]) ~= numel(w) + error('weighting does not match steering information') + end + counter = 0; + for ii = 1:numel(stf) + for j = 1:stf(ii).numOfRays + for k = 1:stf(ii).numOfBixelsPerRay(j) + counter = counter + 1; + stf(ii).ray(j).weight(k) = w(counter); + end + end + end +end + if exist('param','var') if ~isfield(param,'logLevel') param.logLevel = 4; @@ -53,7 +73,6 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) param.logLevel = 4; end -% if ~isfield(param,'outputPath') param.outputPath = mfilename('fullpath'); end @@ -72,6 +91,13 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) end %% load DICOM imported patient +cst = []; +ct = []; +pln = []; +dij = []; +resultGUI = []; +stf = []; + if exist('matPatientPath', 'var') && ~isempty(matPatientPath) && exist('matPatientPath','file') == 2 load(matPatientPath) else @@ -84,6 +110,8 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) end end +stf = copyWeightsToStf(stf, resultGUI.w); + % check if nominal workspace is complete if ~(exist('ct','var') && exist('cst','var') && exist('stf','var') && exist('pln','var') && exist('resultGUI','var')) matRad_dispToConsole('Nominal workspace for sampling is incomplete.\n',param,'error'); @@ -120,10 +148,12 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) [treatmentSimulation, scenContainer, pln, resultGUInomScen, nomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,multScen,param); elseif strcmp(multScen.deepestCorrelationBreak, 'beam') scenContainerParts = cell(numel(stf),1); + nomScenParts = scenContainerParts; for i = 1:numel(stf) stfPart = stf(i); plnPart = pln; plnPart.isoCenter = pln.isoCenter(i,:); + plnPart.numOfBeams = 1; if strcmp(multScen.rangeShiftCorrelationBreak, 'beam') && strcmp(multScen.isoShiftCorrelationBreak, 'beam') multScen = multScen.matRad_recreateInstance; elseif strcmp(multScen.isoShiftCorrelationBreak, 'beam') @@ -134,10 +164,12 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) erorr('Incosistency'); end - [~, scenContainerPart, pln, ~, ~] = matRad_sampling(ct,stfPart,cst,plnPart,resultGUI.w,structSel,multScen,param); + [~, scenContainerPart, ~, ~, nomScenPart] = matRad_sampling(ct, stfPart, cst, plnPart, [],structSel, multScen, param); + nomScenParts{i} = nomScenPart; scenContainerParts{i} = scenContainerPart; end + [treatmentSimulation, scenContainer, nomScen] = mergeSplittedScen(scenContainerParts, nomScenParts); end param.computationTime = toc; @@ -174,3 +206,4 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) % else % system(executeLatex); % end +end From 4b6f462966deb29a28745eb9f7cdc1e3810f4c3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 5 Mar 2018 15:53:22 +0100 Subject: [PATCH 08/13] Reduce dependency on pln (use stf instead), bugfix in rot matrix --- matRad_calcParticleDose.m | 6 +++--- tools/samplingAnalysis/matRad_calcStudy.m | 9 +++++++-- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/matRad_calcParticleDose.m b/matRad_calcParticleDose.m index fae07aee3..692796dfc 100644 --- a/matRad_calcParticleDose.m +++ b/matRad_calcParticleDose.m @@ -58,7 +58,7 @@ end % meta information for dij -dij.numOfBeams = pln.numOfBeams; +dij.numOfBeams = numel(stf); dij.numOfVoxels = pln.numOfVoxels; dij.resolution = ct.resolution; dij.dimensions = pln.voxelDimensions; @@ -263,7 +263,7 @@ % compute SSDs stf = matRad_computeSSD(stf,ct,ctScen); - 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'); @@ -287,7 +287,7 @@ % transformation of the coordinate system need double transpose % rotation around Z axis (gantry) - rotMat_system_T = matRad_getRotationMatrix(pln.gantryAngles(i),pln.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; diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index 0ed6de6c3..731921065 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -33,7 +33,10 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [treatmentSimulation, scenContainer, nominalScenario] = mergeSplittedScen(scenContainerParts, nomScenParts) - % merge nominal scenario + % Merge nominal scenario by cumulating on first part. Remember that all + % scenario class is subclass of abstract handle class, therefore + % nomScenParts{1} changes. + nominalScenario = nomScenParts{1}; for k = 2:numel(nomScenParts) nominalScenario.cumulateDose(nomScenParts{k}.dose); @@ -148,10 +151,12 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) [treatmentSimulation, scenContainer, pln, resultGUInomScen, nomScen] = matRad_sampling(ct,stf,cst,pln,resultGUI.w,structSel,multScen,param); elseif strcmp(multScen.deepestCorrelationBreak, 'beam') scenContainerParts = cell(numel(stf),1); - nomScenParts = scenContainerParts; + nomScenParts = cell(numel(stf),1); for i = 1:numel(stf) stfPart = stf(i); plnPart = pln; + plnPart.couchAngles = pln.couchAngles(i); + plnPart.gantryAngles = pln.gantryAngles(i); plnPart.isoCenter = pln.isoCenter(i,:); plnPart.numOfBeams = 1; if strcmp(multScen.rangeShiftCorrelationBreak, 'beam') && strcmp(multScen.isoShiftCorrelationBreak, 'beam') From feb8eda9b4898de61da307bd1d4f13f010f42a19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Rapahel=20Mu=CC=88ller?= Date: Tue, 6 Mar 2018 12:50:26 +0100 Subject: [PATCH 09/13] Clear up --- tools/samplingAnalysis/MatRadScenario.m | 4 +-- tools/samplingAnalysis/MatRadSimulation.m | 1 - tools/samplingAnalysis/matRad_calcStudy.m | 30 ++++++++++++++++------- 3 files changed, 23 insertions(+), 12 deletions(-) diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 43ce49663..9333ccc97 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -84,8 +84,8 @@ function plotQiDVH(obj) end end - function plotDoseSlice(obj, ax, ct, plane, slice,legendOn) - matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice,[],[],jet,[],[],[],[],[],legendOn); + function plotDoseSlice(obj, ax, ct, plane, slice, doseWindow, legendOn) + matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice,[],[],jet,[],doseWindow,[],[],[],legendOn); end function singleStructDVH = getSingleStructDVH(obj, voi) diff --git a/tools/samplingAnalysis/MatRadSimulation.m b/tools/samplingAnalysis/MatRadSimulation.m index 3e9bca1dd..e82bd04f8 100644 --- a/tools/samplingAnalysis/MatRadSimulation.m +++ b/tools/samplingAnalysis/MatRadSimulation.m @@ -180,7 +180,6 @@ function plotDVHband(obj, listOfStruct) if ~obj.statisticsComputed warning('Please run statistical analysis first.'); else - for i = 1:numel(listOfStruct) for j = 1:numel(obj.dvhStatistics) if strcmp(obj.nominalScenario.dvh(j).name, listOfStruct{i}) diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index 731921065..0e574b3ec 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -32,7 +32,7 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -function [treatmentSimulation, scenContainer, nominalScenario] = mergeSplittedScen(scenContainerParts, nomScenParts) +function [treatmentSimulation, nominalScenario] = mergeSplittedScen(scenContainerParts, nomScenParts) % Merge nominal scenario by cumulating on first part. Remember that all % scenario class is subclass of abstract handle class, therefore % nomScenParts{1} changes. @@ -41,15 +41,15 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) for k = 2:numel(nomScenParts) nominalScenario.cumulateDose(nomScenParts{k}.dose); end - scenContainer = scenContainerParts{1}; - for j = 1:numel(scenContainer) + sc = scenContainerParts{1}; + for j = 1:numel(sc) for k = 2:numel(scenContainerParts) - scenContainer{j}.cumulateDose(scenContainerParts{k}{j}.dose); + sc{j}.cumulateDose(scenContainerParts{k}{j}.dose); end end treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nominalScenario, ct, cst, pln, pln.numOfFractions); - for j = 1:numel(scenContainer) - treatmentSimulation.initNewScen(scenContainer{j}) + for j = 1:numel(sc) + treatmentSimulation.initNewScen(sc{j}) end end @@ -104,9 +104,19 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) if exist('matPatientPath', 'var') && ~isempty(matPatientPath) && exist('matPatientPath','file') == 2 load(matPatientPath) else + % search for .mat files listOfMat = dir('*.mat'); + % exclude all resultSampling_*.mat + listOfMat = {listOfMat(:).name}; + deleteIx = false(1,numel(listOfMat)); + for i = 1:numel(listOfMat) + if strncmp(listOfMat{i}, 'resultSampling', 14) + deleteIx(i) = true; + end + end + listOfMat(deleteIx) = []; if numel(listOfMat) == 1 - load(listOfMat.name); + load(listOfMat{1}); else matRad_dispToConsole('Ambigous set of .mat files in the current folder (i.e. more than one possible patient or already results available).\n',param,'error'); return @@ -174,12 +184,14 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) scenContainerParts{i} = scenContainerPart; end - [treatmentSimulation, scenContainer, nomScen] = mergeSplittedScen(scenContainerParts, nomScenParts); + [treatmentSimulation, nomScen] = mergeSplittedScen(scenContainerParts, nomScenParts); end param.computationTime = toc; param.reportPath = fullfile('report','data'); -filename = 'resultSampling'; +t = datetime('now','Format','yyyy-MM-dd''T''HHmmss'); +filename = ['resultSampling_', char(t)]; +clear scenContainer scenContainerPart scenContainerParts % redundant vars save(filename, '-v7.3'); %% perform analysis From fd14bda27e74a7d3e93c0c5c76200659018cc128 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Rapahel=20Mu=CC=88ller?= Date: Mon, 12 Mar 2018 13:28:56 +0100 Subject: [PATCH 10/13] Store information when cumulating dose --- tools/samplingAnalysis/MatRadScenario.m | 14 +++++++++++++- tools/samplingAnalysis/matRad_calcStudy.m | 13 ++++++++++--- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 9333ccc97..550ca664a 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -13,6 +13,7 @@ shift; relRangeShift; absRangeShift; + isoShift; ctShiftIdentifier; radiationQuantity; doseCubeDim; @@ -61,8 +62,12 @@ addlistener(obj,'dose','PostSet',@obj.handleChangeOfDose); end % eof constructor - function obj = cumulateDose(obj, dose) + function obj = cumulateDose(obj, dose, ctShiftIdentifier, shift, absRangeShift, relRangeShift) obj.dose = obj.dose + dose; + obj.ctShiftIdentifier = [obj.ctShiftIdentifier ctShiftIdentifier]; + obj.shift = [obj.shift; shift]; + obj.absRangeShift = [obj.absRangeShift absRangeShift]; + obj.relRangeShift = [obj.relRangeShift relRangeShift]; end function obj = calcQiDVH(obj, cst, pln, dvhType, doseGrid, refGy, refVol) @@ -140,6 +145,13 @@ function plotDoseSlice(obj, ax, ct, plane, slice, doseWindow, legendOn) end end + function shift = set.shift(obj, v) + if isempty(v) + obj.shift = [0 0 0]; + else + obj.shift = v; + end + end end % eof methods diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index 0e574b3ec..2fd7618ae 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -39,12 +39,16 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) nominalScenario = nomScenParts{1}; for k = 2:numel(nomScenParts) - nominalScenario.cumulateDose(nomScenParts{k}.dose); + nominalScenario.cumulateDose(nomScenParts{k}.dose, 1, [0 0 0], 0, 0); end sc = scenContainerParts{1}; for j = 1:numel(sc) for k = 2:numel(scenContainerParts) - sc{j}.cumulateDose(scenContainerParts{k}{j}.dose); + sc{j}.cumulateDose(scenContainerParts{k}{j}.dose, ... + 1, ... % ct shift identifier + scenContainerParts{k}{j}.shift, ... + scenContainerParts{k}{j}.absRangeShift, ... + scenContainerParts{k}{j}.relRangeShift); end end treatmentSimulation = MatRadSimulation(pln.bioParam.quantityVis, nominalScenario, ct, cst, pln, pln.numOfFractions); @@ -79,6 +83,9 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) if ~isfield(param,'outputPath') param.outputPath = mfilename('fullpath'); end +if ~isfield(param,'fileSuffix') + param.fileSuffix = ''; +end if ~isfield(param, 'percentiles') param.percentiles = [0.01 0.05 0.125 0.25 0.5 0.75 0.875 0.95 0.99]; end @@ -190,7 +197,7 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) param.reportPath = fullfile('report','data'); t = datetime('now','Format','yyyy-MM-dd''T''HHmmss'); -filename = ['resultSampling_', char(t)]; +filename = ['resultSampling_', char(t), param.fileSuffix]; clear scenContainer scenContainerPart scenContainerParts % redundant vars save(filename, '-v7.3'); From 93ca8eb124398e22f80218b55b921faa4913a116 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 19 Mar 2018 14:31:35 +0100 Subject: [PATCH 11/13] Add daemon to continously evaluate patients. Minor bugfixes. --- matRad_calcQualityIndicators.m | 2 +- matRad_simulationAutoRunDaemon.m | 102 ++++++++++++++++++++++ tools/samplingAnalysis/MatRadSimulation.m | 44 +++++++--- 3 files changed, 133 insertions(+), 15 deletions(-) create mode 100644 matRad_simulationAutoRunDaemon.m diff --git a/matRad_calcQualityIndicators.m b/matRad_calcQualityIndicators.m index 503d4c85b..34a54926f 100644 --- a/matRad_calcQualityIndicators.m +++ b/matRad_calcQualityIndicators.m @@ -137,7 +137,7 @@ for j = 1:numel(listOfFields) qi(i).(listOfFields{j}) = NaN; end - qi(i).VOIname = cst{i,2}; + qi(i).name = cst{i,2}; end end diff --git a/matRad_simulationAutoRunDaemon.m b/matRad_simulationAutoRunDaemon.m new file mode 100644 index 000000000..002637a50 --- /dev/null +++ b/matRad_simulationAutoRunDaemon.m @@ -0,0 +1,102 @@ +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% matrad_simulationAutoRunDaemon enables continuos evaluation of simulation +% +% call +% matrad_simulationAutoRunDaemon +% +% usage +% place this file in the parent folder of your simulations and create a +% directory 'queue'. Place patients .mat in a subfolder as well as a +% setupStudy configuration file. This daemon will continue to look for +% new simulation folders until the script is terminated by the user. +% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 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/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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +directoryOfQueue = 'queue'; +marker = {'_RUNNING', '_FINISHED', '_ERROR'}; +numWorker = 4; + +currPath = pwd; +computeQueue = true; + +while computeQueue + listOfPatients = dir(directoryOfQueue); + listOfPatients(1:2) = []; % remove '.' and '..' + + % if no patient is in queue wait a bit longer + if numel(listOfPatients) == 0 + fprintf('No patient found. Waiting ... \n'); + pause(300); + continue + end + + % reduce list to non-running and non-finished + i = 1; + while true + if (~isempty(strfind(listOfPatients(i).name, marker{1})) ... + || ~isempty(strfind(listOfPatients(i).name, marker{2})) ... + || ~isempty(strfind(listOfPatients(i).name, marker{3}))) + listOfPatients(i) = []; + i = 1; % jump to beginning of list + else + i = i + 1; % step forward + end + if i >= numel(listOfPatients) + 1 + break + end + end + + % if all patients are running or finished + if numel(listOfPatients) == 0 + fprintf('All simulations finished or already running. Waiting ... \n'); + pause(300); + continue + end + + % pick a patient randomly + currPatientIx = randi(numel(listOfPatients)); + currPatientPath = fullfile(listOfPatients(currPatientIx).folder, listOfPatients(currPatientIx).name); + tempPatientPath = [currPatientPath marker{1}]; + movefile(currPatientPath, tempPatientPath); + + % get or start parpool + p = gcp('nocreate'); + if numel(p) == 0 + try + parpool(numWorker); + catch + warning('Could not start parpool.'); + end + end + + cd(tempPatientPath) + try + setupStudy + % after completion change directory back and rename + cd(currPath); + finishedPatientPath = [currPatientPath marker{2}]; + movefile(tempPatientPath, finishedPatientPath); + catch + warning('Error during computation of simulation.'); + cd(currPath); + finishedPatientPath = [currPatientPath '_ERROR']; + movefile(tempPatientPath, finishedPatientPath); + end + + % safety wait to not overdo it + pause(10) +end diff --git a/tools/samplingAnalysis/MatRadSimulation.m b/tools/samplingAnalysis/MatRadSimulation.m index e82bd04f8..b2f43e975 100644 --- a/tools/samplingAnalysis/MatRadSimulation.m +++ b/tools/samplingAnalysis/MatRadSimulation.m @@ -138,6 +138,13 @@ obj.statisticsComputed = true; end + function obj = runReducedAnalysis(obj) + obj.fillDoseContainer(); + obj.computeMeanStdCube(); + + obj.statisticsComputed = true; + end + function weights_unnormalised = get.weights_unnormalised(obj) weights_unnormalised = NaN * ones(obj.numOfScen, 1); for i = 1:obj.numOfScen @@ -193,7 +200,16 @@ function plotDVHband(obj, listOfStruct) end end end - end + + % consistency checker / repair functions + function checkRepairCst(obj) + fprintf('Trying to resolve utf8 character encoding. Recheck structure names.\n'); + for i = 1:size(obj.cst,1) + obj.cst{i,2} = regexprep(obj.cst{i,2}, '[^a-zA-Z0-9]', ' '); + end + end + + end methods (Static) function S = wMean(X,w) @@ -224,6 +240,19 @@ function computeAllDvhQi(obj, refVol, refGy) obj.scenContainer{i}.calcQiDVH(obj.cst, obj.pln, 'cum', doseGrid, refGy, refVol); end obj.dvhDoseGrid = doseGrid; + end + + function computeMeanStdCube(obj) + meanCube = zeros(obj.doseCubeDim); + stdCube = zeros(obj.doseCubeDim); + + ix = obj.subIx; + + meanCube(ix) = (sum(obj.doseContainer * diag(obj.weights),2)); + stdCube(ix) = std(obj.doseContainer, obj.weights,2); + + obj.meanDoseScenario = MatRadScenario(meanCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); + obj.stdDoseScenario = MatRadScenario(stdCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); end function computeStructureWiseDvhQi(obj) @@ -309,19 +338,6 @@ function computeStructureWiseStatistics(obj) end end - function computeMeanStdCube(obj) - meanCube = zeros(obj.doseCubeDim); - stdCube = zeros(obj.doseCubeDim); - - ix = obj.subIx; - - meanCube(ix) = (sum(obj.doseContainer * diag(obj.weights),2)); - stdCube(ix) = std(obj.doseContainer, obj.weights,2); - - obj.meanDoseScenario = MatRadScenario(meanCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); - obj.stdDoseScenario = MatRadScenario(stdCube, obj.subIx, obj.radiationQuantity, NaN, NaN, NaN, NaN, obj.doseCubeDim, NaN); - end - function computeGammaCube(obj) gammaCube = zeros(obj.doseCubeDim); gammaCube = matRad_gammaIndex(obj.nominalScenario.dose,obj.meanDoseScenario.dose,[obj.ct.resolution.x obj.ct.resolution.y obj.ct.resolution.z],[obj.gammaDoseAgreement obj.gammaDistAgreement]); From 302cd380df7a1388858e0ce7726957c8743a0fe7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 26 Mar 2018 12:39:28 +0200 Subject: [PATCH 12/13] Add help text for shift combination, update bioModel generation --- tools/samplingAnalysis/matRad_calcStudy.m | 8 +++++--- tools/samplingAnalysis/setupStudy_template.m | 3 +++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/tools/samplingAnalysis/matRad_calcStudy.m b/tools/samplingAnalysis/matRad_calcStudy.m index 3fa0a913c..3c9c9c858 100644 --- a/tools/samplingAnalysis/matRad_calcStudy.m +++ b/tools/samplingAnalysis/matRad_calcStudy.m @@ -82,11 +82,13 @@ function matRad_calcStudy(structSel,multScen,matPatientPath,param) % calculate RBExDose if ~isfield(pln, 'bioParam') if strcmp(pln.radiationMode, 'protons') - pln.bioOptimization = 'constRBE_RBExD'; + pln.bioOptimization = 'RBExD'; + pln.model = 'constRBE'; elseif strcmp(pln.radiationMode, 'carbon') - pln.bioOptimization = 'LEM_RBExD'; + pln.bioOptimization = 'RBExD'; + pln.model = 'LEM'; end - pln.bioParam = matRad_bioModel(pln.radiationMode,pln.bioOptimization); + pln.bioParam = matRad_bioModel(pln.radiationMode, pln.bioOptimization, pln.model); end diff --git a/tools/samplingAnalysis/setupStudy_template.m b/tools/samplingAnalysis/setupStudy_template.m index a9fd7e742..541b06f4a 100644 --- a/tools/samplingAnalysis/setupStudy_template.m +++ b/tools/samplingAnalysis/setupStudy_template.m @@ -9,6 +9,9 @@ multScen.shiftGenType = 'equidistant'; % equidistant: equidistant shifts, sampled: sample shifts from normal distribution multScen.shiftCombType = 'individual'; % individual: no combination of shift scenarios; number of shift scenarios is sum(multScen.numOfShiftScen) % permuted: create every possible shift combination; number of shift scenarios is 8,27,64 ... + % combined: number of shifts in each directions must be the same. only use when shifts are sampled + + % b) define range error scenarios multScen.numOfRangeShiftScen = 30; % number of absolute and/or relative range scnearios. % if absolute and relative range scenarios are defined then multScen.rangeCombType defines the resulting number of range scenarios From 8841c0834d5316895944a896a9864df31732f594 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucas-Raphael=20M=C3=BCller?= Date: Mon, 26 Mar 2018 14:01:38 +0200 Subject: [PATCH 13/13] Add colormap label --- tools/samplingAnalysis/MatRadScenario.m | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tools/samplingAnalysis/MatRadScenario.m b/tools/samplingAnalysis/MatRadScenario.m index 550ca664a..55ccee6e5 100644 --- a/tools/samplingAnalysis/MatRadScenario.m +++ b/tools/samplingAnalysis/MatRadScenario.m @@ -90,7 +90,14 @@ function plotQiDVH(obj) end function plotDoseSlice(obj, ax, ct, plane, slice, doseWindow, legendOn) - matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice,[],[],jet,[],doseWindow,[],[],[],legendOn); + if ~exist('doseWindow', 'var') || isempty(doseWindow) + doseWindow = [0 max(obj.dose(:))]; + end + if ~exist('legendOn', 'var') || isempty(legendOn) + legendOn = false; + end + colorMapLabel = obj.radiationQuantity; + matRad_plotSliceWrapper(ax,ct,obj.cst,1,obj.dose,plane,slice,[],[],jet,[],doseWindow,[],[],colorMapLabel,legendOn); end function singleStructDVH = getSingleStructDVH(obj, voi)