NLP_optimizer_v2.m 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387
  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_v2(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, beamlet_weights)
  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. if nargin<2
  24. load('WiscPlan_preferences.mat')
  25. [Pat_path] = uigetdir([WiscPlan_preferences.patientDataPath ], 'Select Patient folder');
  26. [Goal_file,Goal_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select OptGoal file');
  27. path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
  28. path2goal = [Goal_path, Goal_file];
  29. else
  30. Pat_path = varargin{1};
  31. path2geometry = [Pat_path, '\matlab_files\Geometry.mat'];
  32. path2goal = varargin{2};
  33. end
  34. str = inputdlg({'N of iterations for initial calc', 'N of iterations for full calc', ...
  35. 'Use pre-existing NLP_result to initiate? (y/n)'}, 'input', [1,35], {'10000', '500000', 'n'});
  36. N_fcallback1 = str2double(str{1}); % 10000 is a good guesstimate
  37. N_fcallback2 = str2double(str{2}); % 500000 is a good guesstimate
  38. pre_beamWeights = str{3};
  39. switch pre_beamWeights
  40. case 'y'
  41. [NLP_file,NLP_path,indx] = uigetfile([Pat_path '\matlab_files\*.mat'], 'Select NLP_result file');
  42. load([NLP_path, NLP_file])
  43. w_beamlets = NLP_result.weights;
  44. load([Pat_path, '\all_beams.mat'])
  45. if numel(all_beams) ~= numel(w_beamlets)
  46. error('Provided weight number does not match beamlet number!')
  47. end
  48. case 'n'
  49. disp('Initial beam weights will be calculated.')
  50. end
  51. %% PROGRAM STARTS HERE
  52. % - no tocar lo que hay debajo -
  53. fprintf('starting NLP optimization process... \n')
  54. % % -- LOAD GEOMETRY, GOALS, BEAMLETS --
  55. load(path2geometry)
  56. load(path2goal)
  57. [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, Pat_path);
  58. %% -- OPTIMIZATION TARGETS --
  59. % -- make the optimization optGoal structure --
  60. for i_goal = 1:size(OptGoals.goals,1)
  61. answer = inputdlg(['# of supervoxels for "' OptGoals.data{i_goal}.name '" with ' num2str(numel(OptGoals.data{i_goal}.ROI_idx)) ' vox: ("0" to skip)'])
  62. switch answer{1}
  63. case '0'
  64. % if not supervoxel, just select provided ROI_idx
  65. optGoal{i_goal} = OptGoals.data{i_goal};
  66. optGoal{i_goal}.beamlets_pruned = sparse(beamlets(optGoal{i_goal}.ROI_idx, :));
  67. optGoal_beam{i_goal} = OptGoals.data{i_goal};
  68. optGoal_beam{i_goal}.beamlets_pruned = sparse(beamlets_joined(optGoal{i_goal}.ROI_idx, :));
  69. otherwise
  70. % -- if supervoxel, merge given columns
  71. % - make supervoxel map
  72. mask = zeros(OptGoals.data{i_goal}.imgDim);
  73. mask(OptGoals.data{i_goal}.ROI_idx) = 1;
  74. superMask = superpix_group(mask, str2double(answer{1}));
  75. superVoxList = unique(superMask);
  76. superVoxList = superVoxList(superVoxList>0);
  77. optGoal{i_goal} = OptGoals.data{i_goal};
  78. optGoal{i_goal}.ROI_idx_old = optGoal{i_goal}.ROI_idx; % copy old index data
  79. optGoal{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
  80. optGoal{i_goal}.opt_weight = optGoal{i_goal}.opt_weight * numel(optGoal{i_goal}.ROI_idx_old)/numel(optGoal{i_goal}.ROI_idx);
  81. optGoal_beam{i_goal} = OptGoals.data{i_goal};
  82. optGoal_beam{i_goal}.ROI_idx_old = optGoal_beam{i_goal}.ROI_idx;
  83. optGoal_beam{i_goal}.ROI_idx = zeros(numel(superVoxList), 1);
  84. optGoal_beam{i_goal}.opt_weight = optGoal_beam{i_goal}.opt_weight * numel(optGoal_beam{i_goal}.ROI_idx_old)/numel(optGoal_beam{i_goal}.ROI_idx);
  85. h_w1 = waitbar(0, 'merging superboxels');
  86. for i_supVox = 1:numel(superVoxList)
  87. supVox_idx = superVoxList(i_supVox);
  88. waitbar(i_supVox/numel(superVoxList), h_w1)
  89. idxList = find(superMask == supVox_idx);
  90. optGoal{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets(idxList, :),1));
  91. optGoal_beam{i_goal}.beamlets_pruned(i_supVox,:) = sparse(mean(beamlets_joined(idxList, :),1));
  92. % -- make new indeces
  93. optGoal{i_goal}.ROI_idx(i_supVox) = idxList(1);
  94. optGoal_beam{i_goal}.ROI_idx(i_supVox) = idxList(1);
  95. end
  96. close(h_w1)
  97. end
  98. end
  99. % -- make them robust --
  100. RO_params=0;
  101. optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
  102. optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  103. % -- CALLBACK OPTIMIZATION FUNCTION --
  104. fun1 = @(x) get_penalty(x, optGoal_beam);
  105. fun2 = @(x) get_penalty(x, optGoal);
  106. % -- OPTIMIZATION PARAMETERS --
  107. % define optimization parameters
  108. A = [];
  109. b = [];
  110. Aeq = [];
  111. beq = [];
  112. lb = zeros(1, numBeamlet);
  113. lb_beam = zeros(1, numBeam);
  114. ub = [];
  115. nonlcon = [];
  116. % define opt limits, and make it fmincon progress
  117. options = optimoptions('fmincon');
  118. options.MaxFunctionEvaluations = N_fcallback1;
  119. options.Display = 'iter';
  120. options.PlotFcn = 'optimplotfval';
  121. % options.UseParallel = true;
  122. options.UseParallel = false;
  123. % options.OptimalityTolerance = 1e-9;
  124. %% -- INITIALIZE BEAMLET WEIGHTS --
  125. switch pre_beamWeights
  126. case 'y'
  127. % should have been assigned previously.
  128. disp('Provided beamlet weights used for initial comparison')
  129. case 'n'
  130. % if initial beamlet weights are not provided, get quick estimate
  131. fprintf('\n running initial optimizer:')
  132. % initialize beamlet weights, OR
  133. w0_beams = ones(numBeam,1);
  134. w0_beams = mean(optGoal_beam{1}.D_final(optGoal{1}.ROI_idx) ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
  135. w0_beams = double(w0_beams);
  136. % -- GET BEAM WEIGHTS --
  137. tic
  138. w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb_beam,ub,nonlcon,options);
  139. fprintf(' done!:')
  140. t=toc;
  141. disp(['Optimization time for beams = ',num2str(t)]);
  142. w_beamlets = ones(numBeamlet,1);
  143. numBeam=numel(unique(beam_i_list));
  144. for beam_i = 1:numBeam % assign weights to beamlets
  145. % beamlets from same beam get same initial weights
  146. w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
  147. end
  148. end
  149. %% FULL OPTIMIZATION
  150. % -- GET FULL BEAMLET WEIGHTS --
  151. options.MaxFunctionEvaluations = N_fcallback2;
  152. tic
  153. fprintf('\n running full optimizer:')
  154. w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
  155. fprintf(' done!:')
  156. t=toc;
  157. disp(['Optimization time for beamlets = ',num2str(t)]);
  158. %% evaluate the results
  159. D_full = reshape(beamlets * w_fin, size(Geometry.data));
  160. %% save outputs
  161. NLP_result.dose = D_full;
  162. NLP_result.weights = w_fin;
  163. save([Pat_path, '\matlab_files\NLP_result.mat'], 'NLP_result');
  164. plot_DVH(Geometry, D_full)
  165. colorwash(Geometry.data, D_full, [-500, 500], [0, 35]);
  166. end
  167. %% support functions
  168. % ---- PENALTY FUNCTION ----
  169. function penalty = get_penalty(x, optGoal)
  170. % this function gets called by the optimizer. It checks the penalty for
  171. % all the robust implementation and returns the worst result.
  172. NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
  173. fobj = zeros(NumScenarios,1);
  174. sc_i = 1;
  175. for nrs_i = 1:optGoal{1}.NbrRandScenarios
  176. for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  177. for rgs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
  178. fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i);
  179. sc_i = sc_i + 1;
  180. end
  181. end
  182. end
  183. % take the worst case penalty of evaluated scenarios
  184. penalty=max(fobj);
  185. end
  186. % ------ supp: penalty for single scenario ------
  187. function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
  188. penalty = 0;
  189. % for each condition
  190. for goal_i = 1:numel(optGoal)
  191. switch optGoal{goal_i}.function
  192. % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
  193. case 'min'
  194. % penalize if achieved dose is lower than target dose
  195. d_penalty = 1.0e0 * sum(max(0, ...
  196. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  197. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)));
  198. case 'max'
  199. % penalize if achieved dose is higher than target dose
  200. d_penalty = 1.0e0 * sum(max(0, ...
  201. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  202. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target)));
  203. case 'min_sq'
  204. % penalize if achieved dose is lower than target dose
  205. temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  206. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  207. d_penalty = 1.0e0 * sum(temp1.*temp1);
  208. case 'max_sq'
  209. % penalize if achieved dose is higher than target dose
  210. temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  211. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  212. d_penalty = 1.0e0 * sum(temp1.*temp1);
  213. case 'LeastSquare'
  214. % penalize with sum of squares any deviation from target
  215. % dose
  216. temp1 = (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) - ...
  217. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target;
  218. d_penalty = 1.0e0* sum(temp1.^2);
  219. case 'min_perc_Volume'
  220. % penalize by amount of volume under threshold
  221. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  222. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) > 0)) ...
  223. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
  224. d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
  225. case 'max_perc_Volume'
  226. % penalize by amount of volume under threshold
  227. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  228. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) < 0)) ...
  229. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
  230. d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
  231. end
  232. penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
  233. end
  234. end
  235. % ---- MAKE ROI ROBUST ----
  236. function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  237. % take regular optimal goal and translate it into several robust cases
  238. % RO_params - should have the information below
  239. % nrs - random scenarios
  240. % sss - system setup scenarios
  241. % rgs - random range scenarios
  242. % X - X>0 moves image right
  243. % Y - Y>0 moves image down
  244. % Z - in/out.
  245. shift_mag = 1; % vox of shift
  246. nrs_scene_list={[0,0,0]};
  247. % ----====#### CHANGE ROBUSTNESS HERE ####====----
  248. sss_scene_list={[0,0,0]};
  249. % sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0], [0,0,-1], [0,0,1]};
  250. % sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
  251. % [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
  252. % ----====#### CHANGE ROBUSTNESS HERE ####====----
  253. rgs_scene_list={[0,0,0]};
  254. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
  255. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
  256. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
  257. for i = 1:numel(optGoal)
  258. optGoal{i}.NbrRandScenarios =numel(nrs_scene_list);
  259. optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
  260. optGoal{i}.NbrRangeScenarios =numel(rgs_scene_list);
  261. end
  262. for goal_i = 1:numel(optGoal)
  263. % get target
  264. idx=optGoal{goal_i}.ROI_idx;
  265. targetImg1=zeros(optGoal{goal_i}.imgDim);
  266. targetImg1(idx)=1;
  267. % get beamlets
  268. for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios % num. of random scenarios
  269. % modify target and beamlets
  270. targetImg2=targetImg1;
  271. % beamlets stay the same
  272. for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  273. % modify target and beamlets
  274. [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
  275. % beamlets stay the same
  276. for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios % range scenario = rgs
  277. % modify target and beamlets
  278. targetImg4=targetImg3;
  279. % beamlets stay the same
  280. %% make new target and beamlets
  281. ROI_idx=[];
  282. ROI_idx=find(targetImg4>0);
  283. target = optGoal{goal_i}.D_final(idxValid);
  284. beamlets_pruned = beamlets(ROI_idx, :);
  285. % save to optGoal output
  286. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx = ROI_idx;
  287. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned = beamlets_pruned;
  288. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target = target;
  289. end
  290. end
  291. end
  292. end
  293. end
  294. % ------ supp: RO case SSS ------
  295. function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
  296. % translate the target image
  297. targetImg3 = imtranslate(targetImg2,sss_scene_shift);
  298. % now we need to figure out if any target voxels fell out during the
  299. % shift
  300. imgValid = imtranslate(targetImg3,-sss_scene_shift);
  301. imgInvalid = (targetImg2-imgValid);
  302. idx_1 = find(targetImg2);
  303. idx_2 = find(imgInvalid);
  304. [idxValid,ia] = setdiff(idx_1,idx_2);
  305. [C,ia, ib] = intersect(idx_1,idxValid);
  306. end