123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103 |
- function [D_full, w_fin] = NLP_beamlet_optimizer
- fprintf('starting NLP optimization process')
- patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
- load([patient_dir '\matlab_files\Geometry.mat'])
- beamlet_batch_filename = [patient_dir '\beamlet_batch_files\' 'beamletbatch0.bin'];
- beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
- fprintf('.')
- targetROI_ind = 1;
- numVox = numel(Geometry.data);
- numBeam = size(beamlet_cell_array,2);
- ROI_idx=Geometry.ROIS{1, 2}.ind;
- fprintf('.')
- target = sparse(zeros(numVox,1));
- target(ROI_idx) = 60;
- beamlets = sparse(zeros(numVox, numBeam));
- fprintf('.')
- wbar1 = waitbar(0, 'Creating beamlet array');
- for beam_i=1:numBeam
-
- idx=beamlet_cell_array{1, beam_i}.non_zero_indices;
- beamlets(idx, beam_i) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
-
- waitbar(beam_i/numBeam)
- end
- close(wbar1)
- fprintf('.')
- pruneIdx= find(target~=0);
- beamlets_pruned= beamlets(pruneIdx, :);
- target_pruned = full(target(pruneIdx));
- w0 = ones(numBeam,1);
- w0 = mean(target_pruned./ (beamlets_pruned*w0)) * w0;
- fun = @(x) eval_f(x, beamlets_pruned, target_pruned);
- A = [];
- b = [];
- Aeq = [];
- beq = [];
- lb = zeros(1, numBeam);
- ub = [];
- nonlcon = [];
- options = optimoptions('fmincon');
- options.MaxFunctionEvaluations = 100000;
- options.Display = 'iter';
- fprintf('--')
- tic
- x = fmincon(fun,w0,A,b,Aeq,beq,lb,ub,nonlcon,options);
- t=toc;
- disp(['Optimization run time = ',num2str(t)]);
- D_full = reshape(beamlets * x, size(Geometry.data));
- w_fin=x;
- NLP_result.dose = D_full;
- NLP_result.weights = w_fin;
- save('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\NLP_result.mat', 'NLP_result');
- end
- function penalty = eval_f(x, beamlets, target)
-
-
- penalty = 0;
- penalty = penalty + sum(max(0, target-(beamlets * x)));
-
- penalty = penalty + sum(max(0, (beamlets * x)-(target+5)));
- end
|