diff --git a/basedata/helium_Generic.mat b/basedata/helium_Generic.mat index ef14fda20..373555832 100644 Binary files a/basedata/helium_Generic.mat and b/basedata/helium_Generic.mat differ diff --git a/helium_Generic.mat b/helium_Generic.mat deleted file mode 100644 index 373555832..000000000 Binary files a/helium_Generic.mat and /dev/null differ diff --git a/matRadGUI.m b/matRadGUI.m index 2d6fec8c2..f01160ce1 100644 --- a/matRadGUI.m +++ b/matRadGUI.m @@ -191,6 +191,7 @@ end set(handles.popUpMachine,'String',handles.Machines); +%TODO: replace with class crawling multScenDummy = matRad_multScen([],'nomScen'); set(handles.popupmenuScenGen,'String',multScenDummy.AvailableScenCreationTYPE); @@ -310,9 +311,9 @@ function initViewSliceSlider(handles) setPln(handles); end -catch +catch ME handles.State = 0; - handles = showError(handles,'GUI OpeningFunc: Could not set or get pln'); + handles = showError(handles,sprintf('GUI OpeningFunc: Could not set or get pln: %s',ME.message)); end % check for dij structure @@ -3668,8 +3669,8 @@ function sliderOpacity_CreateFcn(hObject, eventdata, handles) cursorText{end+1,1} = ['Cube Index: ' mat2str(cubeIx)]; %Space Coordinates coords = zeros(1,3); - coords(1) = cubePos(2)*ct.resolution.y; - coords(2) = cubePos(1)*ct.resolution.x; + coords(1) = cubePos(2)*ct.resolution.x; + coords(2) = cubePos(1)*ct.resolution.y; coords(3) = cubePos(3)*ct.resolution.z; cursorText{end+1,1} = ['Space Coordinates: ' mat2str(coords,5) ' mm']; @@ -4301,7 +4302,7 @@ function popupmenuScenGen_Callback(hObject, eventdata, handles) % hObject handle to popupmenuScenGen (see GCBO) contents = cellstr(get(hObject,'String')); -if handles.State > 1 +if handles.State >= 1 ct = evalin('base','ct'); pln = evalin('base','pln'); pln.scenGenType = contents{get(hObject,'Value')}; diff --git a/matRad_calcDoseInitBeam.m b/matRad_calcDoseInitBeam.m index a191f041b..cd4b3677e 100644 --- a/matRad_calcDoseInitBeam.m +++ b/matRad_calcDoseInitBeam.m @@ -52,13 +52,10 @@ % interpolate radiological depth cube to dose grid resolution radDepthVdoseGrid = matRad_interpRadDepth(ct,VctGrid,VdoseGrid,dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z,radDepthVctGrid); - + if exist('radDepthsMat', 'var') - for ctScen = 1:ct.numOfCtScen - % interpolate radiological depth cube used for fine sampling to dose grid resolution - radDepthsMat{ctScen} = matRad_interp3(dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z, radDepthsMat{1}, ... - dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'); - end + % interpolate radiological depth cube used for fine sampling to dose grid resolution + radDepthsMat = cellfun(@(radDepthCube) matRad_interp3(dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z,radDepthCube,dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'),radDepthsMat,'UniformOutput',false); end % limit rotated coordinates to positions where ray tracing is availabe diff --git a/matRad_calcParticleDose.m b/matRad_calcParticleDose.m index 615961dfa..904c890ac 100644 --- a/matRad_calcParticleDose.m +++ b/matRad_calcParticleDose.m @@ -55,19 +55,14 @@ for ctScen = 1:pln.multScen.numOfCtScen for shiftScen = 1:pln.multScen.totNumShiftScen - for rangeShiftScen = 1:pln.multScen.totNumRangeScen - + for rangeShiftScen = 1:pln.multScen.totNumRangeScen if pln.multScen.scenMask(ctScen,shiftScen,rangeShiftScen) dij.mAlphaDose{ctScen,shiftScen,rangeShiftScen} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1); dij.mSqrtBetaDose{ctScen,shiftScen,rangeShiftScen} = spalloc(dij.doseGrid.numOfVoxels,numOfColumnsDij,1); - end - - end - - end - - end - + end + end + end + end end if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'calcLET') @@ -181,7 +176,6 @@ end % end is LEM model end - % lateral cutoff for raytracing and geo calculations if ~isfield(pln,'propDoseCalc') || ~isfield(pln.propDoseCalc,'geometricCutOff') effectiveLateralCutoff = matRad_cfg.propDoseCalc.defaultGeometricCutOff; @@ -191,8 +185,6 @@ pln.propDoseCalc.lateralCutOff = matRad_cfg.propDoseCalc.defaultLateralCutOff; end - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %loop over all shift scenarios %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -229,17 +221,16 @@ if ~isempty(stf(i).ray(j).energy) % find index of maximum used energy (round to keV for numerical reasons - energyIx = max(round2(stf(i).ray(j).energy,4)) == round2([machine.data.energy],4); maxLateralCutoffDoseCalc = max(machine.data(energyIx).LatCutOff.CutOff); - + % calculate initial sigma for all bixel on current ray sigmaIniRay = matRad_calcSigmaIni(machine.data,stf(i).ray(j),stf(i).ray(j).SSD); - if strcmp(pbCalcMode, 'fineSampling') - % Ray tracing for beam i and ray j + % Ray tracing for beam i and ray j with explicit + % lateral distances for fine sampling [ix,~,~,~,latDistsX,latDistsZ] = matRad_calcGeoDists(rot_coordsVdoseGrid, ... stf(i).sourcePoint_bev, ... stf(i).ray(j).targetPoint_bev, ... @@ -263,23 +254,22 @@ end end else - % Ray tracing for beam i and ray j + % Ray tracing for beam i and ray j without explicitly + % obtaining lateral distances [ix,currRadialDist_sq,~,~,~,~] = matRad_calcGeoDists(rot_coordsVdoseGrid, ... stf(i).sourcePoint_bev, ... stf(i).ray(j).targetPoint_bev, ... machine.meta.SAD, ... find(~isnan(radDepthVdoseGrid{1})), ... maxLateralCutoffDoseCalc); - - radDepths = radDepthVdoseGrid{1}(ix); end % just use tissue classes of voxels found by ray tracer if pln.bioParam.bioOpt vTissueIndex_j = vTissueIndex(ix,:); - end + end - for k = 1:stf(i).numOfBixelsPerRay(j) % loop over all bixels per ray + for k = 1:stf(i).numOfBixelsPerRay(j) % loop over all bixels per ray counter = counter + 1; bixelsPerBeam = bixelsPerBeam + 1; @@ -289,14 +279,13 @@ if mod(bixelsPerBeam,max(1,round(stf(i).totalNumOfBixels/200))) == 0 matRad_progress(bixelsPerBeam/max(1,round(stf(i).totalNumOfBixels/200)),... floor(stf(i).totalNumOfBixels/max(1,round(stf(i).totalNumOfBixels/200)))); - end - - % update waitbar only 100 times if it is not closed - if mod(counter,round(dij.totalNumOfBixels/100)) == 0 && ishandle(figureWait) - waitbar(counter/dij.totalNumOfBixels,figureWait); - end + end end - + % update waitbar only 100 times if it is not closed + if mod(counter,round(dij.totalNumOfBixels/100)) == 0 && ishandle(figureWait) + waitbar(counter/dij.totalNumOfBixels,figureWait); + end + % remember beam and bixel number if ~calcDoseDirect dij.beamNum(counter) = i; @@ -323,7 +312,7 @@ dij.maxMU(counter,1) = maxMU; dij.numParticlesPerMU(counter,1) = numParticlesPerMU; end - + % find energy index in base data energyIx = find(round2(stf(i).ray(j).energy(k),4) == round2([machine.data.energy],4)); @@ -343,83 +332,96 @@ % that is within the limits of the base data set (-> machine.data(i).dephts). By this means, we only allow % interpolations in matRad_calcParticleDoseBixel() and avoid extrapolations. offsetRadDepth = machine.data(energyIx).offset - (stf(i).ray(j).rangeShifter(k).eqThickness + dR); - - %Fine Sampling coordinate projections - if strcmp(pbCalcMode, 'fineSampling') - % calculate projected coordinates for fine sampling of - % each beamlet + + % calculate projected coordinates for fine sampling of + % each beamlet + if strcmp(pbCalcMode, 'fineSampling') projCoords = matRad_projectOnComponents(VdoseGrid(ix), size(radDepthsMat{1}), stf(i).sourcePoint_bev,... stf(i).ray(j).targetPoint_bev, stf(i).isoCenter,... [dij.doseGrid.resolution.x dij.doseGrid.resolution.y dij.doseGrid.resolution.z],... -posX(:,k), -posZ(:,k), rotMat_system_T); end + % We do now loop over scenarios that alter voxel + % values, e.g. range scenarios or ct phases, as we can + % vectorize computations more efficiently than when + % making this an outer loop for ctScen = 1:pln.multScen.numOfCtScen + if any(any(pln.multScen.scenMask(ctScen,:,:))) %We don't need it if no scenario for this ct scenario is relevant + % precomputations for fine-sampling + if strcmp(pbCalcMode, 'fineSampling') + % compute radial distances relative to pencil beam + % component + currRadialDist_sq = reshape(bsxfun(@plus,latDistsX,posX(:,k)'),[],1,numOfSub(k)).^2 + reshape(bsxfun(@plus,latDistsZ,posZ(:,k)'),[],1,numOfSub(k)).^2; + + % interpolate radiological depths at projected + % coordinates + radDepths = interp3(radDepthsMat{ctScen},projCoords(:,1,:)./dij.doseGrid.resolution.x,... + projCoords(:,2,:)./dij.doseGrid.resolution.y,projCoords(:,3,:)./dij.doseGrid.resolution.z,'nearest'); + else + radDepths = radDepthVdoseGrid{ctScen}(ix); + end + end for rangeShiftScen = 1:pln.multScen.totNumRangeScen - if pln.multScen.scenMask(ctScen,shiftScen,rangeShiftScen) - - %Fine Sampling radiological depths - if strcmp(pbCalcMode, 'fineSampling') - - % interpolate radiological depths at projected - % coordinates - radDepths = interp3(radDepthsMat{ctScen},projCoords(:,1,:)./dij.doseGrid.resolution.x,... - projCoords(:,2,:)./dij.doseGrid.resolution.y,projCoords(:,3,:)./dij.doseGrid.resolution.z,'nearest'); - - % compute radial distances relative to pencil beam - % component - currRadialDist_sq = reshape(bsxfun(@plus,latDistsX,posX(:,k)'),[],1,numOfSub(k)).^2 + reshape(bsxfun(@plus,latDistsZ,posZ(:,k)'),[],1,numOfSub(k)).^2; - else - radDepths = radDepthVdoseGrid{ctScen}(ix); - end + if pln.multScen.scenMask(ctScen,shiftScen,rangeShiftScen) % manipulate radDepthCube for range scenarios if pln.multScen.relRangeShift(rangeShiftScen) ~= 0 || pln.multScen.absRangeShift(rangeShiftScen) ~= 0 - radDepths = radDepths + ... - radDepths*pln.multScen.relRangeShift(rangeShiftScen) +... % rel range shift + currRadDepths = radDepths * (1+pln.multScen.relRangeShift(rangeShiftScen)) +... % rel range shift pln.multScen.absRangeShift(rangeShiftScen); % absolute range shift - radDepths(radDepths < 0) = 0; + currRadDepths(currRadDepths < 0) = 0; + else + currRadDepths = radDepths; end - - - + % find depth depended lateral cut off if cutOffLevel >= 1 - currIx = radDepths <= machine.data(energyIx).depths(end) + offsetRadDepth; + currIx = currRadDepths <= machine.data(energyIx).depths(end) + offsetRadDepth; elseif cutOffLevel < 1 && cutOffLevel > 0 % perform rough 2D clipping - currIx = radDepths <= machine.data(energyIx).depths(end) + offsetRadDepth & ... + currIx = currRadDepths <= machine.data(energyIx).depths(end) + offsetRadDepth & ... currRadialDist_sq <= max(machine.data(energyIx).LatCutOff.CutOff.^2); % peform fine 2D clipping if length(machine.data(energyIx).LatCutOff.CutOff) > 1 currIx(currIx) = matRad_interp1((machine.data(energyIx).LatCutOff.depths + offsetRadDepth)',... - (machine.data(energyIx).LatCutOff.CutOff.^2)', radDepths(currIx)) >= currRadialDist_sq(currIx); + (machine.data(energyIx).LatCutOff.CutOff.^2)', currRadDepths(currIx)) >= currRadialDist_sq(currIx); end else matRad_cfg.dispError('Lateral Cut-Off must be a value between 0 and 1!') end - + % empty bixels may happen during recalculation of error % scenarios -> skip to next bixel - if ~any(currIx) + if ~any(currIx) + %Create empty container entries for + %this bixel doseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(dij.doseGrid.numOfVoxels,1); + if isfield(dij,'mLETDose') + letDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(dij.doseGrid.numOfVoxels,1); + end + if pln.bioParam.bioOpt + alphaDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(VdoseGrid(ix(currIx)),1,bixelAlpha.*bixelDose,dij.doseGrid.numOfVoxels,1); + betaDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(VdoseGrid(ix(currIx)),1,sqrt(bixelBeta).*bixelDose,dij.doseGrid.numOfVoxels,1); + end + + %skip bixel continue; end if pln.propDoseCalc.airOffsetCorrection - currRadDepths = radDepths(currIx) + stf(i).ray(j).rangeShifter(k).eqThickness + dR; + currRadDepths(currIx) = currRadDepths(currIx) + stf(i).ray(j).rangeShifter(k).eqThickness + dR; %sanity check due to negative corrections currRadDepths(currRadDepths < 0) = 0; else - currRadDepths = radDepths(currIx) + stf(i).ray(j).rangeShifter(k).eqThickness; + currRadDepths(currIx) = currRadDepths(currIx) + stf(i).ray(j).rangeShifter(k).eqThickness; end - - % calculate initial focus sigma + + % select correct initial focus sigma squared sigmaIni_sq = sigmaIniRay(k)^2; - + % consider range shifter for protons if applicable if stf(i).ray(j).rangeShifter(k).eqThickness > 0 && strcmp(pln.radiationMode,'protons') @@ -430,9 +432,8 @@ % add to initial sigma in quadrature sigmaIni_sq = sigmaIni_sq + sigmaRashi^2; - end - + if strcmp(pbCalcMode, 'fineSampling') % initialise empty dose array totalDose = zeros(size(currIx,1),1); @@ -447,7 +448,7 @@ for c = 1:numOfSub(k) tmpDose = zeros(size(currIx,1),1); bixelDose = finalWeight(c,k).*matRad_calcParticleDoseBixel(... - radDepths(currIx(:,:,c),1,c), ... + currRadDepths(currIx(:,:,c),1,c), ... currRadialDist_sq(currIx(:,:,c),:,c), ... sigmaSub(k)^2, ... machine.data(energyIx)); @@ -457,7 +458,7 @@ if isfield(dij,'mLETDose') tmpLET = zeros(size(currIx,1),1); - tmpLET(currIx(:,:,c)) = matRad_interp1(depths,machine.data(energyIx).LET,radDepths(currIx(:,:,c),1,c)); + tmpLET(currIx(:,:,c)) = matRad_interp1(depths,machine.data(energyIx).LET,currRadDepths(currIx(:,:,c),1,c)); totalLET = totalLET + tmpLET; end end @@ -467,38 +468,36 @@ letDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(VdoseGrid(ix),1,totalDose.*totalLET,dij.doseGrid.numOfVoxels,1); end else - % calculate particle dose for bixel k on ray j of beam i bixelDose = matRad_calcParticleDoseBixel(... - currRadDepths, ... + currRadDepths(currIx), ... currRadialDist_sq(currIx), ... sigmaIni_sq, ... machine.data(energyIx)); - % dij sampling is exluded for particles until we investigated the influence of voxel sampling for particles %relDoseThreshold = 0.02; % sample dose values beyond the relative dose %Type = 'dose'; - %[currIx,bixelDose] = matRad_DijSampling(currIx,bixelDose,radDepths(currIx),currRadialDist_sq(currIx),Type,relDoseThreshold); + %[currIx,bixelDose] = matRad_DijSampling(currIx,bixelDose,radDepths(currIx),radialDist_sq(currIx),Type,relDoseThreshold); - % save dose for every bixel in cell array + % Save dose for every bixel in cell array doseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(VdoseGrid(ix(currIx)),1,bixelDose,dij.doseGrid.numOfVoxels,1); - if isfield(dij,'mLETDose') % calculate particle LET for bixel k on ray j of beam i - depths = machine.data(energyIx).depths + machine.data(energyIx).offset; - bixelLET = matRad_interp1(depths,machine.data(energyIx).LET,currRadDepths); + depths = machine.data(energyIx).depths + machine.data(energyIx).offset; + bixelLET = matRad_interp1(depths,machine.data(energyIx).LET,currRadDepths(currIx)); bixelLET(isnan(bixelLET)) = 0; - % save LET for every bixel in cell array + + % Save LET for every bixel in cell array letDoseTmpContainer{mod(counter-1,numOfBixelsContainer)+1,ctScen,shiftScen,rangeShiftScen} = sparse(VdoseGrid(ix(currIx)),1,bixelLET.*bixelDose,dij.doseGrid.numOfVoxels,1); end end - + % save alpha_p and beta_p radiosensititvy parameter for every bixel in cell array if pln.bioParam.bioOpt - [bixelAlpha,bixelBeta] = pln.bioParam.calcLQParameter(currRadDepths,machine.data(energyIx),vTissueIndex_j(currIx,:),dij.ax(VdoseGrid(ix(currIx))),... + [bixelAlpha,bixelBeta] = pln.bioParam.calcLQParameter(currRadDepths(currIx),machine.data(energyIx),vTissueIndex_j(currIx,:),dij.ax(VdoseGrid(ix(currIx))),... dij.bx(VdoseGrid(ix(currIx))),... dij.abx(VdoseGrid(ix(currIx)))); @@ -517,8 +516,6 @@ end % end bixels per ray - - end end % end ray loop diff --git a/matRad_fluenceOptimization.m b/matRad_fluenceOptimization.m index 58edbf5c9..27132404d 100644 --- a/matRad_fluenceOptimization.m +++ b/matRad_fluenceOptimization.m @@ -97,11 +97,11 @@ wOnes = ones(dij.totalNumOfBixels,1); % calculate initial beam intensities wInit -matRad_cfg.dispInfo('Getting initial weights... '); +matRad_cfg.dispInfo('Estimating initial weights... '); if exist('wInit','var') %do nothing as wInit was passed to the function - matRad_cfg.dispInfo('use as given!\n'); -elseif strcmp(pln.bioParam.model,'constRBE') && strcmp(pln.radiationMode,'protons') + matRad_cfg.dispInfo('chosen provided wInit!\n'); +elseif strcmp(pln.bioParam.model,'constRBE') && strcmp(pln.radiationMode,'protons') % check if a constant RBE is defined - if not use 1.1 if ~isfield(dij,'RBE') dij.RBE = 1.1; @@ -126,7 +126,7 @@ for j = 1:size(cst{i,6},2) % check if prescribed doses are in a valid domain if any(cst{i,6}{j}.getDoseParameters() > 5) && isequal(cst{i,3},'TARGET') - matRad_cfg.dispError('Reference dose > 10 Gy[RBE] for target. Biological optimization outside the valid domain of the base data. Reduce dose prescription or use more fractions.\n'); + matRad_cfg.dispWarning('Reference dose > 10 Gy[RBE] for target. Biological optimization outside the valid domain of the base data. Reduce dose prescription or use more fractions.\n'); end end @@ -164,8 +164,8 @@ 4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*doseTmp(V))); wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(doseTmp(V))))* wOnes; end - -else + matRad_cfg.dispInfo('chosen weights adapted to biological dose calculation!\n'); +else doseTmp = dij.physicalDose{1}*wOnes; bixelWeight = (doseTarget)/mean(doseTmp(V)); wInit = wOnes * bixelWeight; diff --git a/matRad_generateStf.m b/matRad_generateStf.m index be0f10e37..ffff8d953 100644 --- a/matRad_generateStf.m +++ b/matRad_generateStf.m @@ -293,10 +293,12 @@ Counter = 0; + %Here we iterate through scenarios to check the required + %energies w.r.t lateral position. + %TODO: iterate over the linear scenario mask instead? for CtScen = 1:pln.multScen.numOfCtScen for ShiftScen = 1:pln.multScen.totNumShiftScen for RangeShiftScen = 1:pln.multScen.totNumRangeScen - if pln.multScen.scenMask(CtScen,ShiftScen,RangeShiftScen) Counter = Counter+1; @@ -322,8 +324,7 @@ %This captures the case that the first %relevant voxel is a target voxel targetEntry(Counter,1:length(entryIx)) = (radDepths(entryIx) + radDepths(entryIx+1)) ./ 2; - targetExit(Counter,1:length(exitIx)) = (radDepths(exitIx) + radDepths(exitIx+1)) ./ 2; - + targetExit(Counter,1:length(exitIx)) = (radDepths(exitIx) + radDepths(exitIx+1)) ./ 2; end end @@ -458,7 +459,7 @@ stf(i).numOfRays = size(stf(i).ray,2); % post processing for particle remove energy slices - if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'carbon') + if strcmp(stf(i).radiationMode,'protons') || strcmp(stf(i).radiationMode,'helium') || strcmp(stf(i).radiationMode,'carbon') % get minimum energy per field minEnergy = min([stf(i).ray.energy]); diff --git a/matRad_getIsoCenter.m b/matRad_getIsoCenter.m index 249f79f88..a62ce658b 100644 --- a/matRad_getIsoCenter.m +++ b/matRad_getIsoCenter.m @@ -39,13 +39,17 @@ V = []; %Check if any constraints/Objectives have been defined yet -noObjOrConst = all(cellfun(@isempty,cst(:,6))); +if size(cst,2) >= 6 + noObjOrConst = all(cellfun(@isempty,cst(:,6))); +else + noObjOrConst = true; +end % Save target indices in V variable. for i = 1:size(cst,1) % We only let a target contribute if it has an objective/constraint or % if we do not have specified objectives/constraints at all so far - if isequal(cst{i,3},'TARGET') && (~isempty(cst{i,6}) || noObjOrConst) + if isequal(cst{i,3},'TARGET') && (noObjOrConst || ~isempty(cst{i,6})) V = [V; cst{i,4}{1}]; end end diff --git a/matRad_multScen.m b/matRad_multScen.m index 255fa3fba..4ce1628da 100644 --- a/matRad_multScen.m +++ b/matRad_multScen.m @@ -1,764 +1,52 @@ -classdef matRad_multScen - % matRad_multScen - % This class creates all required biological model parameters according to - % a given radiation modatlity and a given bio model identifier string. - % - % constructor call - % pln.multScen = matRad_multScen(ct,TYPE,uctModel,combineScenarios); - % - % e.g. pln.multScen = matRad_multScen(ct,'nomScen'); - % - % input - % ct: ct cube - % TYPE: string to denote a scenario creation method - % 'nomScen' create only the nominal scenario - % 'wcScen' create worst case scenarios - % 'impScen' create important/grid scenarios - % 'rndScen' create random scenarios - % - % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % - % Copyright 2017 the matRad development team. - % - % This file is part of the matRad project. It is subject to the license - % terms in the LICENSE file found in the top-level directory of this - % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part - % of the matRad project, including this file, may be copied, modified, - % propagated, or distributed except according to the terms contained in the - % LICENSE file. - % - % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - %% properties - - % public properties which can be changed outside this class - properties - - TYPE; % denotes the sampling type which cen be one of the following - % 'nomScen' create only the nominal scenario - % 'wcScen' create worst case scenarios - % 'impScen' create important/grid scenarios - % 'rndScen' create random scenarios - - % Uncertainty Model - rangeRelSD = 3.5; % given in % - rangeAbsSD = 1; % given in [mm] - shiftSD = [2.25 2.25 2.25]; % given in [mm] - wcFactor = 1.5; % Multiplier to compute the worst case / maximum shifts - - - shiftGenType; % 'equidistant': equidistant shifts, 'sampled': sample shifts from normal distribution - 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 . - - 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 - - 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 - - includeNomScen; % boolean to determine if the nominal scenario should be included - - numOfShiftScen; % number of shifts in x y and z direction - numOfRangeShiftScen; % number of absolute and/or relative range scnearios. - end - - %These properties are accessible but not editable - properties (SetAccess = private) - - % ct scenarios - numOfCtScen = 1; % number of imported ct scenarios - - % shift scenarios - - shiftSize; % 3x1 vector to define maximal shift in [mm] % (e.g. abdominal cases 5mm otherwise 3mm) - - % b) define range error scenarios - - % if absolute and relative range scenarios are defined then rangeCombType defines the resulting number of range scenarios - maxAbsRangeShift; % maximum absolute over and undershoot in mm - maxRelRangeShift; % maximum relative over and undershoot in % - - - % these parameters will be filled according to the choosen scenario type - isoShift; - relRangeShift; - absRangeShift; - - totNumShiftScen; % total number of shift scenarios in x,y and z direction - totNumRangeScen; % total number of range and absolute range scenarios - totNumScen; % total number of samples - - scenForProb; % matrix for probability calculation - each row denotes one scenario - scenProb; % probability of each scenario stored in a vector - scenMask; - linearMask; - end - - % constant properties which are visible outside of matRad_multScen - properties(Constant = true) - - AvailableScenCreationTYPE = {'nomScen','wcScen','impScen','rndScen'}; - end - - % constant deafult properties which are only visible within matRad_multScen - properties (Constant = true, Access = private) - 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 - - % 'nomScen' default parameters for nominal scenario - numOfShiftScen_nomScen = [0 0 0]; % number of shifts in x y and z direction - %shiftSize_nomScen = [0 0 0]; % given in [mm] - shiftGenType_nomScen = 'equidistant'; % equidistant: equidistant shifts, - shiftCombType_nomScen = 'individual'; % individual: no combination of shift scenarios; - numOfRangeShiftScen_nomScen = 0 % number of absolute and/or relative range scnearios. - maxAbsRangeShift_nomScen = 0; % maximum absolute over and undershoot in mm - maxRelRangeShift_nomScen = 0; % maximum relative over and undershoot in % - rangeCombType_nomScen = 'combined'; % combine absolute and relative range scenarios - rangeGenType_nomScen = 'equidistant'; % equidistant: equidistant range shifts, - scenCombType_nomScen = 'individual'; % individual: no combination of range and setup scenarios, - includeNomScen_nomScen = true; - - % 'wcScen' default parameters for worst case scenarios - numOfShiftScen_wcScen = [2 2 2]; % number of shifts in x y and z direction - %shiftSize_wcScen = [4 4 4]; % given in [mm] - shiftGenType_wcScen = 'equidistant'; % equidistant: equidistant shifts - shiftCombType_wcScen = 'individual'; % individual: no combination of shift scenarios - numOfRangeShiftScen_wcScen = 2; % number of absolute and/or relative range scnearios. - maxAbsRangeShift_wcScen = 1; % maximum absolute over and undershoot in mm - maxRelRangeShift_wcScen = 3.5; % maximum relative over and undershoot in % - rangeCombType_wcScen = 'combined'; % combine absolute and relative range scenarios - rangeGenType_wcScen = 'equidistant'; % equidistant: equidistant range shifts - scenCombType_wcScen = 'individual'; % individual: no combination of range and setup scenarios - includeNomScen_wcScen = true; % include nominal scenario - - % 'impScen' create important/grid scenarios - numOfShiftScen_impScen = [10 10 10]; % number of shifts in x y and z direction - %shiftSize_impScen = [0 0 0]; % given in [mm] - shiftGenType_impScen = 'equidistant'; % equidistant: equidistant shifts - shiftCombType_impScen = 'individual'; % individual: no combination of shift scenarios - numOfRangeShiftScen_impScen = 10; % number of absolute and/or relative range scnearios - maxAbsRangeShift_impScen = 1; % maximum absolute over and undershoot in mm - maxRelRangeShift_impScen = 3.5; % maximum relative over and undershoot in % - rangeCombType_impScen = 'combined'; % combine absolute and relative range scenarios - rangeGenType_impScen = 'equidistant'; % equidistant: equidistant range shifts - scenCombType_impScen = 'individual'; % individual: no combination of range and setup scenarios, - includeNomScen_impScen = false; % exclude nominal scenario - - % 'rndScen' default parameters for random sampling - numOfShiftScen_rndScen = [25 25 25]; % number of shifts in x y and z direction - %shiftSize_rndScen = [3 3 3]; % given in [mm] - shiftGenType_rndScen = 'sampled'; % sample shifts from normal distribution - shiftCombType_rndScen = 'combined'; % individual: no combination of shift scenarios; - numOfRangeShiftScen_rndScen = 25; % number of absolute and/or relative range scnearios. - maxAbsRangeShift_rndScen = NaN; % maximum absolute over and undershoot in mm - maxRelRangeShift_rndScen = NaN; % maximum relative over and undershoot in % - rangeCombType_rndScen = 'combined'; % combine absolute and relative range scenarios - rangeGenType_rndScen = 'sampled'; % sampled: sample range shifts from normal distribution - scenCombType_rndScen = 'combined'; % combine range and setup scenarios if their scenario number is consistent - includeNomScen_rndScen = false; % exclude nominal scenario - end - - properties (Access = private) - lockUpdate = false; - lockInit = false; - end - - %% methods - - % public methods go here - methods - - %Attach a listener to a TYPE - function attachListener(obj) - addlistener(obj,'TYPE','PostSet',@PropLis.propChange); - end - - % default constructor - function this = matRad_multScen(ct,TYPE) - - matRad_cfg = MatRad_Config.instance(); - if isempty(ct) - this.numOfCtScen = 1; - else - this.numOfCtScen = ct.numOfCtScen; - end - - this.lockInit = true; - if exist('TYPE','var') && ~isempty(TYPE) - if sum(strcmp(this.AvailableScenCreationTYPE,TYPE))>0 - this.TYPE = TYPE; - else - matRad_cfg.dispWarning('matRad_multScen: Unknown TYPE - using the nominal scenario now'); - this.TYPE = this.DEFAULT_TYPE; - end - else - this.TYPE = this.DEFAULT_TYPE; - end - this.lockInit = false; - - this = this.initialize(); - - end % end constructor - - function newInstance = extractSingleNomScen(this,ctScen,i) - newInstance = this; - newInstance.TYPE = 'nomScen'; - newInstance.lockInit = true; - newInstance.lockUpdate = true; - - newInstance.scenForProb = this.scenForProb(i,:); - newInstance.relRangeShift = this.scenForProb(i,5); - newInstance.absRangeShift = this.scenForProb(i,4); - newInstance.isoShift = this.scenForProb(i,1:3); - newInstance.totNumShiftScen = 1; - newInstance.totNumRangeScen = 1; - newInstance.numOfCtScen = ctScen; - newInstance.scenMask = 1; - newInstance.linearMask = 1; - newInstance.scenProb = 1; - - newInstance.lockInit = false; - newInstance.lockUpdate = false; - end - - % create valid instance of an object - function this = matRad_createValidInstance(this) - this = setMultScen(this); - this = calcScenProb(this); - end - - %% Setter functions - function this = set.TYPE(this,value) - if ischar(value) && any(strcmp(value,this.AvailableScenCreationTYPE)) - this.TYPE = value; - this = this.initialize(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid scenario TYPE!'); - end - end - - function this = set.rangeRelSD(this,value) - if isnumeric(value) && isscalar(value) && value >= 0 - this.rangeRelSD = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for relative range uncertainty!'); - end - end - - function this = set.rangeAbsSD(this,value) - if isnumeric(value) && isscalar(value) && value >= 0 - this.rangeAbsSD = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for absolute range uncertainty!'); - end - end - - function this = set.shiftSD(this,value) - if isnumeric(value) && isvector(value) && numel(value) == 3 && all(value >= 0) - this.shiftSD = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for setup uncertainty!'); - end - end - - function this = set.wcFactor(this,value) - if isnumeric(value) && isscalar(value) && value >= 0 - this.wcFactor = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for worst case factor!'); - end - end - - function this = set.numOfShiftScen(this,value) - if isnumeric(value) && isvector(value) && numel(value) == 3 && all(value >= 0) && all(mod(value,1) == 0) - this.numOfShiftScen = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for shift scenario number!'); - end - end - - function this = set.numOfRangeShiftScen(this,value) - if isnumeric(value) && isscalar(value) && value >= 0 && mod(value,1) == 0 - this.numOfRangeShiftScen = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for number of range over/undershoots!'); - end - end - - function this = set.shiftGenType(this,value) - if ischar(value) && any(strcmp(value,{'equidistant','sampled'})) - this.shiftGenType = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for shift generation type!'); - end - end - - function this = set.rangeGenType(this,value) - if ischar(value) && any(strcmp(value,{'equidistant','sampled'})) - this.rangeGenType = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for range shift generation type!'); - end - end - - function this = set.shiftCombType(this,value) - if ischar(value) && any(strcmp(value,{'individual','combined','permuted'})) - this.shiftCombType = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for shift combination!'); - end - end - - function this = set.scenCombType(this,value) - if ischar(value) && any(strcmp(value,{'individual','combined','permuted'})) - this.scenCombType = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for scenario combination!'); - end - end - - function this = set.rangeCombType(this,value) - if ischar(value) && any(strcmp(value,{'individual','combined'})) - this.rangeCombType = value; - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for absolute & relative range combination!'); - end - end - - function this = set.includeNomScen(this,value) - if isnumeric(value) || islogical(value) - this.includeNomScen = logical(value); - this = this.updateScenariosFromUncertainties(); - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Invalid value for inclusion of nominal scenario!'); - end - end - - function listAllScenarios(this) - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispInfo('Listing all scenarios...\n'); - matRad_cfg.dispInfo('\t#\txShift\tyShift\tzShift\tabsRng\trelRng\tprob.\n'); - for s = 1:size(this.scenForProb,1) - str = num2str(this.scenForProb(s,:),'\t%.3f'); - matRad_cfg.dispInfo('\t%d\t%s\t%.3f\n',s,str,this.scenProb(s)); - end - end - end % end public methods - - %private methods go here - methods(Access = private) - %% - % set parameters according to the choosing scenario generation type - function this = initialize(this) - if this.lockInit - return; - end - - - this.lockUpdate = true; - - this.numOfShiftScen = this.(['numOfShiftScen_' this.TYPE]); - this.shiftGenType = this.(['shiftGenType_' this.TYPE]); - this.shiftCombType = this.(['shiftCombType_' this.TYPE]); - this.numOfRangeShiftScen = this.(['numOfRangeShiftScen_' this.TYPE]); - - this.rangeCombType = this.(['rangeCombType_' this.TYPE]); - this.rangeGenType = this.(['rangeGenType_' this.TYPE]); - - this.scenCombType = this.(['scenCombType_' this.TYPE]); - this.includeNomScen = this.(['includeNomScen_' this.TYPE]); - - this.lockUpdate = false; - - this = this.updateScenariosFromUncertainties(); - end - - %% - function this = updateScenariosFromUncertainties(this) - if this.lockUpdate - return; - end - switch this.TYPE - case 'wcScen' - this.shiftSize = this.shiftSD * this.wcFactor; - this.maxAbsRangeShift = this.rangeAbsSD * this.wcFactor; - this.maxRelRangeShift = this.rangeRelSD * this.wcFactor; - case 'nomScen' - this.shiftSize = [0 0 0]; - this.maxAbsRangeShift = 0; - this.maxRelRangeShift = 0; - case 'impScen' - this.shiftSize = this.shiftSD * this.wcFactor; - this.maxAbsRangeShift = this.rangeAbsSD * this.wcFactor; - this.maxRelRangeShift = this.rangeRelSD * this.wcFactor; - case 'rndScen' - this.shiftSize = [NaN NaN NaN]; - this.maxAbsRangeShift = this.rangeAbsSD * 3; - this.maxRelRangeShift = this.rangeRelSD * 3; - end - - this = this.setMultScen(); - this = this.calcScenProb(); - end - - %% - % creates individual treatment planning scenarios - function this = setMultScen(this) - matRad_cfg = MatRad_Config.instance(); - if this.includeNomScen - nomScen = 0; - else - nomScen = []; - end - - %% calc setup shift scenarios - switch this.shiftGenType - case 'equidistant' - % create grid vectors - isoShiftVec{1} = [0 linspace(-this.shiftSize(1), this.shiftSize(1), this.numOfShiftScen(1))]; - isoShiftVec{2} = [0 linspace(-this.shiftSize(2), this.shiftSize(2), this.numOfShiftScen(2))]; - isoShiftVec{3} = [0 linspace(-this.shiftSize(3), this.shiftSize(3), this.numOfShiftScen(3))]; - case 'sampled' - meanP = zeros(1,3); % mean (parameter) - if matRad_cfg.isMatlab - rng('shuffle'); - end - isoShiftVec{1} = [0 this.shiftSD(1) .* randn(1, this.numOfShiftScen(1)) + meanP(1)]; - isoShiftVec{2} = [0 this.shiftSD(2) .* randn(1, this.numOfShiftScen(2)) + meanP(2)]; - isoShiftVec{3} = [0 this.shiftSD(3) .* randn(1, this.numOfShiftScen(3)) + meanP(3)]; - otherwise - matRad_cfg.dispError('did not expect that!'); - end - - %Remove double shifts - isoShiftVec = cellfun(@unique,isoShiftVec,'UniformOutput',false); - - % create scenMaskIso for isoShifts - numIso(1) = numel(isoShiftVec{1}); - numIso(2) = numel(isoShiftVec{2}); - numIso(3) = numel(isoShiftVec{3}); - - scenMaskIso = false(numIso(1), numIso(2), numIso(3)); - - switch this.shiftCombType - case 'individual' - scenMaskIso(:,1,1) = true; % x shifts - scenMaskIso(1,:,1) = true; % y shifts - scenMaskIso(1,1,:) = true; % z shifts - case 'permuted' - scenMaskIso(:,:,:) = true; - case 'combined' - % determine that matrix is cubic - if isequal(numIso(1), numIso(2), numIso(3)) - for i = 1:numIso(1) - scenMaskIso(i,i,i) = true; - end - else - this.shiftCombType = 'individual'; - matRad_cfg.dispWarning('Numnber of isoShifts in every direction has to be equal in order to perform direct combination. Performing individually instead.'); - % call the function itself to get a working combination - [this] = setMultScen(this); - end - otherwise - matRad_cfg.dispError('Uncaught exception. Probably TYPO.','error'); - end - - % create list of increasing integers with referenced scenario - [xIso, yIso, zIso] = ind2sub(size(scenMaskIso),find(scenMaskIso)); - - matchMaskIso = cell(numel(xIso),2); - for i = 1:numel(xIso) - matchMaskIso{i,1} = i; - matchMaskIso{i,2} = [xIso(i) yIso(i) zIso(i)]; - end - - % create isoShift vector based on the matrix and matching - this.isoShift = zeros(size(matchMaskIso,1),3); - if numel(isoShiftVec{1}) + numel(isoShiftVec{2}) + numel(isoShiftVec{3}) > 0 - for i = 1:size(matchMaskIso,1) - matchPos = num2cell(matchMaskIso{i,2}); - if ~isequal([isoShiftVec{1}(matchPos{1}) isoShiftVec{2}(matchPos{2}) isoShiftVec{3}(matchPos{3})],[0 0 0]) || i == 1 - this.isoShift(i,:) = [isoShiftVec{1}(matchPos{1}) isoShiftVec{2}(matchPos{2}) isoShiftVec{3}(matchPos{3})] * ... - scenMaskIso(matchPos{:}); - end - end - end - - if ~this.includeNomScen - if size(this.isoShift,1)>1 - this.isoShift = this.isoShift(2:end,:); % cut away the first (=nominal) scenario - else - this.isoShift = []; - end - end - - this.totNumShiftScen = size(this.isoShift,1); % total number of shift scenarios in x,y and z direction - - %% calc range shift scenarios - switch this.rangeGenType - case 'equidistant' - this.relRangeShift = [nomScen linspace(-this.maxRelRangeShift, this.maxRelRangeShift, this.numOfRangeShiftScen)]; - this.absRangeShift = [nomScen linspace(-this.maxAbsRangeShift, this.maxAbsRangeShift, this.numOfRangeShiftScen)]; - case 'sampled' - % relRange - std = this.rangeRelSD; meanP = 0; - if matRad_cfg.isMatlab - rng('shuffle'); - end - this.relRangeShift = [nomScen std .* randn(1, this.numOfRangeShiftScen) + meanP]; - % absRange - std = this.rangeAbsSD; meanP = 0; - this.absRangeShift = [nomScen std .* randn(1, this.numOfRangeShiftScen) + meanP]; - otherwise - matRad_cfg.dispError('Not a valid type of generating data.'); - end - - numOfRelRangeShift = numel(this.relRangeShift); - numOfAbsRangeShift = numel(this.absRangeShift); - - if ~isequal(numOfRelRangeShift,numOfAbsRangeShift) - matRad_cfg.dispError('Number of relative and absolute range shifts must not differ.'); - else - this.totNumRangeScen = numOfRelRangeShift; - end - - this.relRangeShift = (this.relRangeShift./100); % consider [%] - - % check how absolute range error scenarios and relative range error scenarios should be combined - switch this.rangeCombType - case 'individual' - rangeShift = zeros(length(this.relRangeShift)*2 ,2); - for i = 1:length(this.relRangeShift) - rangeShift((i * 2)-1,1) = this.absRangeShift(i); - rangeShift((i * 2) ,2) = this.relRangeShift(i); - end - - this.totNumRangeScen = (numOfRelRangeShift*2)-1; - this.absRangeShift = [this.absRangeShift zeros(1,length(this.relRangeShift(this.relRangeShift~=0)))]; - this.relRangeShift = [zeros(1,length(this.relRangeShift)) this.relRangeShift(this.relRangeShift~=0)] ; - - case 'combined' - rangeShift = zeros(length(this.relRangeShift),2); - for i = 1:1:length(this.relRangeShift) - rangeShift(i,1) = this.absRangeShift(i); - rangeShift(i,2) = this.relRangeShift(i); - end - - otherwise - end - - % combine setup and range scenarios according to scenCombType - switch this.scenCombType - - case 'individual' % combine setup and range scenarios individually - - nIso = size(this.isoShift,1); - nRange = size(rangeShift,1); - - - % range errors should come last - if this.includeNomScen - nIso = nIso - 1; - nRange = nRange -1; - - this.scenForProb = zeros(1,5); %Nominal Scenario - - rangeScen = zeros(nRange,5); - rangeScen(:,4:5) = rangeShift(2:end,:); - - isoScen = zeros(nIso,5); - isoScen(:,1:3) = this.isoShift(2:end,:); - this.scenForProb = [this.scenForProb; isoScen; rangeScen]; - else - rangeScen = zeros(nRange,5); - isoScen = zeros(nIso,5); - rangeScen(:,4:5) = rangeShift; - isoScen(:,1:3) = this.isoShift; - - this.scenForProb = [isoScen; rangeScen]; - end - - case 'permuted' - - this.scenForProb = zeros(size(this.isoShift,1) * size(rangeShift,1),5); - Cnt = 1; - for j = 1:size(rangeShift,1) - for i = 1:size(this.isoShift,1) - this.scenForProb(Cnt,:) = [this.isoShift(i,:) rangeShift(j,:)]; - Cnt = Cnt + 1; - end - end - - case 'combined' - - if size(this.isoShift,1) == size(rangeShift,1) && this.totNumShiftScen > 0 && this.totNumRangeScen > 0 - this.scenForProb = zeros(size(this.isoShift,1),5); - this.scenForProb(1:end,1:3) = this.isoShift; - this.scenForProb(1:end,4:5) = rangeShift; - else - matRad_cfg.dispWarning('number of setup and range scenarios MUST be the same \n'); - this.scenCombType = 'individual'; - multScen = setMultScen(this); - this.scenForProb = multScen.scenForProb; - end - end - - % sanity check - UniqueRowScenForProb = unique(this.scenForProb,'rows'); - - if size(UniqueRowScenForProb,1) ~= size(this.scenForProb,1) && size(UniqueRowScenForProb,1)>1 - matRad_cfg.dispWarning('Some scenarios seem to be defined multiple times'); - end - - %% setup and fill combinatorics mask - % 1st dim: ct scenarios, - % 2nd dim: shift scenarios, - % 3rd dim: range scenarios - this.scenMask = false(this.numOfCtScen, this.totNumShiftScen, this.totNumRangeScen); - this.scenMask(:,1,1) = true; % ct scenarios - - % switch between combination modes here - % only makes scence when numOfShiftScen>0 and numOfRangeShiftScen>0; - if this.totNumShiftScen > 0 && this.totNumRangeScen > 0 - switch this.scenCombType - case 'individual' - % get all setup scenarios - % [~,ixUnq] = unique(this.scenForProb(:,1:3),'rows','stable'); - uq = this.scenForProb(1,1:3); - ixUnq = [1]; - for col = 2: size(this.scenForProb(:,1:3),1) - flag=1; - for colq = 1: size(uq,1) - if (sum(this.scenForProb(col,1:3) ==uq(colq,:)) == 3 ) - flag=0; - end - end - if flag - ixUnq = [ixUnq; col]; - end - end - this.scenMask = false(this.numOfCtScen, length(ixUnq), this.totNumRangeScen); - this.scenMask(:,1,1) = true; % ct scenarios - this.scenMask(1,:,1) = true; % iso shift scenarios - this.scenMask(1,1,:) = true; % range shift scenarios - case 'permuted' - this.scenMask(:,:,:) = true; - case 'combined' - % check if scenForProb is cubic (ignore ct scen) - if so fill diagonal - if isequal(this.totNumShiftScen, this.totNumRangeScen) - for i = 1:size(this.scenForProb, 1) - this.scenMask(1,i,i) = true; - end - else - this.shiftCombType = 'individual'; - matRad_cfg.dispWarning('number of setup and range scenarios MUST be the same \n'); - this = setMultScen(this); - end - end - elseif this.totNumShiftScen == 0 - this.scenMask(:,1,:) = true; - elseif this.totNumRangeScen == 0 - this.scenMask(:,:,1) = true; - 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)); - this.linearMask = cell2mat(x); - this.totNumScen = sum(this.scenMask(:)); - - end % end of setMultScen - - - %% - % provides different ways of calculating the probability - % of occurance of individual scenarios - function this = calcScenProb(this) - - mu = [0 0 0 0 0]; - sigma = [this.shiftSD this.rangeAbsSD this.rangeRelSD/100]; - - if isequal(this.DEFAULT_probDist,'normDist') - - this.scenProb = 1; - - if isequal(this.DEFAULT_TypeProp,'probBins') - - for i = 1:length(mu) - samplePosSorted = sort(unique(this.scenForProb(:,i))); - if numel(samplePosSorted) == 1 || sigma(i) == 0 - continue; - end - binWidth = (samplePosSorted(2) - samplePosSorted(1)); - lowerBinLevel = this.scenForProb(:,i) - 0.5*binWidth; - upperBinLevel = this.scenForProb(:,i) + 0.5*binWidth; - - this.scenProb = this.scenProb.*0.5.*(erf((upperBinLevel-mu(i))/(sqrt(2)*sigma(i)))-erf((lowerBinLevel-mu(i))/(sqrt(2)*sigma(i)))); - end - - elseif isequal(this.DEFAULT_TypeProp,'pointwise') - for i = 1:length(mu) - this.scenProb = this.scenProb .* (1/sqrt(2*pi*sigma(i)^2)*exp(-((this.scenForProb(:,i)-mu(i)).^2./(2*sigma(i)^2)))); - end - end - - % normalize probabilities since we use only a subset of - % the 3D grid - this.scenProb = this.scenProb./sum(this.scenProb); - - elseif isequal(probDist,'equalProb') - - numScen = size(samplePos,1); - this.scenProb = repmat(1/numScen,1,numScen); - - else - matRad_cfg = MatRad_Config.instance(); - matRad_cfg.dispError('Until now, only normally distributed scenarios implemented'); - end - - - end % end of calcScenProb - - - - - end % end private methods - - - -end % end of matRad_multScen class - - +function multScen = matRad_multScen(ct,scenarioModel) +% matRad_multScen +% This function replaces the deprecated matRad_multScen class. It creates +% the respective matRad_ScenarioModel instance for the specific scenario +% model chosen with standard parameters. +% +% call +% pln.multScen = matRad_multScen(ct,scenarioModel); +% +% e.g. pln.multScen = matRad_multScen(ct,'nomScen'); +% +% input +% ct: ct cube +% scenarioModel: string to denote a scenario creation method +% 'nomScen' create only the nominal scenario +% 'wcScen' create worst case scenarios +% 'impScen' create important/grid scenarios +% 'rndScen' create random scenarios +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +matRad_cfg = MatRad_Config.instance(); +matRad_cfg.dispWarning('The matRad_multScen function will be deprecated soon!\nCheck out the new Scenario Models in the scenarios folder.'); + +switch scenarioModel + case 'nomScen' + multScen = matRad_NominalScenario(ct); + case 'wcScen' + multScen = matRad_WorstCaseScenarios(ct); + case 'impScen' + multScen = matRad_ImportanceScenarios(ct); + case 'rndScen' + multScen = matRad_RandomScenarios(ct); + otherwise + matRad_cfg.dispError('''%s'' not known as scenario type!'); +end + +end diff --git a/matRad_sampling.m b/matRad_sampling.m index b0cf97331..3505d1a3b 100644 --- a/matRad_sampling.m +++ b/matRad_sampling.m @@ -51,9 +51,14 @@ pln.multScen = multScen; else % create random scenarios for sampling - pln.multScen = matRad_multScen(ctSamp,'rndScen'); % 'impSamp' or 'wcSamp' - pln.multScen.numOfShiftScen = matRad_cfg.defaults.samplingScenarios * ones(3,1); - pln.multScen.numOfRangeShiftScen = matRad_cfg.defaults.samplingScenarios; + % Deprecated Code + %pln.multScen = matRad_multScen(ctSamp,'rndScen'); % 'impSamp' or 'wcSamp' + %pln.multScen.numOfShiftScen = matRad_cfg.defaults.samplingScenarios * ones(3,1); + %pln.multScen.numOfRangeShiftScen = matRad_cfg.defaults.samplingScenarios; + %pln.multScen.totNumScen = matRad_cfg.defaults.samplingScenarios; %This yields an error now + + pln.multScen = matRad_RandomScenarios(ctSamp); + pln.multScen.nSamples = matRad_cfg.defaults.samplingScenarios; end matRad_cfg.dispInfo('Using %d samples in total \n',pln.multScen.totNumScen); @@ -141,7 +146,7 @@ FlagParforProgressDisp = false; end - for i = 1:pln.multScen.totNumScen + parfor i = 1:pln.multScen.totNumScen % create nominal scenario plnSamp = pln; diff --git a/matRad_setOverlapPriorities.m b/matRad_setOverlapPriorities.m index 6b6d9123d..5a389e3a7 100644 --- a/matRad_setOverlapPriorities.m +++ b/matRad_setOverlapPriorities.m @@ -61,7 +61,7 @@ cst{j,4}{i} = idx; if isempty(cst{j,4}{i}) && ~isempty(cst{j,6}) - matRad_cfg.dispWarning([cst{j,2} ': Objective(s) and/or constraints for inverse planning defined ' ... + matRad_cfg.dispWarning(['\n' cst{j,2} ': Objective(s) and/or constraints for inverse planning defined ' ... 'but structure overlapped by structure with higher overlap priority.' ... 'Objective(s) will not be considered during optimization']); end diff --git a/optimization/optimizer/compile_ipopt_minGW_octave520.sh b/optimization/optimizer/compile_ipopt_minGW_octave520.sh index aabcc12c7..8be00d39d 100644 --- a/optimization/optimizer/compile_ipopt_minGW_octave520.sh +++ b/optimization/optimizer/compile_ipopt_minGW_octave520.sh @@ -22,8 +22,13 @@ cd $IPOPTDIR/ThirdParty/Blas ./get.Blas cd $IPOPTDIR/ThirdParty/Lapack ./get.Lapack +# For Metis, we need to clone the updated git with the new get url +rm -rf $IPOPTDIR/ThirdParty/Metis +git clone https://github.com/coin-or-tools/ThirdParty-Metis.git --depth 1 --branch releases/1.3.10 $IPOPTDIR/ThirdParty/Metis cd $IPOPTDIR/ThirdParty/Metis ./get.Metis +rm -rf $IPOPTDIR/ThirdParty/Mumps +git clone https://github.com/coin-or-tools/ThirdParty-Mumps.git --depth 1 --branch releases/1.6.3 $IPOPTDIR/ThirdParty/Mumps cd $IPOPTDIR/ThirdParty/Mumps ./get.Mumps @@ -39,16 +44,16 @@ cp $MINGW_PREFIX/lib/gcc/$MINGW_CHOST/$GCC_VERSION/libquadmath.* /usr/lib/gcc/$M mkdir build cd build -../configure --prefix=$IPOPTINSTALLDIR --disable-shared --enable-static +../configure --prefix=$IPOPTINSTALLDIR --disable-shared --enable-static ADD_FFLAGS=-fallow-argument-mismatch make make install cd ../.. -# If everything worked, you should see some (static) libraries when doing ls /usr/local/lib +# If everything worked, you should see some (static) libraries when doing ls $IPOPTINSTALLDIR # we can get the mex interface -git clone https://github.com/ebertolazzi/mexIPOPT +git clone https://github.com/ebertolazzi/mexIPOPT --depth 1 --branch 1.0.0 # and compile it. mkoctfile --mex -ImexIPOPT/src -I$IPOPTINSTALLDIR/include/coin mexIPOPT/src/ipopt.cc mexIPOPT/src/IpoptInterfaceCommon.cc -v -DMATLAB_MEXFILE -DHAVE_CSTDDEF -DIPOPT_INTERFACE_MISSING_COPY_N -lipopt -lcoinmumps -lcoinmetis -lcoinlapack -lcoinblas -lgfortran -L$IPOPTINSTALLDIR/lib diff --git a/scenarios/matRad_ImportanceScenarios.m b/scenarios/matRad_ImportanceScenarios.m new file mode 100644 index 000000000..3d3b286de --- /dev/null +++ b/scenarios/matRad_ImportanceScenarios.m @@ -0,0 +1,278 @@ +classdef matRad_ImportanceScenarios < matRad_ScenarioModel +% matRad_ImportanceScenarios +% Implements gridded importance scenarios, i.e., weighted according to a +% probability distribution. It is not advised to create "combined" grids +% with a large number of grid-points as the curse of dimensionality will +% quickly break memory requirements when putting this in to dose influence +% matrix computation. +% +% +% constructor +% matRad_ImportanceScenarios() +% matRad_ImportanceScenarios(ct) +% +% input +% ct: ct cube +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + 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 + + properties (SetAccess=protected) + name = 'impScen'; + end + + properties (Constant) + validCombinationTypes = {'all','none','shift'}; + end + + methods + function this = matRad_ImportanceScenarios(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 + + 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)); + 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 + 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 new file mode 100644 index 000000000..99a08d969 --- /dev/null +++ b/scenarios/matRad_NominalScenario.m @@ -0,0 +1,95 @@ +classdef matRad_NominalScenario < matRad_ScenarioModel +% matRad_RandomScenarios +% Implements a single nominal planning scenario +% +% constructor +% matRad_NominalScenario() +% matRad_NominalScenario(ct) +% +% input +% ct: ct cube +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (SetAccess = protected) + name = 'nomScen'; + end + + methods + function this = matRad_NominalScenario(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 scenarios = updateScenarios(this) + %Scenario Probability from pdf - here it is one since only one + %scenario exist + %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; + + %Scenario weight + this.scenWeight = 1; + + %set variables + this.totNumShiftScen = 1; + this.totNumRangeScen = 1; + this.totNumScen = this.numOfCtScen; + + %Individual shifts + this.relRangeShift = 0; + this.absRangeShift = 0; + this.isoShift = [0 0 0]; + + this.maxAbsRangeShift = max(this.absRangeShift); + this.maxRelRangeShift = max(this.absRangeShift); + + %Mask for scenario selection + this.scenMask = true(this.numOfCtScen,this.totNumShiftScen,this.totNumRangeScen); + + %generic code + [x{1}, x{2}, x{3}] = ind2sub(size(this.scenMask),find(this.scenMask)); + this.linearMask = cell2mat(x); + totNumScen = sum(this.scenMask(:)); + + %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((this.scenForProb/cs).^2, 2)) / prod(diag(cs)); + this.scenWeight = this.scenProb./sum(this.scenProb); + + %Return variable + scenarios = [0 0 0 0 0]; + + if totNumScen ~= this.totNumScen + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispWarning('Check Implementation of Total Scenario computation - given %d but found %d!',this.totNumScen,totNumScen); + this.totNumScen = totNumScen; + end + end + + end +end + diff --git a/scenarios/matRad_RandomScenarios.m b/scenarios/matRad_RandomScenarios.m new file mode 100644 index 000000000..0a018df68 --- /dev/null +++ b/scenarios/matRad_RandomScenarios.m @@ -0,0 +1,162 @@ +classdef matRad_RandomScenarios < matRad_ScenarioModel +% matRad_RandomScenarios +% Implements randomly sampled scenarios +% +% constructor +% matRad_RandomScenarios() +% matRad_RandomScenarios(ct) +% +% input +% ct: ct cube +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties + includeNominalScenario = false; %Forces inclusion of the nominal scenario + nSamples = 10; %Standard number of random samples + end + + %Deprecated Properties that were used + properties (Dependent) + numOfShiftScen; + numOfRangeShiftScen; + end + + properties (SetAccess=protected) + name = 'rndScen'; + end + + methods + function this = matRad_RandomScenarios(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 + + %% Setters & Update + function set.nSamples(this,nSamples) + valid = isnumeric(nSamples) && isscalar(nSamples) && mod(nSamples,1) == 0 && nSamples > 0; + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for nSamples! Needs to be a positive integer!'); + end + this.nSamples = nSamples; + this.updateScenarios(); + end + + function set.numOfShiftScen(this,numOfShiftScen) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The property numOfShiftScen of the scenario class will soon be deprecated! Use nSamples instead'); + + %That's for downwards compatibility + if ~isscalar(numOfShiftScen) + numOfShiftScen = unique(numOfShiftScen); + end + + this.nSamples = numOfShiftScen; + end + + function value = get.numOfShiftScen(this) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The property numOfShiftScen of the scenario class will soon be deprecated! Use nSamples instead'); + value = this.nSamples; + end + + function set.numOfRangeShiftScen(this,numOfRangeShiftScen) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The property numOfRangeShiftScen of the scenario class will soon be deprecated! Use nSamples instead'); + this.nSamples = numOfRangeShiftScen; + end + + function value = get.numOfRangeShiftScen(this) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The property numOfRangeShiftScen of the scenario class will soon be deprecated! Use nSamples instead'); + value = this.nSamples; + end + + function scenarios = updateScenarios(this) + matRad_cfg = MatRad_Config.instance(); + + %Multivariate Normal Sampling + Sigma = diag([this.shiftSD,this.rangeAbsSD,this.rangeRelSD./100].^2); + d = size(Sigma,1); + [cs,p] = chol(Sigma); + + %The lower part is there but commented for completeness, we do + %not need it since we know we that sigma is PSD + %if p ~= 0 %for the positive semi-definite case, we don't check for negative definite + % % [V,L] = eig(Sigma); + % cs = sqrt(L) * V'; + %end + %transform normal samples, mean is always zero + scenarios = randn(this.nSamples,d) * cs; + + if this.includeNominalScenario + %We include the nominal scenario by just replacing the + %first one to keep the number of scenarios the same + scenarios(1,:) = 0; + end + + %Scenario Probability from pdf + this.scenForProb = scenarios; + this.scenProb = (2*pi)^(-d/2) * exp(-0.5*sum((scenarios/cs).^2, 2)) / prod(diag(cs)); + + %Scenario weight + this.scenWeight = ones(this.nSamples,1)./this.nSamples; %equal weights since they have been randomly sampled (not entirely true if the Nominal scenario was forced) + + %set variables + this.totNumShiftScen = this.nSamples; + this.totNumRangeScen = this.nSamples; + this.totNumScen = this.nSamples; %check because of CT scenarios + %this.totNumCtScen = + %this.numOfShiftScen = [nSamples,nSamples,nSamples]; + %this.numOfRangeShiftScen = nSamples; + + %Individual shifts + this.relRangeShift = scenarios(:,5); + this.absRangeShift = scenarios(:,4); + this.isoShift = scenarios(:,1:3); + + this.maxAbsRangeShift = max(this.absRangeShift); + this.maxRelRangeShift = max(this.relRangeShift); + + %Mask for scenario selection + this.scenMask = false(this.numOfCtScen,this.totNumShiftScen,this.totNumRangeScen); + + for sCt = 1:this.numOfCtScen + this.scenMask(sCt,:,:) = diag(true(this.nSamples,1)); + end + + [x{1}, x{2}, x{3}] = ind2sub(size(this.scenMask),find(this.scenMask)); + this.linearMask = cell2mat(x); + totNumScen = sum(this.scenMask(:)); + + if totNumScen ~= this.totNumScen + matRad_cfg.dispWarning('Check Implementation of Total Scenario computation - given %d but found %d!',this.totNumScen,totNumScen); + this.totNumScen = totNumScen; + end + end + end + + +end + diff --git a/scenarios/matRad_ScenarioModel.m b/scenarios/matRad_ScenarioModel.m new file mode 100644 index 000000000..c403836c8 --- /dev/null +++ b/scenarios/matRad_ScenarioModel.m @@ -0,0 +1,175 @@ +classdef matRad_ScenarioModel < handle +% matRad_ScenarioModel +% This is an abstract interface class to define Scenario Models for use in +% robust treatment planning and uncertainty analysis. +% Subclasses should at least implement the update() function to generate +% their own scenarios. +% +% constructor (Abstract) +% matRad_ScenarioModel() +% matRad_ScenarioModel(ct) +% +% input +% ct: ct cube +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + properties (AbortSet = true) %We use AbortSet = true here to avoid updates when + %Uncertainty model + 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 + end + + properties (Abstract,SetAccess=protected) + name + end + + properties (SetAccess = protected) + numOfCtScen; % total number of CT scenarios + + + % these parameters will be filled according to the choosen scenario type + isoShift; + relRangeShift; + absRangeShift; + + maxAbsRangeShift; + maxRelRangeShift; + + totNumShiftScen; % total number of shift scenarios in x,y and z direction + totNumRangeScen; % total number of range and absolute range scenarios + totNumScen; % total number of samples + + scenForProb; % matrix for probability calculation - each row denotes one scenario + scenProb; % probability of each scenario stored in a vector (according to uncertainty model) + scenWeight; % weight of scenario relative to the underlying uncertainty model (depends on how scenarios are chosen / sampled) + scenMask; + linearMask; + end + + methods + function this = matRad_ScenarioModel(ct) + if nargin == 0 || isempty(ct) + this.numOfCtScen = 1; + else + this.numOfCtScen = ct.numOfCtScen; + end + + %TODO: We could do this here automatically in the constructur, but + %Octave 5 has a bug here and throws an error + %this.updateScenarios(); + end + + %% SETTERS & UPDATE + function set.rangeRelSD(this,rangeRelSD) + valid = isnumeric(rangeRelSD) && isscalar(rangeRelSD) && rangeRelSD >= 0; + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for rangeRelSD! Needs to be a real positive scalar!'); + end + this.rangeRelSD = rangeRelSD; + this.updateScenarios(); + end + + function set.rangeAbsSD(this,rangeAbsSD) + valid = isnumeric(rangeAbsSD) && isscalar(rangeAbsSD) && rangeAbsSD >= 0; + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for rangeAbsSD! Needs to be a real positive scalar!'); + end + this.rangeAbsSD = rangeAbsSD; + this.updateScenarios(); + end + + function set.shiftSD(this,shiftSD) + valid = isnumeric(shiftSD) && isrow(shiftSD) && numel(shiftSD) == 3; + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for shiftSD! Needs to be 3-element numeric row vector!'); + end + this.shiftSD = shiftSD; + this.updateScenarios(); + end + + function set.wcSigma(this,wcSigma) + valid = isnumeric(wcSigma) && isscalar(wcSigma) && wcSigma >= 0; + if ~valid + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('Invalid value for wcSigma! Needs to be a real positive scalar!'); + end + this.wcSigma = wcSigma; + this.updateScenarios(); + end + + + + function scenarios = updateScenarios(this) + %This function will always update the scenarios given the + %current property settings + + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispError('This abstract function needs to be implemented!'); + end + + function newInstance = extractSingleScenario(this,scenIdx) + newInstance = matRad_NominalScenario(); + + newInstance.scenForProb = this.scenForProb(scenIdx,:); + newInstance.relRangeShift = this.scenForProb(scenIdx,5); + newInstance.absRangeShift = this.scenForProb(scenIdx,4); + newInstance.isoShift = this.scenForProb(scenIdx,1:3); + newInstance.scenProb = this.scenProb(scenIdx); + newInstance.scenWeight = this.scenWeight(scenIdx); + newInstance.numOfCtScen = this.numOfCtScen; + end + + %% Deprecated functions / properties + function newInstance = extractSingleNomScen(this,~,scenIdx) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The function extractSingleNomScen of the scenario class will soon be deprecated! Use extractSingleScenario instead!'); + newInstance = this.extractSingleScenario(scenIdx); + end + + function t = TYPE(this) + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The property TYPE of the scenario class will soon be deprecated!'); + t = this.name; + end + end + + methods (Static) + %{ + %TODO: implement automatic collection of available scenario classes + + function metaScenarioModels = getAvailableModels() + matRad_cfg = MatRad_Config.instance(); + + %Use the root folder and the scenarios folder only + folders = {matRad_cfg.matRadRoot,mfilename("fullpath")}; + + % + end + %} + + function types = AvailableScenCreationTYPE() + matRad_cfg = MatRad_Config.instance(); + matRad_cfg.dispDeprecationWarning('The function/property AvailableScenarioCreationTYPE of the scenario class will soon be deprecated!'); + %Hardcoded for compatability with matRad_multScen + types = {'nomScen','wcScen','impScen','rndScen'}; + end + end +end + diff --git a/scenarios/matRad_WorstCaseScenarios.m b/scenarios/matRad_WorstCaseScenarios.m new file mode 100644 index 000000000..8f2103161 --- /dev/null +++ b/scenarios/matRad_WorstCaseScenarios.m @@ -0,0 +1,125 @@ +classdef matRad_WorstCaseScenarios < matRad_ScenarioModel +% matRad_WorstCaseScenarios +% Implements single worst-case shifts per dimension.% +% +% constructor +% matRad_WorstCaseScenarios() +% matRad_WorstCaseScenarios(ct) +% +% input +% ct: ct cube +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSE.md. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the +% LICENSE file. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + 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'}; + end + + methods + function this = matRad_WorstCaseScenarios(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.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; + else + nPoints = 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); + end + end + + +end