matRad_fluenceOptimization

 matRad inverse planning wrapper function
 
 call
   [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln)

 input
   dij:        matRad dij struct
   cst:        matRad cst struct
   pln:        matRad pln struct

 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.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln)
0002 % matRad inverse planning wrapper function
0003 %
0004 % call
0005 %   [resultGUI,optimizer] = matRad_fluenceOptimization(dij,cst,pln)
0006 %
0007 % input
0008 %   dij:        matRad dij struct
0009 %   cst:        matRad cst struct
0010 %   pln:        matRad pln struct
0011 %
0012 % output
0013 %   resultGUI:  struct containing optimized fluence vector, dose, and (for
0014 %               biological optimization) RBE-weighted dose etc.
0015 %   optimizer:  Used Optimizer Object
0016 %
0017 % References
0018 %   -
0019 %
0020 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 %
0022 % Copyright 2016 the matRad development team.
0023 %
0024 % This file is part of the matRad project. It is subject to the license
0025 % terms in the LICENSE file found in the top-level directory of this
0026 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0027 % of the matRad project, including this file, may be copied, modified,
0028 % propagated, or distributed except according to the terms contained in the
0029 % LICENSE file.
0030 %
0031 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 
0033 matRad_cfg = MatRad_Config.instance();
0034 
0035 % issue warning if biological optimization impossible
0036 if sum(strcmp(pln.propOpt.bioOptimization,{'LEMIV_effect','LEMIV_RBExD'}))>0 && (~isfield(dij,'mAlphaDose') || ~isfield(dij,'mSqrtBetaDose')) && strcmp(pln.radiationMode,'carbon')
0037     warndlg('Alpha and beta matrices for effect based and RBE optimization not available - physical optimization is carried out instead.');
0038     pln.propOpt.bioOptimization = 'none';
0039 end
0040 
0041 % consider VOI priorities
0042 cst  = matRad_setOverlapPriorities(cst);
0043 
0044 
0045 % check & adjust objectives and constraints internally for fractionation
0046 for i = 1:size(cst,1)
0047     %Compatibility Layer for old objective format
0048     if isstruct(cst{i,6})
0049         cst{i,6} = arrayfun(@matRad_DoseOptimizationFunction.convertOldOptimizationStruct,cst{i,6},'UniformOutput',false);
0050     end
0051     for j = 1:numel(cst{i,6})
0052         
0053         obj = cst{i,6}{j};        
0054         
0055         %In case it is a default saved struct, convert to object
0056         %Also intrinsically checks that we have a valid optimization
0057         %objective or constraint function in the end
0058         if ~isa(obj,'matRad_DoseOptimizationFunction')
0059             try
0060                 obj = matRad_DoseOptimizationFunction.createInstanceFromStruct(obj);
0061             catch
0062                 matRad_cfg.dispError('cst{%d,6}{%d} is not a valid Objective/constraint! Remove or Replace and try again!',i,j);
0063             end
0064         end
0065         
0066         obj = obj.setDoseParameters(obj.getDoseParameters()/pln.numOfFractions);
0067         
0068         cst{i,6}{j} = obj;        
0069     end
0070 end
0071 
0072 % resizing cst to dose cube resolution
0073 cst = matRad_resizeCstToGrid(cst,dij.ctGrid.x,dij.ctGrid.y,dij.ctGrid.z,...
0074                                  dij.doseGrid.x,dij.doseGrid.y,dij.doseGrid.z);
0075 
0076 % find target indices and described dose(s) for weight vector
0077 % initialization
0078 V          = [];
0079 doseTarget = [];
0080 ixTarget   = [];
0081 
0082 for i = 1:size(cst,1)
0083     if isequal(cst{i,3},'TARGET') && ~isempty(cst{i,6})
0084         V = [V;cst{i,4}{1}];
0085         
0086         %Iterate through objectives/constraints
0087         fDoses = [];
0088         for fObjCell = cst{i,6}
0089             dParams = fObjCell{1}.getDoseParameters();
0090             %Don't care for Inf constraints
0091             dParams = dParams(isfinite(dParams));
0092             %Add do dose list
0093             fDoses = [fDoses dParams];
0094         end
0095                 
0096         
0097         doseTarget = [doseTarget fDoses];
0098         ixTarget   = [ixTarget i*ones(1,length(fDoses))];
0099     end
0100 end
0101 [doseTarget,i] = max(doseTarget);
0102 ixTarget       = ixTarget(i);
0103 wOnes          = ones(dij.totalNumOfBixels,1);
0104 
0105 % modified settings for photon dao
0106 if pln.propOpt.runDAO && strcmp(pln.radiationMode,'photons')
0107 %    options.ipopt.max_iter = 50;
0108 %    options.ipopt.acceptable_obj_change_tol     = 7e-3; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled
0109 
0110 end
0111 % calculate initial beam intensities wInit
0112 if  strcmp(pln.propOpt.bioOptimization,'const_RBExD') && strcmp(pln.radiationMode,'protons')
0113     
0114     % check if a constant RBE is defined - if not use 1.1
0115     if ~isfield(dij,'RBE')
0116         dij.RBE = 1.1;
0117     end
0118     bixelWeight =  (doseTarget)/(dij.RBE * mean(dij.physicalDose{1}(V,:)*wOnes)); 
0119     wInit       = wOnes * bixelWeight;
0120         
0121 elseif (strcmp(pln.propOpt.bioOptimization,'LEMIV_effect') || strcmp(pln.propOpt.bioOptimization,'LEMIV_RBExD')) ... 
0122                                 && strcmp(pln.radiationMode,'carbon')
0123                             
0124     % retrieve photon LQM parameter
0125     [ax,bx] = matRad_getPhotonLQMParameters(cst,dij.doseGrid.numOfVoxels,1);
0126 
0127     if ~isequal(dij.ax(dij.ax~=0),ax(dij.ax~=0)) || ...
0128        ~isequal(dij.bx(dij.bx~=0),bx(dij.bx~=0))
0129          matRad_cfg.dispError('Inconsistent biological parameter - please recalculate dose influence matrix!\n');
0130     end
0131 
0132     for i = 1:size(cst,1)
0133         
0134         for j = 1:size(cst{i,6},2)
0135             % check if prescribed doses are in a valid domain
0136             if any(cst{i,6}{j}.getDoseParameters() > 5) && isequal(cst{i,3},'TARGET')
0137                 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');
0138             end
0139             
0140         end
0141     end
0142     
0143     dij.ixDose  = dij.bx~=0; 
0144         
0145     if isequal(pln.propOpt.bioOptimization,'LEMIV_effect')
0146         
0147            effectTarget = cst{ixTarget,5}.alphaX * doseTarget + cst{ixTarget,5}.betaX * doseTarget^2;
0148            p            = (sum(dij.mAlphaDose{1}(V,:)*wOnes)) / (sum((dij.mSqrtBetaDose{1}(V,:) * wOnes).^2));
0149            q            = -(effectTarget * length(V)) / (sum((dij.mSqrtBetaDose{1}(V,:) * wOnes).^2));
0150            wInit        = -(p/2) + sqrt((p^2)/4 -q) * wOnes;
0151 
0152     elseif isequal(pln.propOpt.bioOptimization,'LEMIV_RBExD')
0153         
0154            %pre-calculations
0155            dij.gamma              = zeros(dij.doseGrid.numOfVoxels,1);   
0156            dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose)); 
0157             
0158            % calculate current in target
0159            CurrEffectTarget = (dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^2);
0160            % ensure a underestimated biological effective dose
0161            TolEstBio        = 1.2;
0162            % calculate maximal RBE in target
0163            maxCurrRBE = max(-cst{ixTarget,5}.alphaX + sqrt(cst{ixTarget,5}.alphaX^2 + ...
0164                         4*cst{ixTarget,5}.betaX.*CurrEffectTarget)./(2*cst{ixTarget,5}.betaX*(dij.physicalDose{1}(V,:)*wOnes)));
0165            wInit    =  ((doseTarget)/(TolEstBio*maxCurrRBE*max(dij.physicalDose{1}(V,:)*wOnes)))* wOnes;
0166     end
0167     
0168 else 
0169     bixelWeight =  (doseTarget)/(mean(dij.physicalDose{1}(V,:)*wOnes)); 
0170     wInit       = wOnes * bixelWeight;
0171     pln.propOpt.bioOptimization = 'none';
0172 end
0173 
0174 % set optimization options
0175 options.radMod          = pln.radiationMode;
0176 options.bioOpt          = pln.propOpt.bioOptimization;
0177 options.ID              = [pln.radiationMode '_' pln.propOpt.bioOptimization];
0178 options.numOfScenarios  = dij.numOfScenarios;
0179 
0180 %Select Projection
0181 
0182 switch pln.propOpt.bioOptimization
0183     case 'LEMIV_effect'
0184         backProjection = matRad_EffectProjection;
0185     case 'const_RBExD'
0186         backProjection = matRad_ConstantRBEProjection;
0187     case 'LEMIV_RBExD'
0188         backProjection = matRad_VariableRBEProjection;
0189     case 'none'
0190         backProjection = matRad_DoseProjection;
0191     otherwise
0192         warning(['Did not recognize bioloigcal setting ''' pln.probOpt.bioOptimization '''!\nUsing physical dose optimization!']);
0193         backProjection = matRad_DoseProjection;
0194 end
0195         
0196 
0197 %backProjection = matRad_DoseProjection();
0198 
0199 optiProb = matRad_OptimizationProblem(backProjection);
0200 
0201 %optimizer = matRad_OptimizerIPOPT;
0202 
0203 if ~isfield(pln.propOpt,'optimizer')
0204     pln.propOpt.optimizer = 'IPOPT';
0205 end
0206 
0207 switch pln.propOpt.optimizer
0208     case 'IPOPT'
0209         optimizer = matRad_OptimizerIPOPT;
0210     case 'fmincon'
0211         optimizer = matRad_OptimizerFmincon;
0212     otherwise
0213         warning(['Optimizer ''' pln.propOpt.optimizer ''' not known! Fallback to IPOPT!']);
0214         optimizer = matRad_OptimizerIPOPT;
0215 end
0216         
0217 %optimizer = matRad_OptimizerFmincon;
0218 
0219 optimizer = optimizer.optimize(wInit,optiProb,dij,cst);
0220 
0221 wOpt = optimizer.wResult;
0222 info = optimizer.resultInfo;
0223 
0224 resultGUI = matRad_calcCubes(wOpt,dij);
0225 resultGUI.wUnsequenced = wOpt;
0226 resultGUI.usedOptimizer = optimizer;
0227 resultGUI.info = info;
0228 
0229 % unblock mex files
0230 clear mex