forked from e0404/matRad
-
Notifications
You must be signed in to change notification settings - Fork 0
/
matRad_fluenceOptimization.m
206 lines (168 loc) · 7.92 KB
/
matRad_fluenceOptimization.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
function [resultGUI,info] = matRad_fluenceOptimization(dij,cst,pln)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% matRad inverse planning wrapper function
%
% call
% [resultGUI,info] = 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.
% info: struct containing information about optimization
%
% 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.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% issue warning if biological optimization impossible
if sum(strcmp(pln.bioOptimization,{'LEMIV_effect','LEMIV_RBExD'}))>0 && (~isfield(dij,'mAlphaDose') || ~isfield(dij,'mSqrtBetaDose')) && strcmp(pln.radiationMode,'carbon')
warndlg('Alpha and beta matrices for effect based and RBE optimization not available - physical optimization is carried out instead.');
pln.bioOptimization = 'none';
end
if ~isdeployed % only if _not_ running as standalone
% add path for optimization functions
matRadRootDir = fileparts(mfilename('fullpath'));
addpath(fullfile(matRadRootDir,'optimization'))
% get handle to Matlab command window
mde = com.mathworks.mde.desk.MLDesktop.getInstance;
cw = mde.getClient('Command Window');
xCmdWndView = cw.getComponent(0).getViewport.getComponent(0);
h_cw = handle(xCmdWndView,'CallbackProperties');
% set Key Pressed Callback of Matlab command window
set(h_cw, 'KeyPressedCallback', @matRad_CWKeyPressedCallback);
end
% initialize global variables for optimizer
global matRad_global_x;
global matRad_global_d;
global matRad_STRG_C_Pressed;
global matRad_objective_function_value;
matRad_global_x = NaN * ones(dij.totalNumOfBixels,1);
matRad_global_d = NaN * ones(dij.numOfVoxels,1);
matRad_STRG_C_Pressed = false;
matRad_objective_function_value = [];
% consider VOI priorities
cst = matRad_setOverlapPriorities(cst);
% adjust objectives and constraints internally for fractionation
for i = 1:size(cst,1)
for j = 1:size(cst{i,6},1)
cst{i,6}(j).dose = cst{i,6}(j).dose/pln.numOfFractions;
end
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}];
doseTarget = [doseTarget cst{i,6}.dose];
ixTarget = [ixTarget i*ones(1,length([cst{i,6}.dose]))];
end
end
[doseTarget,i] = max(doseTarget);
ixTarget = ixTarget(i);
wOnes = ones(dij.totalNumOfBixels,1);
% set the IPOPT options.
matRad_ipoptOptions;
% modified settings for photon dao
if pln.runDAO && strcmp(pln.radiationMode,'photons')
% options.ipopt.max_iter = 50;
% options.ipopt.acceptable_obj_change_tol = 7e-3; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled
end
% set bounds on optimization variables
options.lb = zeros(1,dij.totalNumOfBixels); % Lower bound on the variables.
options.ub = inf * ones(1,dij.totalNumOfBixels); % Upper bound on the variables.
funcs.iterfunc = @(iter,objective,paramter) matRad_IpoptIterFunc(iter,objective,paramter,options.ipopt.max_iter);
% calculate initial beam intensities wInit
if strcmp(pln.bioOptimization,'const_RBExD') && 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
bixelWeight = (doseTarget)/(dij.RBE * mean(dij.physicalDose{1}(V,:)*wOnes));
wInit = wOnes * bixelWeight;
elseif (strcmp(pln.bioOptimization,'LEMIV_effect') || strcmp(pln.bioOptimization,'LEMIV_RBExD')) ...
&& strcmp(pln.radiationMode,'carbon')
% check if you are running a supported rad
dij.ax = zeros(dij.numOfVoxels,1);
dij.bx = zeros(dij.numOfVoxels,1);
for i = 1:size(cst,1)
if isequal(cst{i,3},'OAR') || isequal(cst{i,3},'TARGET')
dij.ax(cst{i,4}{1}) = cst{i,5}.alphaX;
dij.bx(cst{i,4}{1}) = cst{i,5}.betaX;
end
for j = 1:size(cst{i,6},2)
% check if prescribed doses are in a valid domain
if cst{i,6}(j).dose > 5 && isequal(cst{i,3},'TARGET')
error('Reference dose > 5Gy[RBE] for target. Biological optimization outside the valid domain of the base data. Reduce dose prescription or use more fractions.');
end
end
end
dij.ixDose = dij.bx~=0;
if isequal(pln.bioOptimization,'LEMIV_effect')
effectTarget = cst{ixTarget,5}.alphaX * doseTarget + cst{ixTarget,5}.betaX * doseTarget^2;
p = (sum(dij.mAlphaDose{1}(V,:)*wOnes)) / (sum((dij.mSqrtBetaDose{1}(V,:) * wOnes).^2));
q = -(effectTarget * length(V)) / (sum((dij.mSqrtBetaDose{1}(V,:) * wOnes).^2));
wInit = -(p/2) + sqrt((p^2)/4 -q) * wOnes;
elseif isequal(pln.bioOptimization,'LEMIV_RBExD')
%pre-calculations
dij.gamma = zeros(dij.numOfVoxels,1);
dij.gamma(dij.ixDose) = dij.ax(dij.ixDose)./(2*dij.bx(dij.ixDose));
% calculate current in target
CurrEffectTarget = (dij.mAlphaDose{1}(V,:)*wOnes + (dij.mSqrtBetaDose{1}(V,:)*wOnes).^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*(dij.physicalDose{1}(V,:)*wOnes)));
wInit = ((doseTarget)/(TolEstBio*maxCurrRBE*max(dij.physicalDose{1}(V,:)*wOnes)))* wOnes;
end
else
bixelWeight = (doseTarget)/(mean(dij.physicalDose{1}(V,:)*wOnes));
wInit = wOnes * bixelWeight;
pln.bioOptimization = 'none';
end
% set optimization options
options.radMod = pln.radiationMode;
options.bioOpt = pln.bioOptimization;
options.ID = [pln.radiationMode '_' pln.bioOptimization];
options.numOfScenarios = dij.numOfScenarios;
% set callback functions.
[options.cl,options.cu] = matRad_getConstBoundsWrapper(cst,options);
funcs.objective = @(x) matRad_objFuncWrapper(x,dij,cst,options);
funcs.constraints = @(x) matRad_constFuncWrapper(x,dij,cst,options);
funcs.gradient = @(x) matRad_gradFuncWrapper(x,dij,cst,options);
funcs.jacobian = @(x) matRad_jacobFuncWrapper(x,dij,cst,options);
funcs.jacobianstructure = @( ) matRad_getJacobStruct(dij,cst);
% Run IPOPT.
[wOpt, info] = ipopt(wInit,funcs,options);
% calc dose and reshape from 1D vector to 2D array
fprintf('Calculating final cubes...\n');
resultGUI = matRad_calcCubes(wOpt,dij,cst);
resultGUI.wUnsequenced = wOpt;
% unset Key Pressed Callback of Matlab command window
if ~isdeployed
set(h_cw, 'KeyPressedCallback',' ');
end
% clear global variables
clearvars -global matRad_global_x matRad_global_d matRad_objective_function_value matRad_STRG_C_Pressed;
% unblock mex files
clear mex