1
0

NLP_beamlet_optimizer.m 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366
  1. % [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
  2. % orthoslice(D_full, [0,70])
  3. % colorwash(Geometry.data, D_full, [500, 1500], [0,70])
  4. function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
  5. % This function performs the beamlet optimization
  6. % Inputs: none. Still needs to find "Geometry" and "beamlets" from Wiscplan
  7. % Outputs: full dose image dose: D_full, optimal beamlet weights: w_fin
  8. %
  9. % Made by Peter Ferjancic 1. May 2018
  10. % Last updated: 14. August 2018
  11. % Inspired by Ana Barrigan's REGGUI optimization procedures
  12. N_fcallback1 = 5000;
  13. N_fcallback2 = 200000;
  14. patient = 'gbm_005';
  15. switch patient
  16. case 'patient'
  17. patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
  18. case 'tomoPhantom'
  19. patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
  20. case 'phantom_HD'
  21. patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom';
  22. case 'doggo'
  23. patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_dog5_3';
  24. case 'gbm_005'
  25. patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
  26. otherwise
  27. error('invalid case')
  28. end
  29. merge_beamlets(4, patient_dir);
  30. %% PROGRAM STARTS HERE
  31. % - no tocar lo que hay debajo -
  32. fprintf('starting NLP optimization process... ')
  33. % -- LOAD GEOMETRY AND BEAMLETS --
  34. load([patient_dir '\matlab_files\Geometry.mat'])
  35. % beamlet_batch_filename = [patient_dir '\beamlet_batch_files\' 'beamletbatch0.bin'];
  36. beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
  37. beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
  38. fprintf('\n loaded beamlets...')
  39. % -- SET INITIAL PARAMETERS --
  40. numVox = numel(Geometry.data);
  41. numBeamlet = size(beamlet_cell_array,2);
  42. %% -- BEAMLET DOSE DELIVERY --
  43. beamlets = get_beamlets(beamlet_cell_array, numVox);
  44. % show_joint_beamlets(beamlets, size(Geometry.data), 7:9)
  45. fprintf('\n beamlet array constructed...')
  46. % - merge beamlets into beams -
  47. load([patient_dir '\all_beams.mat'])
  48. beamletOrigin=[0 0 0];
  49. beam_i=0;
  50. beam_i_list=[];
  51. for beamlet_i = 1:numel(all_beams)
  52. newBeamletOrigin = all_beams{beamlet_i}.ip;
  53. if any(newBeamletOrigin ~= beamletOrigin)
  54. beam_i = beam_i+1;
  55. beamletOrigin = newBeamletOrigin;
  56. end
  57. beam_i_list=[beam_i_list, beam_i];
  58. end
  59. beamlets_joined=[];
  60. target_joined=[];
  61. wbar2 = waitbar(0, 'merging beamlets into beams');
  62. numBeam=numel(unique(beam_i_list));
  63. for beam_i = unique(beam_i_list)
  64. beam_full = sum(beamlets(:,beam_i_list == beam_i), 2);
  65. beamlets_joined(:,end+1) = beam_full;
  66. waitbar(beam_i/numBeam, wbar2)
  67. end
  68. close(wbar2)
  69. %% -- OPTIMIZATION TARGETS --
  70. make_ROI_goals(Geometry, beamlets, beamlets_joined, patient);
  71. [optGoal, optGoal_beam, optGoal_idx, targetMinMax_idx] = get_ROI_goals(patient);
  72. % -- make them robust --
  73. RO_params=0;
  74. optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
  75. optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  76. %% -- INITIALIZE BEAMLET WEIGHTS --
  77. w0_beams = ones(numBeam,1);
  78. w0_beams = mean(optGoal_beam{1}.target ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
  79. w0_beams = w0_beams + (2*rand(size(w0_beams))-1) *0.1 .*w0_beams; % random perturbation
  80. % -- CALLBACK OPTIMIZATION FUNCTION --
  81. fun1 = @(x) get_penalty(x, optGoal_beam);
  82. fun2 = @(x) get_penalty(x, optGoal);
  83. % -- OPTIMIZATION PARAMETERS --
  84. % define optimization parameters
  85. A = [];
  86. b = [];
  87. Aeq = [];
  88. beq = [];
  89. lb = zeros(1, numBeamlet);
  90. lb_beam = zeros(1, numBeamlet);
  91. ub = [];
  92. nonlcon = [];
  93. % define opt limits, and make it fmincon progress
  94. options = optimoptions('fmincon');
  95. options.MaxFunctionEvaluations = N_fcallback1;
  96. options.Display = 'iter';
  97. options.PlotFcn = 'optimplotfval';
  98. % options.UseParallel = true;
  99. fprintf('\n running initial optimizer:')
  100. %% Run the optimization
  101. % -- GET FULL BEAM WEIGHTS --
  102. tic
  103. w_beam = fmincon(fun1,w0_beams,A,b,Aeq,beq,lb_beam,ub,nonlcon,options);
  104. % t=toc;
  105. % disp(['Optimization time for beams = ',num2str(t)]);
  106. w_beamlets = ones(numBeamlet,1);
  107. numBeam=numel(unique(beam_i_list));
  108. for beam_i = 1:numBeam % assign weights to beamlets
  109. % beamlets from same beam get same initial weights
  110. w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
  111. end
  112. w_beamlets = w_beamlets + (2*rand(size(w_beamlets))-1) *0.1 .*w_beamlets; % small random perturbation
  113. % -- GET FULL BEAMLET WEIGHTS --
  114. options.MaxFunctionEvaluations = N_fcallback2;
  115. % tic
  116. w_fin = fmincon(fun2,w_beamlets,A,b,Aeq,beq,lb,ub,nonlcon,options);
  117. t=toc;
  118. disp(['Optimization time for beamlets = ',num2str(t)]);
  119. %% evaluate the results
  120. D_full = reshape(beamlets * w_fin, size(Geometry.data));
  121. %% save outputs
  122. NLP_result.dose = D_full;
  123. NLP_result.weights = w_fin;
  124. save([patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');
  125. plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
  126. colorwash(Geometry.data, D_full);
  127. % plot_DVH_robust(D_full, optGoal, optGoal_idx)
  128. end
  129. %% support functions
  130. % ---- PENALTY FUNCTION ----
  131. function penalty = get_penalty(x, optGoal)
  132. % this function gets called by the optimizer. It checks the penalty for
  133. % all the robust implementation and returns the worst result.
  134. NumScenarios = optGoal{1}.NbrRandScenarios * optGoal{1}.NbrSystSetUpScenarios * optGoal{1}.NbrRangeScenarios;
  135. fobj = zeros(NumScenarios,1);
  136. sc_i = 1;
  137. for nrs_i = 1:optGoal{1}.NbrRandScenarios
  138. for sss_i = 1 :optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  139. for rrs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
  140. fobj(sc_i)=eval_f(x, optGoal, nrs_i, sss_i, rrs_i);
  141. sc_i = sc_i + 1;
  142. end
  143. end
  144. end
  145. % take the worst case penalty of evaluated scenarios
  146. penalty=max(fobj);
  147. end
  148. % ------ supp: penalty for single scenario ------
  149. function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
  150. penalty = 0;
  151. % for each condition
  152. for goal_i = 1:numel(optGoal)
  153. switch optGoal{goal_i}.function
  154. % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
  155. case 'min'
  156. % penalize if achieved dose is lower than target dose
  157. d_penalty = 1.0e0 * sum(max(0, ...
  158. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
  159. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)));
  160. case 'max'
  161. % penalize if achieved dose is higher than target dose
  162. d_penalty = 1.0e0 * sum(max(0, ...
  163. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)-...
  164. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target)));
  165. case 'min_sq'
  166. % penalize if achieved dose is higher than target dose
  167. temp1=min(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)-...
  168. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target));
  169. d_penalty = 1.0e0 * sum(temp1.^2);
  170. case 'max_sq'
  171. % penalize if achieved dose is higher than target dose
  172. temp1=max(0, (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x)-...
  173. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target));
  174. d_penalty = 1.0e0 * sum(temp1.^2);
  175. case 'LeastSquare'
  176. % penalize with sum of squares any deviation from target
  177. % dose
  178. d_penalty = 1.0e-1* sum(((...
  179. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) - ...
  180. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target).^2);
  181. case 'min_perc_Volume'
  182. % penalize by amount of volume under threshold
  183. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
  184. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) > 0)) ...
  185. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
  186. d_penalty = 3.0e5 * min(perc_vox-0.05, 0)
  187. case 'max_perc_Volume'
  188. % penalize by amount of volume under threshold
  189. perc_vox = numel(find((optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target) -...
  190. (optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned * x) < 0)) ...
  191. / numel(optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target);
  192. d_penalty = 3.0e4 * min(perc_vox-0.05, 0)
  193. end
  194. penalty = penalty + d_penalty * optGoal{goal_i}.opt_weight;
  195. end
  196. end
  197. % ---- GET BEAMLETS ----
  198. function beamlets = get_beamlets(beamlet_cell_array, numVox);
  199. wbar1 = waitbar(0, 'Creating beamlet array');
  200. numBeam = size(beamlet_cell_array,2);
  201. batchSize=100;
  202. beamlets = sparse(0, 0);
  203. for beam_i=1:numBeam
  204. % for each beam define how much dose it delivers on each voxel
  205. idx=beamlet_cell_array{1, beam_i}.non_zero_indices;
  206. % break the beamlets into multiple batches
  207. if rem(beam_i, batchSize)==1;
  208. beamlet_batch = sparse(numVox, batchSize);
  209. beam_i_temp=1;
  210. end
  211. beamlet_batch(idx, beam_i_temp) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
  212. waitbar(beam_i/numBeam, wbar1, ['Adding beamlet array: #', num2str(beam_i)])
  213. % add the batch to full set when filled
  214. if rem(beam_i, batchSize)==0;
  215. beamlets =[beamlets, beamlet_batch];
  216. end
  217. % crop and add the batch to full set when completed
  218. if beam_i==numBeam;
  219. beamlet_batch=beamlet_batch(:, 1:beam_i_temp);
  220. beamlets =[beamlets, beamlet_batch];
  221. end
  222. beam_i_temp=beam_i_temp+1;
  223. end
  224. close(wbar1)
  225. end
  226. function show_joint_beamlets(beamlets, IMGsize, Beam_list)
  227. % this function overlays and plots multiple beamlets. The goal is to
  228. % check whether some beamlets are part of the same beam manually.
  229. beam=sum(beamlets(:, Beam_list), 2);
  230. beamImg=reshape(full(beam), IMGsize);
  231. orthoslice(beamImg)
  232. end
  233. % ---- MAKE ROI ROBUST ----
  234. function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
  235. % take regular optimal goal and translate it into several robust cases
  236. % RO_params - should have the information below
  237. % nrs - random scenarios
  238. % sss - system setup scenarios
  239. % rrs - random range scenarios
  240. % X - X>0 moves image right
  241. % Y - Y>0 moves image down
  242. % Z - in/out.
  243. shift_mag = 1; % vox of shift
  244. nrs_scene_list={[0,0,0]};
  245. % sss_scene_list={[0,0,0]};
  246. sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0]};
  247. % sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
  248. % [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
  249. rrs_scene_list={[0,0,0]};
  250. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_DP_target.nrrd');
  251. % [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
  252. for i = 1:numel(optGoal)
  253. optGoal{i}.NbrRandScenarios =numel(nrs_scene_list);
  254. optGoal{i}.NbrSystSetUpScenarios=numel(sss_scene_list);
  255. optGoal{i}.NbrRangeScenarios =numel(rrs_scene_list);
  256. end
  257. for goal_i = 1:numel(optGoal)
  258. % get target
  259. idx=optGoal{goal_i}.ROI_idx;
  260. targetImg1=zeros(optGoal{goal_i}.imgDim);
  261. targetImg1(idx)=1;
  262. % get beamlets
  263. for nrs_i = 1:optGoal{goal_i}.NbrRandScenarios % num. of random scenarios
  264. % modify target and beamlets
  265. targetImg2=targetImg1;
  266. % beamlets stay the same
  267. for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios % syst. setup scenarios = sss
  268. % modify target and beamlets
  269. targetImg3=get_RO_sss(targetImg2, sss_scene_list{sss_i});
  270. % beamlets stay the same
  271. for rrs_i = 1:optGoal{goal_i}.NbrRangeScenarios % range scenario = rrs
  272. % modify target and beamlets
  273. targetImg4=targetImg3;
  274. % beamlets stay the same
  275. %% make new target and beamlets
  276. ROI_idx=[];
  277. ROI_idx=find(targetImg4>0);
  278. if isfield(optGoal{goal_i}, 'target_alpha')
  279. target = double(optGoal{goal_i}.target_alpha * targetIn(ROI_idx));
  280. disp('exists')
  281. else
  282. target = ones(numel(ROI_idx), 1) * optGoal{goal_i}.D_final;
  283. disp('not exists')
  284. end
  285. beamlets_pruned = beamlets(ROI_idx, :);
  286. % save to optGoal output
  287. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.ROI_idx = ROI_idx;
  288. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.beamlets_pruned = beamlets_pruned;
  289. optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.target = target;
  290. end
  291. end
  292. end
  293. end
  294. end
  295. % ------ supp: RO case SSS ------
  296. function targetImg3=get_RO_sss(targetImg2, sss_scene_shift);
  297. % translate the target image
  298. targetImg3 = imtranslate(targetImg2,sss_scene_shift);
  299. end