Added optimization tracking (penalty function plot). Preparation before changing to intake ROI structs for optimization parameters.

NLP_beamlet_optimizer.m

@@ -3,16 +3,29 @@
 % optimization toolbox
 
 function [D_full, w_fin] = NLP_beamlet_optimizer
-fprintf('starting NLP optimization process')
+% This function performs the beamlet optimization
+% Inputs: none. Still needs to find "Geometry" and "beamlets" from Wiscplan
+% Outputs: full dose image dose: D_full, optimal beamlet weights: w_fin
+%
+% Made by Peter Ferjancic May 2nd 2018
+% Inspired by Ana Barrigan's REGGUI optimization procedures
 
-% load the data
+% To-do:
+% - add functionality for multiple ROI optimization/multiple goals
+% - Add robusness aspect (+take worst case scenario, see REGGUI)
+
+
+%% Start code here
+fprintf('starting NLP optimization process: ')
+
+% -- LOAD GEOMETRY AND BEAMLETS --
 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('.')
 
-% set initial parameters
+% -- SET INITIAL PARAMETERS --
 targetROI_ind = 1;
 numVox  = numel(Geometry.data);
 numBeam = size(beamlet_cell_array,2);
@@ -21,15 +34,14 @@ ROI_idx=Geometry.ROIS{1, 2}.ind;
 fprintf('.')
 
 %% create input for optimization
-% ideal target values
+% -- OPTIMIZATION TARGETS --
 target = sparse(zeros(numVox,1));
 target(ROI_idx) = 60; % 60 units of dose to the target, 0 otherwise
 
-% create the map on how beamlets deliver the dose
-beamlets = sparse(zeros(numVox, numBeam));
+% -- BEAMLET DOSE DELIVERY --
+beamlets = sparse(numVox, numBeam);
 fprintf('.')
 wbar1 = waitbar(0, 'Creating beamlet array');
-
 for beam_i=1:numBeam
     % for each beam define how much dose it delivers on each voxel
     idx=beamlet_cell_array{1, beam_i}.non_zero_indices;
@@ -39,22 +51,19 @@ for beam_i=1:numBeam
 end
 close(wbar1)
 fprintf('.')
-
 % prune the beamlets matrix only to areas of interest
 pruneIdx= find(target~=0);
 beamlets_pruned= beamlets(pruneIdx, :);
 target_pruned = full(target(pruneIdx));
 
-% initial beamlet weights
+% -- INITIALIZE BEAMLET WEIGHTS --
 w0 = ones(numBeam,1);
 w0 = mean(target_pruned./ (beamlets_pruned*w0)) * w0;
 
-% optimization function - least squares
-% fun = @(x) sum(((beamlets_pruned * x) - target_pruned).^2);
+% -- CALLBACK OPTIMIZATION FUNCTION --
 fun = @(x) eval_f(x, beamlets_pruned, target_pruned);
 
-
-% all weights*-1 should be < 0 (b)
+% -- OPTIMIZATION PARAMETERS --
 A = [];
 b = [];
 Aeq = [];
@@ -65,10 +74,11 @@ nonlcon = [];
 
 % define optimizer options
 options = optimoptions('fmincon');
-options.MaxFunctionEvaluations = 100000;
+options.MaxFunctionEvaluations = 50000;
 options.Display = 'iter';
+options.PlotFcn = 'optimplotfval';
 
-fprintf('--')
+fprintf('-')
 %% Run the optimization
 tic
 x = fmincon(fun,w0,A,b,Aeq,beq,lb,ub,nonlcon,options);
@@ -78,8 +88,6 @@ disp(['Optimization run time = ',num2str(t)]);
 %% evaluate the results
 D_full = reshape(beamlets * x, size(Geometry.data));
 
-
-
 w_fin=x;
 
 %% save outputs
@@ -89,8 +97,25 @@ NLP_result.weights = w_fin;
 % orthoslice(NLP_result.dose, [0,70])
 
 save('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\NLP_result.mat', 'NLP_result');
+
+plot_DVH(D_full, Geometry, 2, 1)
+
 end
 
+
+%% support functions
+% ---- PENALTY FUNCTION ----
+function penalty = eval_f_new(x, beamlets, target)
+    % sum squares
+%     penalty = sum(((beamlets * x) - target).^2);
+
+    % at least>target
+    penalty = 0;
+    penalty = penalty + sum(max(0, target-(beamlets * x)));
+    
+    penalty = penalty + sum(max(0, (beamlets * x)-(target+5)));
+end
+% ---- the old penalty function ----
 function penalty = eval_f(x, beamlets, target)
     % sum squares
 %     penalty = sum(((beamlets * x) - target).^2);
@@ -101,3 +126,14 @@ function penalty = eval_f(x, beamlets, target)
     
     penalty = penalty + sum(max(0, (beamlets * x)-(target+5)));
 end
+% ---- DVH PLOT FUNCTION ----
+function plot_DVH(dose, Geometry, roi_idx, nfrac)
+    % this function plots the DVH of the given dose
+    [dvh dosebins] = dvhist(dose,Geometry,roi_idx,nfrac);
+    figure
+    plot(dosebins, dvh,'Color', [0.9,0.2,0.2],'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+end
+
+
+
+