added DVH plotting, cleaned up the code a bit.

NLP_beamlet_optimizer.m

@@ -11,11 +11,10 @@ function [D_full, w_fin] = NLP_beamlet_optimizer
 % Inspired by Ana Barrigan's REGGUI optimization procedures
 
 % To-do:
-% - add functionality for multiple ROI optimization/multiple goals
 % - Add robusness aspect (+take worst case scenario, see REGGUI)
 
 
-%% Start code here
+%% Program starts here
 fprintf('starting NLP optimization process: ')
 
 % -- LOAD GEOMETRY AND BEAMLETS --
@@ -26,18 +25,13 @@ beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
 fprintf('.')
 
 % -- SET INITIAL PARAMETERS --
-targetROI_ind = 1;
 numVox  = numel(Geometry.data);
 numBeam = size(beamlet_cell_array,2);
 
-ROI_idx=Geometry.ROIS{1, 2}.ind;
+% ROI_idx=Geometry.ROIS{1, 2}.ind;
 fprintf('.')
 
 %% create input for optimization
-% -- OPTIMIZATION TARGETS --
-target = sparse(zeros(numVox,1));
-target(ROI_idx) = 60; % 60 units of dose to the target, 0 otherwise
-
 % -- BEAMLET DOSE DELIVERY --
 beamlets = sparse(numVox, numBeam);
 fprintf('.')
@@ -46,24 +40,23 @@ 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;
     beamlets(idx, beam_i) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
-    
     waitbar(beam_i/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));
+
+% -- OPTIMIZATION TARGETS --
+optGoal = make_ROI_goals(Geometry, beamlets);
 
 % -- INITIALIZE BEAMLET WEIGHTS --
 w0 = ones(numBeam,1);
-w0 = mean(target_pruned./ (beamlets_pruned*w0)) * w0;
+w0 = mean(optGoal{1}.target ./ (optGoal{1}.beamlets_pruned*w0)) * w0;
 
 % -- CALLBACK OPTIMIZATION FUNCTION --
-fun = @(x) eval_f(x, beamlets_pruned, target_pruned);
+fun = @(x) eval_f(x, optGoal);
 
 % -- OPTIMIZATION PARAMETERS --
+% define optimization parameters
 A = [];
 b = [];
 Aeq = [];
@@ -72,68 +65,139 @@ lb = zeros(1, numBeam);
 ub = [];
 nonlcon = [];
 
-% define optimizer options
+% define opt limits, and make it fmincon progress
 options = optimoptions('fmincon');
-options.MaxFunctionEvaluations = 50000;
+options.MaxFunctionEvaluations = 100000;
 options.Display = 'iter';
 options.PlotFcn = 'optimplotfval';
-
 fprintf('-')
+
 %% Run the optimization
 tic
-x = fmincon(fun,w0,A,b,Aeq,beq,lb,ub,nonlcon,options);
+w_fin = fmincon(fun,w0,A,b,Aeq,beq,lb,ub,nonlcon,options);
 t=toc;
 disp(['Optimization run time = ',num2str(t)]);
 
 %% evaluate the results
-D_full = reshape(beamlets * x, size(Geometry.data));
-
-w_fin=x;
+D_full = reshape(beamlets * w_fin, size(Geometry.data));
 
 %% save outputs
-
 NLP_result.dose = D_full;
 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)
+% plot_DVH(D_full, Geometry, 2, 1)
+optGoal_idx=[1,3];
+plot_DVH(D_full, optGoal, optGoal_idx)
 
 end
 
+% orthoslice(D_full, [0,70])
 
 %% support functions
 % ---- PENALTY FUNCTION ----
-function penalty = eval_f_new(x, beamlets, target)
-    % sum squares
-%     penalty = sum(((beamlets * x) - target).^2);
-
-    % at least>target
+function penalty = eval_f(x, optGoal)
     penalty = 0;
-    penalty = penalty + sum(max(0, target-(beamlets * x)));
-    
-    penalty = penalty + sum(max(0, (beamlets * x)-(target+5)));
+    % for each condition
+    for i = 1:numel(optGoal)
+        switch optGoal{i}.function
+            case 'min'
+                % penalize if achieved dose is lower than target dose
+                d_penalty = sum(max(0, optGoal{i}.target-(optGoal{i}.beamlets_pruned * x)));
+            case 'max'
+                % penalize if achieved dose is higher than target dose
+                d_penalty = sum(max(0, (optGoal{i}.beamlets_pruned * x)-(optGoal{i}.target+5)));
+            case 'LeastSquare'
+                % penalize with sum of squares any deviation from target
+                % dose
+                d_penalty = sum(((beamlets * x) - target).^2);
+        end
+        penalty = penalty + d_penalty * optGoal{i}.opt_weight;
+    end
 end
-% ---- the old penalty function ----
-function penalty = eval_f(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
 % ---- 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);
+function plot_DVH(dose, optGoal, optGoal_idx)
+    % this function plots the DVHs of the given dose
+    nfrac=1;
+    
     figure
-    plot(dosebins, dvh,'Color', [0.9,0.2,0.2],'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+    hold on
+    for roi_idx=optGoal_idx
+        % plot the histogram
+        
+        [dvh dosebins] = get_dvh_data(dose,optGoal{roi_idx}.ROI_idx,nfrac);
+        plot(dosebins, dvh,'Color', optGoal{roi_idx}.dvh_col,'LineStyle', '-','DisplayName', optGoal{roi_idx}.ROI_name);
+    end
+    hold off
+    
+end
+function [dvh dosebins] = get_dvh_data(dose,roi_idx,nfrac);
+    % this function calculates the data for the DVH
+    dosevec = dose(roi_idx);
+    
+    
+    [pdf dosebins] = hist(dosevec, 999);
+    % clip negative bins
+    pdf = pdf(dosebins >= 0);
+    dosebins = dosebins(dosebins >= 0);
+    dvh = fliplr(cumsum(fliplr(pdf))) / numel(dosevec) * 100;
+    % append the last bin
+    dosebins = [dosebins dosebins(end)+0.1];
+    dvh = [dvh 0];
+    
 end
 
+% ---- MAKE ROI GOALS ----
+function optGoal = make_ROI_goals(Geometry, beamlets)
+    optGoal={};
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'Target_min';
+    goal_1.ROI_name = Geometry.ROIS{1, 2}.name;
+    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_1.ROI_idx = ROI_idx;
+    goal_1.D_final = 60;
+    goal_1.function = 'min';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target   = ones(numel(ROI_idx), 1) * goal_1.D_final;
+    goal_1.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
+    goal_1.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_1;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_2.name = 'Target_max';
+    goal_2.ROI_name = Geometry.ROIS{1, 2}.name;
+    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_2.ROI_idx = ROI_idx;
+    goal_2.D_final = 65;
+    goal_2.function = 'max';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
+    goal_2.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
+    goal_2.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_2;
+    % -- END DEFINITION OF GOAL --
+    
+    % -- START DEFINITION OF GOAL --
+    goal_3.name = 'ROI 1_max';
+    goal_3.ROI_name = Geometry.ROIS{1, 3}.name;
+    ROI_idx = Geometry.ROIS{1, 3}.ind;
+    goal_3.ROI_idx = ROI_idx;
+    goal_3.D_final = 0;
+    goal_3.function = 'max';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
+    goal_3.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
+    % assign target
+    optGoal{end+1}=goal_3;
+    % -- END DEFINITION OF GOAL --
+    
+end