Fixed example and implemented preconditioner within fluence optimization

matRad_MixedModality.m

@@ -2,21 +2,21 @@
 matRad_cfg = MatRad_Config.instance();
 matRad_cfg.propOpt.defaultMaxIter = 50000;
 load 'TG119.mat'
-
-% 
+cst{3,6}{2} = struct(DoseObjectives.matRad_MeanDose(100,0,1));
+%% 
 % meta information for treatment plan (1) 
 pln(1).numOfFractions  = 5;
 pln(1).radiationMode   = 'protons';           % either photons / protons / helium / carbon
 pln(1).machine         = 'Generic';
 
 % beam geometry settings
 pln(1).propStf.bixelWidth      = 5; % [mm] / also corresponds to lateral spot spacing for particles
-pln(1).propStf.gantryAngles    = [0]; % [?] ;
+pln(1).propStf.gantryAngles    = [ -45 0 45 ]; % [?] ;
 pln(1).propStf.couchAngles     = zeros(numel(pln(1).propStf.gantryAngles),1); % [?] ; 
 pln(1).propStf.numOfBeams      = numel(pln(1).propStf.gantryAngles);
 pln(1).propStf.isoCenter       = ones(pln(1).propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
 % optimization settings
-pln(1).propDoseCalc.calcLET = 1;
+pln(1).propDoseCalc.calcLET = 0;
 
 pln(1).propOpt.runDAO          = false;      % 1/true: run DAO, 0/false: don't / will be ignored for particles
 pln(1).propOpt.runSequencing   = false;      % 1/true: run sequencing, 0/false: don't / will be ignored for particles and also triggered by runDAO below
@@ -51,7 +51,7 @@
 
 % beam geometry settings
 pln(2).propStf.bixelWidth      = 5; % [mm] / also corresponds to lateral spot spacing for particles
-pln(2).propStf.gantryAngles    = [0:90:359]; % [?] ;
+pln(2).propStf.gantryAngles    = [0:72:359]; % [?] ;
 pln(2).propStf.couchAngles     = zeros(numel(pln(2).propStf.gantryAngles),1);  % [?] ; 
 pln(2).propStf.numOfBeams      = numel(pln(2).propStf.gantryAngles);
 pln(2).propStf.isoCenter       = ones(pln(2).propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
@@ -94,6 +94,8 @@
 stf = matRad_stfWrapper(ct,cst,plnJO);
 % Dij Calculation
 dij = matRad_calcCombiDose(ct,stf,plnJO,cst,false);
+dij.precon = 1;
+% dij = matRad_mixModPreconditioner(dij);
 % Fluence optimization 
 resultGUI = matRad_fluenceOptimizationJO(dij,cst,plnJO);
 

matRad_fluenceOptimizationJO.m

@@ -339,6 +339,9 @@
     optiProb.useLogSumExpForRobOpt = pln.propOpt.useLogSumExpForRobOpt;
 end
 
+if dij.precon
+    dij = matRad_mixModPreconditioner(dij);
+end
 %Get Bounds
 
 if ~isfield(pln.propOpt,'boundMU')
@@ -384,13 +387,18 @@
 wOpt = optimizer.wResult;
 info = optimizer.resultInfo;
 bxidx = 1;
+
 for mod = 1: pln.numOfModalities
+
     wt = [];
     % split the w for current modality
     STrepmat = (~dij.spatioTemp(mod) + dij.spatioTemp(mod)*dij.numOfSTscen(mod));
     wt = reshape(wOpt(bxidx: bxidx+STrepmat*dij.original_Dijs{mod}.totalNumOfBixels-1),[dij.original_Dijs{mod}.totalNumOfBixels,STrepmat]);
 
     resultGUI{mod} = matRad_calcCubes(wt,dij.original_Dijs{mod});
+    if isfield(dij,'preconW') && dij.precon
+       wt = wt.*dij.preconW(mod);     
+    end
     resultGUI{mod}.wUnsequenced = wt;
     resultGUI{mod}.usedOptimizer = optimizer;
     resultGUI{mod}.info = info;

matRad_mixModPreconditioner.m

@@ -0,0 +1,47 @@
+function dij = matRad_mixModPreconditioner(dij) 
+% Dose influence preconditioner for mixed modality plans
+%
+%
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% 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.
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% check if its is a mixed mod dij exists 
+if ~isfield(dij,'original_Dijs')
+    return;
+end
+
+% check to identify what is the preconditioning value 
+
+for i = 1: numel(dij.original_Dijs)
+    preconW(i) = max(mean(dij.original_Dijs{i}.physicalDose{1},1));
+end
+preconW = round(preconW./min(preconW));
+dij.preconW = 1./preconW;
+
+% make change to all dij type structures in one 
+for mod = 1 : numel(dij.original_Dijs)
+    fieldNames = fieldnames(dij.original_Dijs{mod});
+
+    for i = 1 : numel(fieldNames)
+        if iscell(dij.original_Dijs{mod}.(fieldNames{i}))
+
+            for j = 1 : numel(dij.original_Dijs{mod}.(fieldNames{i})) 
+                dij.original_Dijs{mod}.(fieldNames{i}){j} = bsxfun(@times,dij.original_Dijs{mod}.(fieldNames{i}){j}, dij.preconW(mod)); 
+            end
+        end
+    end
+end
+
+end