diff --git a/4D/matRad_calc4dDose.m b/4D/matRad_calc4dDose.m index 1249d54e8..a658faa9c 100644 --- a/4D/matRad_calc4dDose.m +++ b/4D/matRad_calc4dDose.m @@ -55,6 +55,11 @@ timeSequence = []; end + +if any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'})) + [ax,bx] = matRad_getPhotonLQMParameters(cst,numel(resultGUI.physicalDose)); +end + % compute all phases for i = 1:ct.numOfCtScen @@ -70,6 +75,10 @@ elseif any(strcmp(pln.bioParam.model,{'MCN','LEM','WED','HEL'})) resultGUI.phaseAlphaDose{i} = tmpResultGUI.alpha .* tmpResultGUI.physicalDose; resultGUI.phaseSqrtBetaDose{i} = sqrt(tmpResultGUI.beta) .* tmpResultGUI.physicalDose; + ix = ax{i} ~=0; + resultGUI.phaseEffect{i} = resultGUI.phaseAlphaDose{i} + resultGUI.phaseSqrtBetaDose{i}.^2; + resultGUI.phaseRBExD{i} = zeros(ct.cubeDim); + resultGUI.phaseRBExD{i}(ix) = ((sqrt(ax{i}(ix).^2 + 4 .* bx{i}(ix) .* resultGUI.phaseEffect{i}(ix)) - ax{i}(ix))./(2.*bx{i}(ix))); end end @@ -89,12 +98,12 @@ resultGUI.accSqrtBetaDose = matRad_doseAcc(ct,resultGUI.phaseSqrtBetaDose, cst, accType); % only compute where we have biologically defined tissue - ix = any(cell2mat(cellfun(@(ax) ax ~= 0,dij.ax','UniformOutput',false)),2); + ix = (ax{1} ~= 0); resultGUI.accEffect = resultGUI.accAlphaDose + resultGUI.accSqrtBetaDose.^2; resultGUI.accRBExD = zeros(ct.cubeDim); - resultGUI.accRBExD(ix) = ((sqrt(dij.ax(ix).^2 + 4 .* dij.bx(ix) .* resultGUI.accEffect(ix)) - dij.ax(ix))./(2.*dij.bx(ix))); + resultGUI.accRBExD(ix) = ((sqrt(ax{1}(ix).^2 + 4 .* bx{1}(ix) .* resultGUI.accEffect(ix)) - ax{1}(ix))./(2.*bx{1}(ix))); end diff --git a/matRad_calcCubes.m b/matRad_calcCubes.m index 6a1511347..55c5c1ba1 100644 --- a/matRad_calcCubes.m +++ b/matRad_calcCubes.m @@ -48,6 +48,8 @@ beamInfo(dij.numOfBeams+1).suffix = ''; beamInfo(dij.numOfBeams+1).logIx = true(size(resultGUI.w,1),1); +[ctScen,~] = ind2sub(size(dij.physicalDose),scenNum); + %% Physical Dose doseFields = {'physicalDose','doseToWater'}; @@ -112,7 +114,7 @@ wBeam = (resultGUI.w .* beamInfo(i).logIx); % consider biological optimization - ix = dij.bx{scenNum}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0; + ix = dij.bx{ctScen}~=0 & resultGUI.(['physicalDose', beamInfo(i).suffix])(:) > 0; % Calculate effect from alpha- and sqrtBetaDose resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]) = full(dij.(['mAlphaDose' RBE_model{j}]){scenNum} * wBeam + (dij.(['mSqrtBetaDose' RBE_model{j}]){scenNum} * wBeam).^2); @@ -120,7 +122,7 @@ % Calculate RBExD from the effect resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix]) = zeros(size(resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix]))); - resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix) = (sqrt(dij.ax{scenNum}(ix).^2 + 4 .* dij.bx{scenNum}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{scenNum}(ix))./(2.*dij.bx{scenNum}(ix)); + resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])(ix) = (sqrt(dij.ax{ctScen}(ix).^2 + 4 .* dij.bx{ctScen}(ix) .* resultGUI.(['effect', RBE_model{j}, beamInfo(i).suffix])(ix)) - dij.ax{ctScen}(ix))./(2.*dij.bx{ctScen}(ix)); % Divide RBExD with the physicalDose to get the plain RBE cube resultGUI.(['RBE', RBE_model{j}, beamInfo(i).suffix]) = resultGUI.(['RBExD', RBE_model{j}, beamInfo(i).suffix])./resultGUI.(['physicalDose', beamInfo(i).suffix]); diff --git a/matRad_fluenceOptimization.asv b/matRad_fluenceOptimization.asv new file mode 100644 index 000000000..84f8fc1a3 --- /dev/null +++ b/matRad_fluenceOptimization.asv @@ -0,0 +1,343 @@ +function [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln,wInit) +% matRad inverse planning wrapper function +% +% call +% [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln) +% [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln,wInit) +% +% input +% dij: matRad dij struct +% cst: matRad cst struct +% pln: matRad pln struct +% wInit: (optional) custom weights to initialize problems +% +% output +% resultGUI: struct containing optimized fluence vector, dose, and (for +% biological optimization) RBE-weighted dose etc. +% optimizer: Used Optimizer Object +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2016 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(); + +% consider VOI priorities +cst = matRad_setOverlapPriorities(cst); + +% check & adjust objectives and constraints internally for fractionation +for i = 1:size(cst,1) + %Compatibility Layer for old objective format + if isstruct(cst{i,6}) + cst{i,6} = arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct,cst{i,6},'UniformOutput',false); + end + for j = 1:numel(cst{i,6}) + + obj = cst{i,6}{j}; + + %In case it is a default saved struct, convert to object + %Also intrinsically checks that we have a valid optimization + %objective or constraint function in the end + if ~isa(obj,'matRad_DoseOptimizationFunction') + try + obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj); + catch + matRad_cfg.dispError('cst{%d,6}{%d} is not a valid Objective/constraint! Remove or Replace and try again!',i,j); + end + end + + obj = obj.setDoseParameters(obj.getDoseParameters()/pln.numOfFractions); + + cst{i,6}{j} = obj; + end +end + +% resizing cst to dose cube resolution +cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x, dij.ctGrid.y, dij.ctGrid.z,... + dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z); + +% Get rid of voxels that are not interesting for the optimization problem +if ~isfield(pln.propOpt, 'clearUnusedVoxels') + pln.propOpt.clearUnusedVoxels = matRad_cfg.propOpt.defaultClearUnusedVoxels; +end + +if pln.propOpt.clearUnusedVoxels + dij = matRad_clearUnusedVoxelsFromDij(cst, dij); +end + + +% find target indices and described dose(s) for weight vector +% initialization +V = []; +doseTarget = []; +ixTarget = []; + +for i = 1:size(cst,1) + if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6}) + V = [V;cst{i,4}{1}]; + + %Iterate through objectives/constraints + fDoses = []; + for fObjCell = cst{i,6} + dParams = fObjCell{1}.getDoseParameters(); + %Don't care for Inf constraints + dParams = dParams(isfinite(dParams)); + %Add do dose list + fDoses = [fDoses dParams]; + end + + doseTarget = [doseTarget fDoses]; + ixTarget = [ixTarget i*ones(1,length(fDoses))]; + end +end +[doseTarget,i] = max(doseTarget); +ixTarget = ixTarget(i); +wOnes = ones(dij.totalNumOfBixels,1); + +% calculate initial beam intensities wInit +matRad_cfg.dispInfo('Estimating initial weights... '); + +if exist('wInit','var') + %do nothing as wInit was passed to the function + matRad_cfg.dispInfo('chosen provided wInit!\n'); + + % Write ixDose which is needed for the optimizer + if pln.bioParam.bioOpt + dij.ixDose = dij.bx~=0; + + %pre-calculations + dij.gamma = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); + dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose)); + end + +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; + end + + doseTmp = dij.physicalDose{1}*wOnes; + bixelWeight = (doseTarget)/(dij.RBE * mean(doseTmp(V))); + wInit = wOnes * bixelWeight; + matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); + +elseif pln.bioParam.bioOpt + % retrieve photon LQM parameter + [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels); + checkAxBx = cellfun(@(ax1,bx1,ax2,bx2) isequal(ax1(ax1~=0),ax2(ax1~=0)) && isequal(bx1(bx1~=0),bx2(bx1~=0)),dij.ax,dij.bx,ax,bx); + if ~all(checkAxBx) + matRad_cfg.dispError('Inconsistent biological parameters in dij.ax and/or dij.bx - please recalculate dose influence matrix before optimization!\n'); + end + + + + for i = 1:size(cst,1) + + 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 > 5 Gy[RBE] for target. Biological optimization outside the valid domain of the base data. Reduce dose prescription or use more fractions.\n'); + end + + end + end + + for s = 1:numel(dij.bx) + dij.ixDose{s} = dij.bx{s}~=0; + end + + if isequal(pln.bioParam.quantityOpt,'effect') + + effectTarget = cst{ixTarget,5}.alphaX * doseTarget + cst{ixTarget,5}.betaX * doseTarget^2; + aTmp = dij.mAlphaDose{1}*wOnes; + bTmp = dij.mSqrtBetaDose{1} * wOnes; + p = sum(aTmp(V)) / sum(bTmp(V).^2); + q = -(effectTarget * length(V)) / sum(bTmp(V).^2); + + wInit = -(p/2) + sqrt((p^2)/4 -q) * wOnes; + + elseif isequal(pln.bioParam.quantityOpt,'RBExD') + + %pre-calculations + for s = 1:numel(dij.ixDose) + dij.gamma{s} = zeros(dij.doseGrid.numOfVoxels,dij.numOfScenarios); + dij.gamma{s}(dij.ixDose{s}) = dij.ax{s}(dij.ixDose{s})./(2*dij.bx{s}(dij.ixDose{s})); + end + + + % calculate current effect in target + aTmp = dij.mAlphaDose{1}*wOnes; + bTmp = dij.mSqrtBetaDose{1} * wOnes; + doseTmp = dij.physicalDose{1}*wOnes; + + CurrEffectTarget = aTmp(V) + bTmp(V).^2; + % ensure a underestimated biological effective dose + TolEstBio = 1.2; + % calculate maximal RBE in target + maxCurrRBE = max(-cst{ixTarget,5}.alphaX + sqrt(cst{ixTarget,5}.alphaX^2 + ... + 4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*doseTmp(V))); + wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(doseTmp(V))))* wOnes; + end + + 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; + matRad_cfg.dispInfo('chosen uniform weight of %f!\n',bixelWeight); +end + + +%% calculate probabilistic quantities for probabilistic optimization if at least +% one robust objective is defined + +%Check how to use 4D data +if isfield(pln,'propOpt') && isfield(pln.propOpt,'scen4D') + scen4D = pln.propOpt.scen4D; +else + scen4D = 1; %Use only first 4D scenario for optimization +end + +%If "all" provided, use all scenarios +if isequal(scen4D,'all') + scen4D = 1:size(dij.physicalDose,1); +end + +linIxDIJ = find(~cellfun(@isempty,dij.physicalDose(scen4D,:,:)))'; + +%Only select the indexes of the nominal ct Scenarios +linIxDIJ_nominalCT = find(~cellfun(@isempty,dij.physicalDose(scen4D,1,1)))'; + +FLAG_CALC_PROB = false; +FLAG_ROB_OPT = false; + + +for i = 1:size(cst,1) + for j = 1:numel(cst{i,6}) + if strcmp(cst{i,6}{j}.robustness,'PROB') && numel(linIxDIJ) > 1 + FLAG_CALC_PROB = true; + end + if ~strcmp(cst{i,6}{j}.robustness,'none') && numel(linIxDIJ) > 1 + FLAG_ROB_OPT = true; + end + end +end + +if FLAG_CALC_PROB + [dij] = matRad_calculateProbabilisticQuantities(dij,cst,pln); +end + + +% set optimization options +if ~FLAG_ROB_OPT || FLAG_CALC_PROB % if multiple robust objectives are defined for one structure then remove FLAG_CALC_PROB from the if clause + ixForOpt = scen4D; +else + ixForOpt = linIxDIJ; +end + +switch pln.bioParam.quantityOpt + case 'effect' + backProjection = matRad_EffectProjection; + case 'RBExD' + %Capture special case of constant RBE + if strcmp(pln.bioParam.model,'constRBE') + backProjection = matRad_ConstantRBEProjection; + else + backProjection = matRad_VariableRBEProjection; + end + case 'physicalDose' + backProjection = matRad_DoseProjection; + otherwise + warning(['Did not recognize bioloigcal setting ''' pln.probOpt.bioOptimization '''!\nUsing physical dose optimization!']); + backProjection = matRad_DoseProjection; +end + +%Give scenarios used for optimization +backProjection.scenarios = ixForOpt; +backProjection.scenarioProb = pln.multScen.scenProb; +backProjection.nominalCtScenarios = linIxDIJ_nominalCT; +%backProjection.scenDim = pln.multScen + +optiProb = matRad_OptimizationProblem(backProjection); +optiProb.quantityOpt = pln.bioParam.quantityOpt; +if isfield(pln,'propOpt') && isfield(pln.propOpt,'useLogSumExpForRobOpt') + optiProb.useLogSumExpForRobOpt = pln.propOpt.useLogSumExpForRobOpt; +end + +%Get Bounds +if ~isfield(pln.propOpt,'boundMU') + pln.propOpt.boundMU = false; +end + +if pln.propOpt.boundMU + if (isfield(dij,'minMU') || isfield(dij,'maxMU')) && ~isfield(dij,'numParticlesPerMU') + matRad_cfg.dispWarning('Requested MU bounds but number of particles per MU not set! Bounds will not be enforced and standard [0,Inf] will be used instead!'); + elseif ~isfield(dij,'minMU') && ~isfield(dij,'maxMU') + matRad_cfg.dispWarning('Requested MU bounds but machine bounds not defined in dij.minMU & dij.maxMU! Bounds will not be enforced and standard [0,Inf] will be used instead!'); + else + if isfield(dij,'minMU') + optiProb.minimumW = dij.numParticlesPerMU .* dij.minMU / 1e6; + matRad_cfg.dispInfo('Using lower MU bounds provided in dij!\n') + end + + if isfield(dij,'maxMU') + optiProb.maximumW = dij.numParticlesPerMU .* dij.maxMU / 1e6; + matRad_cfg.dispInfo('Using upper MU bounds provided in dij!\n') + end + end +else + matRad_cfg.dispInfo('Using standard MU bounds of [0,Inf]!\n') +end + +if ~isfield(pln.propOpt,'optimizer') + pln.propOpt.optimizer = 'IPOPT'; +end + +switch pln.propOpt.optimizer + case 'IPOPT' + optimizer = matRad_OptimizerIPOPT; + case 'fmincon' + optimizer = matRad_OptimizerFmincon; + otherwise + warning(['Optimizer ''' pln.propOpt.optimizer ''' not known! Fallback to IPOPT!']); + optimizer = matRad_OptimizerIPOPT; +end + +if ~optimizer.IsAvailable() + matRad_cfg.dispError(['Optimizer ''' pln.propOpt.optimizer ''' not available!']); +end + +optimizer = optimizer.optimize(wInit,optiProb,dij,cst); + +wOpt = optimizer.wResult; +info = optimizer.resultInfo; + +resultGUI = matRad_calcCubes(wOpt,dij); +resultGUI.wUnsequenced = wOpt; +resultGUI.usedOptimizer = optimizer; +resultGUI.info = info; + +%Robust quantities +if FLAG_ROB_OPT || numel(ixForOpt) > 1 + Cnt = 1; + for i = find(~cellfun(@isempty,dij.physicalDose))' + tmpResultGUI = matRad_calcCubes(wOpt,dij,i); + resultGUI.([pln.bioParam.quantityVis '_' num2str(Cnt,'%d')]) = tmpResultGUI.(pln.bioParam.quantityVis); + Cnt = Cnt + 1; + end +end + +% unblock mex files +clear mex diff --git a/optimization/projections/matRad_VariableRBEProjection.m b/optimization/projections/matRad_VariableRBEProjection.m index 36dc53177..0e13495fd 100644 --- a/optimization/projections/matRad_VariableRBEProjection.m +++ b/optimization/projections/matRad_VariableRBEProjection.m @@ -21,7 +21,7 @@ function RBExD = computeSingleScenario(obj,dij,scen,w) effect = computeSingleScenario@matRad_EffectProjection(obj,dij,scen,w); %First compute effect RBExD = zeros(dij.doseGrid.numOfVoxels,1); - ctScen = ind2sub(size(dij.physicalDose),scen); + [ctScen,~] = ind2sub(size(dij.physicalDose),scen); RBExD(dij.ixDose{ctScen}) = sqrt((effect(dij.ixDose{ctScen})./dij.bx{ctScen}(dij.ixDose{ctScen}))+(dij.gamma{ctScen}(dij.ixDose{ctScen}).^2)) - dij.gamma{ctScen}(dij.ixDose{ctScen}); end @@ -36,7 +36,7 @@ obj = obj.compute(dij,w); %Get corresponding ct scenario - ctScen = ind2sub(size(dij.physicalDose),scen); + [ctScen,~] = ind2sub(size(dij.physicalDose),scen); %Scaling vor variable RBExD scaledEffect = obj.d{scen} + dij.gamma{ctScen}; @@ -65,7 +65,7 @@ eExp = eExpLinTerm + eExpSqTerm.^2; %Get corresponding ct scenario - ctScen = ind2sub(size(dij.physicalDose),scen); + [ctScen,~] = ind2sub(size(dij.physicalDose),scen); RBExDexp = zeros(dij.doseGrid.numOfVoxels,1); RBExDexp(dij.ixDose{ctScen}) = sqrt((eExp(dij.ixDose{ctScen})./dij.bx{ctScen}(dij.ixDose{ctScen}))+(dij.gamma{ctScen}(dij.ixDose{ctScen}).^2)) - dij.gamma{ctScen}(dij.ixDose{ctScen}); @@ -85,7 +85,7 @@ obj = obj.compute(dij,w); %Get corresponding ct scenario - ctScen = ind2sub(size(dij.physicalDose),scen); + [ctScen,~]= ind2sub(size(dij.physicalDose),scen); %Scaling vor variable RBExD scaledEffect = obj.dExp{scen} + dij.gamma{ctScen};