NLP_beamlet_optimizer.m 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340
  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_beamlet_optimizer(paths)
  8. % This function performs the beamlet optimization
  9. % Inputs: none. Still needs to find "Geometry" and "beamlets" from Wiscplan
  10. % Outputs: full dose image dose: D_full, optimal beamlet weights: w_fin
  11. %
  12. % [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
  13. %
  14. % Made by Peter Ferjancic 1. May 2018
  15. % Last updated: 14. August 2018
  16. % Inspired by Ana Barrigan's REGGUI optimization procedures
  17. %
  18. % paths.patient_dir
  19. % paths.Goal_dir (previously called DP_dir)
  20. % paths.patient
  21. % paths.goalsName
  22. N_fcallback1 = 5000;
  23. N_fcallback2 = 200000;
  24. % paths.patient = 'medivation_01';
  25. % % paths.goalsName = 'ROI_goals_simple';
  26. if false
  27. 'medivation_01'
  28. paths.patient_dir = 'F:\021_WiscPlan_data\Medivation_0391811_crop1';
  29. paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2';
  30. paths.Goal_dir = '\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2';
  31. paths.goalsName = 'ROI_goals_simple';
  32. 'temp'
  33. paths.patient_dir = 'C:\000-temp\DELETEME';
  34. paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis'
  35. 'avastin_009'
  36. paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
  37. paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
  38. 'gbm_005'
  39. % paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
  40. paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam_2';
  41. paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
  42. 'gbm_015'
  43. paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
  44. paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
  45. 'gbm_022'
  46. paths.patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_022';
  47. paths.DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL022\B1\Processed';
  48. 'avastin_009'
  49. error('invalid case')
  50. end
  51. % merge_beamlets(4, patient_dir);
  52. %% PROGRAM STARTS HERE
  53. % - no tocar lo que hay debajo -
  54. fprintf('starting NLP optimization process... \n')
  55. % % -- LOAD GEOMETRY AND BEAMLETS --
  56. load([paths.patient_dir '\matlab_files\Geometry.mat'])
  57. %% -- OPTIMIZATION TARGETS --
  58. load([paths.Goal_dir '\RODP_files\' paths.goalsName '.mat']);
  59. optGoal = ROI_goals.optGoal;
  60. optGoal_beam = ROI_goals.optGoal_beam;
  61. optGoal_idx = ROI_goals.optGoal_idx;
  62. targetMinMax_idx = ROI_goals.targetMinMax_idx;
  63. [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, paths.patient_dir);
  64. % -- make them robust --
  65. RO_params=0;
  66. optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
  67. optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  68. %% -- INITIALIZE BEAMLET WEIGHTS --
  69. w0_beams = ones(numBeam,1);
  70. w0_beams = mean(optGoal_beam{1}.target ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
  71. w0_beams = w0_beams + (2*rand(size(w0_beams))-1) *0.1 .*w0_beams; % random perturbation
  72. w0_beams = double(w0_beams);
  73. % -- CALLBACK OPTIMIZATION FUNCTION --
  74. fun1 = @(x) get_penalty(x, optGoal_beam);
  75. fun2 = @(x) get_penalty(x, optGoal);
  76. % -- OPTIMIZATION PARAMETERS --
  77. % define optimization parameters
  78. A = [];
  79. b = [];
  80. Aeq = [];
  81. beq = [];
  82. lb = zeros(1, numBeamlet);
  83. lb_beam = zeros(1, numBeam);
  84. ub = [];
  85. nonlcon = [];
  86. % define opt limits, and make it fmincon progress
  87. options = optimoptions('fmincon');
  88. options.MaxFunctionEvaluations = N_fcallback1;
  89. options.Display = 'iter';
  90. options.PlotFcn = 'optimplotfval';
  91. % options.UseParallel = true;
  92. % options.OptimalityTolerance = 1e-9;
  93. fprintf('\n running initial optimizer:')
  94. %% Run the optimization
  95. % -- GET FULL BEAM WEIGHTS --
  96. tic
  97. w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb_beam,ub,nonlcon,options);
  98. fprintf(' done!:')
  99. % t=toc;
  100. % disp(['Optimization time for beams = ',num2str(t)]);
  101. w_beamlets = ones(numBeamlet,1);
  102. numBeam=numel(unique(beam_i_list));
  103. for beam_i = 1:numBeam % assign weights to beamlets
  104. % beamlets from same beam get same initial weights
  105. w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
  106. end
  107. w_beamlets = w_beamlets + (2*rand(size(w_beamlets))-1) *0.1 .*w_beamlets; % small random perturbation
  108. w_beamlets = 1.1* w_beamlets; % this just kicks the beamlets up a bit to make it easier for the optimizer to start
  109. % -- GET FULL BEAMLET WEIGHTS --
  110. options.MaxFunctionEvaluations = N_fcallback2;
  111. % tic
  112. fprintf('\n running full optimizer:')
  113. w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
  114. fprintf(' done!:')
  115. t=toc;
  116. disp(['Optimization time for beamlets = ',num2str(t)]);
  117. %% evaluate the results
  118. D_full = reshape(beamlets * w_fin, size(Geometry.data));
  119. %% save outputs
  120. NLP_result.dose = D_full;
  121. NLP_result.weights = w_fin;
  122. save([paths.patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');
  123. plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
  124. colorwash(Geometry.data, D_full, [-500, 500], [0, 110]);
  125. % plot_DVH_robust(D_full, optGoal, optGoal_idx)
  126. end
  127. %% support functions
  128. % ---- PENALTY FUNCTION ----
  129. function penalty = get_penalty(x, optGoal)
  130. % this function gets called by the optimizer. It checks the penalty for
  131. % all the robust implementation and returns the worst result.
  132. NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
  133. fobj = zeros(NumScenarios,1);
  134. sc_i = 1;
  135. for nrs_i = 1:optGoal{1}.NbrRandScenarios
  136. for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  137. for rgs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
  138. fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rgs_i);
  139. sc_i = sc_i + 1;
  140. end
  141. end
  142. end
  143. % take the worst case penalty of evaluated scenarios
  144. penalty=max(fobj);
  145. end
  146. % ------ supp: penalty for single scenario ------
  147. function penalty = eval_f(x, optGoal, nrs_i, sss_i, rgs_i)
  148. penalty = 0;
  149. % for each condition
  150. for goal_i = 1:numel(optGoal)
  151. switch optGoal{goal_i}.function
  152. % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
  153. case 'min'
  154. % penalize if achieved dose is lower than target dose
  155. d_penalty = 1.0e0 * sum(max(0, ...
  156. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  157. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)));
  158. case 'max'
  159. % penalize if achieved dose is higher than target dose
  160. d_penalty = 1.0e0 * sum(max(0, ...
  161. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  162. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target)));
  163. case 'min_sq'
  164. % penalize if achieved dose is higher than target dose
  165. temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  166. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  167. d_penalty = 1.0e0 * sum(temp1.^2);
  168. case 'max_sq'
  169. % penalize if achieved dose is higher than target dose
  170. temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x)-...
  171. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target));
  172. d_penalty = 1.0e0 * sum(temp1.^2);
  173. case 'LeastSquare'
  174. % penalize with sum of squares any deviation from target
  175. % dose
  176. d_penalty = 1.0e-1* sum(((...
  177. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) - ...
  178. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target).^2);
  179. case 'min_perc_Volume'
  180. % penalize by amount of volume under threshold
  181. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  182. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) > 0)) ...
  183. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
  184. d_penalty = 3.0e5 * min(perc_vox-0.05, 0)
  185. case 'max_perc_Volume'
  186. % penalize by amount of volume under threshold
  187. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target) -...
  188. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned * x) < 0)) ...
  189. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target);
  190. d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
  191. end
  192. penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
  193. end
  194. end
  195. % ---- MAKE ROI ROBUST ----
  196. function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  197. % take regular optimal goal and translate it into several robust cases
  198. % RO_params - should have the information below
  199. % nrs - random scenarios
  200. % sss - system setup scenarios
  201. % rgs - random range scenarios
  202. % X - X>0 moves image right
  203. % Y - Y>0 moves image down
  204. % Z - in/out.
  205. shift_mag = 2; % vox of shift
  206. nrs_scene_list={[0,0,0]};
  207. % ----====#### CHANGE ROBUSTNESS HERE ####====----
  208. % sss_scene_list={[0,0,0]};
  209. 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]};
  210. % sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
  211. % [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
  212. % ----====#### CHANGE ROBUSTNESS HERE ####====----
  213. rgs_scene_list={[0,0,0]};
  214. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
  215. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
  216. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
  217. for i = 1:numel(optGoal)
  218. optGoal{i}.NbrRandScenarios =numel(nrs_scene_list);
  219. optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
  220. optGoal{i}.NbrRangeScenarios =numel(rgs_scene_list);
  221. end
  222. for goal_i = 1:numel(optGoal)
  223. % get target
  224. idx=optGoal{goal_i}.ROI_idx;
  225. targetImg1=zeros(optGoal{goal_i}.imgDim);
  226. targetImg1(idx)=1;
  227. % get beamlets
  228. for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios % num. of random scenarios
  229. % modify target and beamlets
  230. targetImg2=targetImg1;
  231. % beamlets stay the same
  232. for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  233. % modify target and beamlets
  234. [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
  235. % beamlets stay the same
  236. for rgs_i = 1:optGoal{goal_i}.NbrRangeScenarios % range scenario = rgs
  237. % modify target and beamlets
  238. targetImg4=targetImg3;
  239. % beamlets stay the same
  240. %% make new target and beamlets
  241. ROI_idx=[];
  242. ROI_idx=find(targetImg4>0);
  243. if isfield(optGoal{goal_i}, 'target_alpha')
  244. target = optGoal{goal_i}.target;
  245. target = target(idxValid);
  246. % disp('exists')
  247. else
  248. target = ones(numel(ROI_idx), 1) .* optGoal{goal_i}.D_final;
  249. % disp('not exists')
  250. end
  251. beamlets_pruned = beamlets(ROI_idx, :);
  252. % save to optGoal output
  253. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.ROI_idx = ROI_idx;
  254. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.beamlets_pruned = beamlets_pruned;
  255. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rgs{rgs_i}.target = target;
  256. end
  257. end
  258. end
  259. end
  260. end
  261. % ------ supp: RO case SSS ------
  262. function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
  263. % translate the target image
  264. targetImg3 = imtranslate(targetImg2,sss_scene_shift);
  265. % now we need to figure out if any target voxels fell out during the
  266. % shift
  267. imgValid = imtranslate(targetImg3,-sss_scene_shift);
  268. imgInvalid = (targetImg2-imgValid);
  269. idx_1 = find(targetImg2);
  270. idx_2 = find(imgInvalid);
  271. [idxValid,ia] = setdiff(idx_1,idx_2);
  272. [C,ia, ib] = intersect(idx_1,idxValid);
  273. end