1
0

NLP_optimizer_v3.m 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569
  1. % paths.patient_dir
  2. % paths.Goal_dir (previously called DP_dir)
  3. % paths.patient
  4. % paths.goalsName
  5. % colorwash(Geometry.data, D_full, [500, 1500], [0,70])
  6. % orthoslice(D_full, [0,70])
  7. function [D_full, w_fin, Geometry, optGoal] = NLP_optimizer_v3(varargin)
  8. % This function performs the beamlet optimization
  9. % [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
  10. %
  11. % Inputs:
  12. % () OR
  13. % (Pat_path, path2goal) OR
  14. % (Pat_path, path2goal, path2NLP_result)
  15. % Pat_path, path2goal = strings to patient folder and optimal goals
  16. % beamlet_weights = initial beamlet weights
  17. %
  18. % Outputs:
  19. % full dose image dose: [D_full, w_fin, Geometry, optGoal]
  20. %
  21. % Made by Peter Ferjancic 1. May 2018
  22. % Last updated: 1. April 2019
  23. dialogue_box = 'no';
  24. if nargin<2
  25. load('WiscPlan_preferences.mat')
  26. [Pat_path] = uigetdir([WiscPlan_preferences.patientDataPath ], 'Select Patient folder');
  27. [Goal_file,Goal_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select OptGoal file');
  28. path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
  29. path2goal = [Goal_path, Goal_file];
  30. elseif nargin ==2
  31. Pat_path = varargin{1};
  32. path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
  33. path2goal = varargin{2};
  34. [Goal_path,Goal_file,ext] = fileparts(path2goal);
  35. else
  36. Pat_path = varargin{1};
  37. path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
  38. path2goal = varargin{2};
  39. [Goal_path,Goal_file,ext] = fileparts(path2goal);
  40. path2NLPres = varargin{3};
  41. load(path2NLPres);
  42. dialogue_box = 'pass';
  43. pre_beamWeights = 'yS';
  44. end
  45. switch dialogue_box
  46. case 'yes'
  47. str = inputdlg({'N of iterations for initial calc', 'N of iterations for full calc', ...
  48. 'Use pre-existing NLP_result to initiate? (y/n)'}, 'input', [1,35], {'100000', '500000', 'n'});
  49. N_fcallback1 = str2double(str{1}); % 100000 is a good guesstimate
  50. N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
  51. pre_beamWeights = str{3};
  52. case 'no'
  53. disp('dialogue box skipped')
  54. N_fcallback1 = 1e5;
  55. N_fcallback2 = 2e6; % 500000;
  56. pre_beamWeights = 'n';
  57. case 'pass'
  58. N_fcallback1 = 1e5;
  59. N_fcallback2 = 2e6; % 500000;
  60. end
  61. switch pre_beamWeights
  62. case 'y'
  63. [NLP_file,NLP_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select NLP_result file');
  64. load([NLP_path, NLP_file])
  65. w_beamlets = NLP_result.weights;
  66. load([Pat_path, '\all_beams.mat'])
  67. if numel(all_beams) ~= numel(w_beamlets)
  68. error('Provided weight number does not match beamlet number!')
  69. end
  70. case 'yS'
  71. w_beamlets = NLP_result.weights;
  72. case 'n'
  73. disp('Initial beam weights will be calculated.')
  74. end
  75. %% PROGRAM STARTS HERE
  76. % - no tocar lo que hay debajo -
  77. fprintf('starting NLP optimization process... \n')
  78. % % -- LOAD GEOMETRY, GOALS, BEAMLETS --
  79. load(path2geometry)
  80. load(path2goal)
  81. [beamlets, numBeamlet] = get_beamlets(Geometry, Pat_path, OptGoals);
  82. % [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, Pat_path);
  83. %% -- OPTIMIZATION TARGETS --
  84. % -- make the optimization optGoal structure --
  85. for i_goal = 1:size(OptGoals.goals,1)
  86. if isfield(OptGoals.data{i_goal}, 'SupVox_num')
  87. SupVox_num = OptGoals.data{i_goal}.SupVox_num;
  88. else
  89. answer = inputdlg(['# of supervoxels for "' OptGoals.data{i_goal}.name '" with ' num2str(numel(OptGoals.data{i_goal}.ROI_idx)) ' vox: ("0" to skip)'])
  90. SupVox_num = str2double(answer{1})
  91. end
  92. optGoal{i_goal} = OptGoals.data{i_goal};
  93. optGoal{i_goal}.sss_scene_list = OptGoals.sss_scene_list;
  94. optGoal{i_goal}.maxModulation = OptGoals.maxModulation;
  95. optGoal{i_goal}.BeamSmoothMax = OptGoals.BeamSmoothMax;
  96. optGoal{i_goal}.optFuncNum = OptGoals.optFuncNum;
  97. % modulation
  98. if strcmp( '-||-' , optGoal{i_goal}.ROI_name)
  99. optGoal{i_goal}.ROI_idx_old = optGoal{i_goal-1}.ROI_idx_old; % copy old index data
  100. optGoal{i_goal}.ROI_idx = optGoal{i_goal-1}.ROI_idx;
  101. optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
  102. optGoal{i_goal}.ROI_idx = optGoal{i_goal-1}.ROI_idx
  103. if isfield(OptGoals.data{i_goal}, 'wgt_map')
  104. optGoal{i_goal}.vox_wgt = optGoal{i_goal-1}.vox_wgt;
  105. end
  106. optGoal{i_goal}.beamlets_pruned = optGoal{i_goal-1}.beamlets_pruned;
  107. else
  108. switch SupVox_num
  109. case 0
  110. % if not supervoxel, just select provided ROI_idx
  111. optGoal{i_goal}.beamlets_pruned = sparse(beamlets(optGoal{i_goal}.ROI_idx, :));
  112. otherwise
  113. % -- if supervoxel, merge given columns
  114. % - make supervoxel map
  115. switch OptGoals.goals{i_goal, 3}
  116. case 'Fixed dose'
  117. mask = zeros(OptGoals.data{i_goal}.imgDim);
  118. mask(OptGoals.data{i_goal}.ROI_idx) = 1;
  119. % group superpixels
  120. superMask = superpix_group(mask, SupVox_num, 'no');
  121. case 'Dose map'
  122. mask = zeros(OptGoals.data{i_goal}.imgDim);
  123. mask(OptGoals.data{i_goal}.ROI_idx) = OptGoals.data{i_goal}.D_final;
  124. mask2 = round(mask/3)*3;
  125. superMask = superpix_group(mask, SupVox_num, 'no');
  126. orthoslice(superMask)
  127. end
  128. superVoxList = unique(superMask);
  129. superVoxList = superVoxList(superVoxList>0);
  130. optGoal{i_goal}.ROI_idx_old = optGoal{i_goal}.ROI_idx; % copy old index data
  131. optGoal{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
  132. optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
  133. if isfield(OptGoals.data{i_goal}, 'wgt_map')
  134. tabula_wgtmap = zeros(size(superMask));
  135. tabula_wgtmap(OptGoals.data{i_goal}.ROI_idx) = OptGoals.data{i_goal}.wgt_map;
  136. end
  137. h_w1 = waitbar(0, 'merging superboxels');
  138. for i_supVox = 1:numel(superVoxList)
  139. waitbar(i_supVox/numel(superVoxList), h_w1)
  140. supVox_idx = superVoxList(i_supVox);
  141. idxList = find(superMask == supVox_idx);
  142. optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
  143. if isfield(OptGoals.data{i_goal}, 'wgt_map')
  144. % optGoal{i_goal}.vox_wgt(i_supVox) = sum(tabula_wgtmap(idxList));
  145. optGoal{i_goal}.vox_wgt(i_supVox) = mean(tabula_wgtmap(idxList));
  146. end
  147. % -- make new indeces
  148. optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
  149. end
  150. close(h_w1)
  151. end
  152. end
  153. end
  154. % -- make them robust --
  155. RO_params=0;
  156. optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  157. % save([Goal_path, Goal_file '_robust.mat'], 'optGoal')
  158. % -- get beamlet indeces --
  159. load([Pat_path, '\all_beams.mat'])
  160. Nbeamlets = all_beams{1}.Mxp; % number of beamlets in a beam - usually 64 or 32
  161. weightTable = zeros(100,Nbeamlets);
  162. for ind_bmlt = 1:numel(all_beams)
  163. bLet_idx.y(ind_bmlt) = floor(all_beams{1, ind_bmlt}.num/Nbeamlets)+1;
  164. bLet_idx.x(ind_bmlt) = rem(all_beams{1, ind_bmlt}.num, Nbeamlets)+1;
  165. weightTable(bLet_idx.y(ind_bmlt),bLet_idx.x(ind_bmlt)) = 1;
  166. end
  167. bLet_idx.idx = find(weightTable>0);
  168. bLet_idx.Nbeamlets = Nbeamlets;
  169. disp('.')
  170. % -- CALLBACK OPTIMIZATION FUNCTION --
  171. fun1 = @(x) get_penalty(x, optGoal_beam, bLet_idx);
  172. fun2 = @(x) get_penalty(x, optGoal, bLet_idx);
  173. % -- OPTIMIZATION PARAMETERS --
  174. % define optimization parameters
  175. A = [];
  176. b = [];
  177. Aeq = [];
  178. beq = [];
  179. lb = zeros(1, numBeamlet);
  180. % lb_beam = zeros(1, numBeam);
  181. ub = [];
  182. nonlcon = [];
  183. % define opt limits, and make it fmincon progress
  184. options = optimoptions('fmincon');
  185. options.MaxFunctionEvaluations = N_fcallback1;
  186. options.Display = 'iter';
  187. options.PlotFcn = 'optimplotfval';
  188. options.MaxIterations = 100;
  189. % options.UseParallel = true;
  190. options.UseParallel = false;
  191. % options.OptimalityTolerance = 1e-9;
  192. %% -- INITIALIZE BEAMLET WEIGHTS --
  193. switch pre_beamWeights
  194. case 'y'
  195. % should have been assigned previously.
  196. disp('Provided beamlet weights used for initial comparison')
  197. case 'n'
  198. % if initial beamlet weights are not provided, get quick estimate
  199. % fprintf('\n running initial optimizer:')
  200. % initialize beamlet weights, OR
  201. w0 = ones(numBeamlet,1);
  202. % w0 = mean(optGoal{1}.D_final(optGoal{1}.ROI_idx) ./ (optGoal{1}.beamlets_pruned*w0+0.1)) * w0; % old
  203. w0 = mean(optGoal{1}.D_final(:)) ./ mean(optGoal{1}.beamlets_pruned*w0+0.05) * w0;
  204. w_beamlets = double(w0);
  205. % -- GET BEAM WEIGHTS --
  206. % tic
  207. % w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb,ub,nonlcon,options);
  208. % fprintf(' done!:')
  209. % t=toc;
  210. % disp(['Optimization time for beams = ',num2str(t)]);
  211. %
  212. % w_beamlets = ones(numBeamlet,1);
  213. % numBeam=numel(unique(beam_i_list));
  214. % for beam_i = 1:numBeam % assign weights to beamlets
  215. % % beamlets from same beam get same initial weights
  216. % w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
  217. % end
  218. end
  219. %% FULL OPTIMIZATION
  220. % -- GET FULL BEAMLET WEIGHTS --
  221. options.MaxFunctionEvaluations = N_fcallback2;
  222. tic
  223. fprintf('\n running full optimizer:')
  224. w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
  225. fprintf(' done!:')
  226. t=toc;
  227. disp(['Optimization time for beamlets = ',num2str(t)]);
  228. %% evaluate the results
  229. D_full = reshape(beamlets * w_fin, size(Geometry.data));
  230. %% save outputs
  231. NLP_result.dose = D_full;
  232. NLP_result.weights = w_fin;
  233. NLP_result.sss_scene_list = optGoal{1}.sss_scene_list;
  234. NLP_result.maxModulation = OptGoals.maxModulation;
  235. NLP_result.BeamSmoothMax = OptGoals.BeamSmoothMax;
  236. save([Pat_path, '\matlab_files\NLP_result_' Goal_file '.mat'], 'NLP_result');
  237. % plot_DVH(Geometry, D_full)
  238. DQM_DVH(Geometry, D_full)
  239. colorwash(Geometry.data, D_full, [500, 1500], [0, 90]);
  240. end
  241. %% support functions
  242. % ---- PENALTY FUNCTION ----
  243. function penalty = get_penalty(x, optGoal, bLet_idx)
  244. % this function gets called by the optimizer. It checks the penalty for
  245. % all the robust implementation and returns the worst result.
  246. NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
  247. fobj = zeros(NumScenarios,numel(optGoal)+2);
  248. sc_i = 1;
  249. for nrs_i = 1:optGoal{1}.NbrRandScenarios
  250. for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  251. for rgs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
  252. % fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx);
  253. fobj(sc_i,:) = eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx);
  254. sc_i = sc_i + 1;
  255. end
  256. end
  257. end
  258. switch optGoal{1}.optFuncNum
  259. case 1 % "RO max"
  260. penalty=max(sum(fobj,2));
  261. case 2 % "RO mean"
  262. penalty=mean(sum(fobj,2));
  263. case 3 % "RO objective-based"
  264. penalty=sum(max(fobj,[],1));
  265. case 4 % "Classical"
  266. penalty=sum(fobj(1,:));
  267. end
  268. % take the worst case penalty of evaluated scenarios
  269. % penalty=max(fobj);
  270. % take the median worst case (stochastic robust)
  271. % penalty=median(fobj);
  272. end
  273. % ------ supp: penalty for single scenario ------
  274. function penArray = eval_f(x, optGoal, nrs_i, sss_i, rgs_i, bLet_idx)
  275. % penalty = 0;
  276. penArray = zeros(1,numel(optGoal)+2); % all the optGoals, plus modulation plus smoothing
  277. % for each condition
  278. for goal_i = 1:numel(optGoal)
  279. switch optGoal{goal_i}.function
  280. % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
  281. case 'min'
  282. % penalize if achieved dose is lower than target dose
  283. d_penalty = 1.0e0 * sum(max(0, ...
  284. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  285. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)));
  286. case 'max'
  287. % penalize if achieved dose is higher than target dose
  288. d_penalty = 1.0e0 * sum(max(0, ...
  289. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  290. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target)));
  291. case 'min_sq'
  292. % penalize if achieved dose is lower than target dose
  293. temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  294. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  295. d_penalty = 1.0e0 * sum(temp1.*temp1);
  296. case 'max_sq'
  297. % penalize if achieved dose is higher than target dose
  298. temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  299. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  300. d_penalty = 1.0e0 * sum(temp1.*temp1);
  301. case 'min_step'
  302. % penalize if achieved dose is lower than target dose
  303. temp1=-min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  304. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  305. d_penalty = 1.0e1 * (sum(temp1.*temp1) + 1.0e3* sum(temp1>0)/numel(temp1));
  306. case 'max_step'
  307. % penalize if achieved dose is higher than target dose
  308. temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  309. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  310. d_penalty = 1.0e1 * (sum(temp1.*temp1) + 1.0e3* sum(temp1>0)/numel(temp1));
  311. case 'LeastSquare'
  312. % penalize with sum of squares any deviation from target
  313. % dose
  314. temp1 = (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) - ...
  315. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target;
  316. d_penalty = 1.0e0* sum(temp1.^2);
  317. case 'min_perc_Volume'
  318. % penalize by amount of volume under threshold
  319. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  320. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) > 0)) ...
  321. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
  322. d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
  323. case 'max_perc_Volume'
  324. % penalize by amount of volume under threshold
  325. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  326. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) < 0)) ...
  327. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
  328. d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
  329. case 'min_sq_voxwgt'
  330. % penalize if achieved dose is lower than target dose
  331. temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  332. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  333. d_penalty = 1.0e0 * sum(temp1.*temp1.* optGoal{goal_i}.vox_wgt');
  334. case 'max_sq_voxwgt'
  335. % penalize if achieved dose is lower than target dose
  336. temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  337. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  338. d_penalty = 1.0e0 * sum(temp1.*temp1.* optGoal{goal_i}.vox_wgt');
  339. end
  340. % penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
  341. penArray(goal_i) = d_penalty * optGoal{goal_i}.opt_weight;
  342. end
  343. %% add modulation penalty
  344. if true
  345. if isfield(optGoal{goal_i}, 'maxModulation')
  346. max_modulation = optGoal{goal_i}.maxModulation;
  347. else
  348. error('no Max modulation parameter enterd!')
  349. max_modulation = 5;
  350. end
  351. mod_pen_weight = 1.0e10;
  352. % calc the penalty
  353. mod_excess = max(0, x-max_modulation*mean(x));
  354. % mod_excess = max(0, x-max_modulation*mean(x(x>1)));
  355. mod_pen1 = mod_pen_weight*(sum(mod_excess) + numel(mod_excess)* any(mod_excess));
  356. % penalty = penalty + mod_pen1;
  357. penArray(end-1) = mod_pen1;
  358. end
  359. %% add penlty for single off beamlets - version 2
  360. if false
  361. if isfield(optGoal{goal_i}, 'BeamSmoothMax')
  362. BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
  363. else
  364. error('no Max beam smooth parameter enterd!')
  365. BeamSmoothMax = 1e6;
  366. end
  367. mod_pen_weight = BeamSmoothMax; %1.0e6
  368. x_down = zeros(size(x));
  369. x_down(2:end) = x(1:end-1);
  370. x_down2 = zeros(size(x));
  371. x_down2(3:end) = x(1:end-2);
  372. x_up = zeros(size(x));
  373. x_up(1:end-1) = x(2:end);
  374. x_up2 = zeros(size(x));
  375. x_up2(1:end-2) = x(3:end);
  376. mod_pen2 = mod_pen_weight*sum((x_down+x_up-2*x).^2 + 0.5*(x_down2+x_up2-2*x).^2)/(mean(x)^(2)*numel(x));
  377. % penalty = penalty + mod_pen2;
  378. penArray(end) = mod_pen2;
  379. end
  380. %% add penlty for single off beamlets - version 3
  381. if true
  382. if isfield(optGoal{goal_i}, 'BeamSmoothMax')
  383. BeamSmoothMax = optGoal{goal_i}.BeamSmoothMax;
  384. else
  385. error('no Max beam smooth parameter enterd!')
  386. BeamSmoothMax = 1e6;
  387. end
  388. mod_pen_weight = BeamSmoothMax; %1.0e6
  389. % make 2D beamletWeight map
  390. maxY = max(bLet_idx.y);
  391. tabula = zeros(maxY, bLet_idx.Nbeamlets);
  392. tabula(bLet_idx.idx) = x;
  393. % for each index calculate laplace
  394. myWeights = 1/12*[1,3,1; 1,-12,1; 1,3,1];
  395. i=2:maxY-1;
  396. j=2:bLet_idx.Nbeamlets-1;
  397. tabula2(i,j) = ...
  398. myWeights(1)*tabula(i-1,j-1) + myWeights(4)*tabula(i-1,j) + myWeights(7)*tabula(i-1,j+1)...
  399. +myWeights(2)*tabula(i,j-1) + myWeights(8)*tabula(i,j+1) ...
  400. +myWeights(3)*tabula(i+1,j-1) + myWeights(6)*tabula(i+1,j) + myWeights(9)*tabula(i+1,j+1)...
  401. +myWeights(5)*tabula(i,j);
  402. tabula2(1,j) = ...
  403. myWeights(2)*tabula(1,j-1) + myWeights(8)*tabula(1,j+1)...
  404. +myWeights(3)*tabula(2,j-1) + myWeights(6)*tabula(2,j) + myWeights(9)*tabula(2,j+1)...
  405. +myWeights(5)*tabula(1,j);
  406. tabula2(maxY,j) =...
  407. myWeights(1)*tabula(maxY-1,j-1) + myWeights(4)*tabula(maxY-1,j) + myWeights(7)*tabula(maxY-1,j+1)...
  408. +myWeights(2)*tabula(maxY,j-1) + myWeights(8)*tabula(maxY,j+1) ...
  409. +myWeights(5)*tabula(maxY,j);
  410. % make sum of squares
  411. mod_pen2 = mod_pen_weight*sum((tabula2(:)).^2)/(mean(x.^2)*numel(x));
  412. % penalty = penalty + mod_pen2;
  413. penArray(end) = mod_pen2;
  414. end
  415. end
  416. % ---- MAKE ROI ROBUST ----
  417. function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  418. % take regular optimal goal and translate it into several robust cases
  419. % RO_params - should have the information below
  420. % nrs - random scenarios
  421. % sss - system setup scenarios
  422. % rgs - random range scenarios
  423. % X - X>0 moves image right
  424. % Y - Y>0 moves image down
  425. % Z - in/out.
  426. shift_X = 2; % vox of shift
  427. shift_Y = 2; % vox of shift
  428. shift_Z = 1; % vox of shift
  429. nrs_scene_list={[0,0,0]};
  430. % ----====#### CHANGE ROBUSTNESS HERE ####====----
  431. if isfield(optGoal{1}, 'sss_scene_list')
  432. sss_scene_list = optGoal{1}.sss_scene_list;
  433. else
  434. error('New OptGoals should no longer reach this code!')
  435. sss_scene_list={[0,0,0], [-shift_Y,0,0], [shift_Y,0,0], [0,-shift_X,0], [0,shift_X,0], [0,0,-shift_Z], [0,0,shift_Z]};
  436. optGoal{1}.sss_scene_list = sss_scene_list;
  437. end
  438. % sss_scene_list={[0,0,0]};
  439. % ----====#### CHANGE ROBUSTNESS HERE ####====----
  440. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
  441. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
  442. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
  443. rgs_scene_list={[0,0,0]};
  444. for i = 1:numel(optGoal)
  445. optGoal{i}.NbrRandScenarios =numel(nrs_scene_list);
  446. optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
  447. optGoal{i}.NbrRangeScenarios =numel(rgs_scene_list);
  448. end
  449. for goal_i = 1:numel(optGoal)
  450. % get target
  451. idx=optGoal{goal_i}.ROI_idx;
  452. targetImg1=zeros(optGoal{goal_i}.imgDim);
  453. targetImg1(idx)=1;
  454. % get beamlets
  455. for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios % num. of random scenarios
  456. % modify target and beamlets
  457. targetImg2=targetImg1;
  458. % beamlets stay the same
  459. for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  460. % modify target and beamlets
  461. [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
  462. % beamlets stay the same
  463. for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios % range scenario = rgs
  464. % modify target and beamlets
  465. targetImg4=targetImg3;
  466. % beamlets stay the same
  467. %% make new target and beamlets
  468. ROI_idx=[];
  469. ROI_idx=find(targetImg4>0);
  470. target = optGoal{goal_i}.D_final(idxValid);
  471. beamlets_pruned = beamlets(ROI_idx, :);
  472. % save to optGoal output
  473. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx = ROI_idx;
  474. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned = beamlets_pruned;
  475. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target = target;
  476. end
  477. end
  478. end
  479. end
  480. end
  481. %% ------ supp: RO case SSS ------
  482. function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
  483. % translate the target image
  484. targetImg3 = imtranslate(targetImg2,sss_scene_shift);
  485. % now we need to figure out if any target voxels fell out during the
  486. % shift
  487. imgValid = imtranslate(targetImg3,-sss_scene_shift);
  488. imgInvalid = (targetImg2-imgValid);
  489. idx_1 = find(targetImg2);
  490. idx_2 = find(imgInvalid);
  491. [idxValid,ia] = setdiff(idx_1,idx_2);
  492. [C,ia, ib] = intersect(idx_1,idxValid);
  493. end