Added new NLP optimizer - goal is to make it robust and capable of taking in multiple ROI requirements.

Compare_optPlans.m

@@ -0,0 +1,24 @@
+
+
+function Compare_optPlans
+    
+
+
+    load('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\Geometry.mat')
+    load('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\NLP_result.mat')
+    load('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\optResults.mat')
+    
+    roi_idx = 1;
+    nfrac = 1;
+    dose_WP = optResults.dose{end};
+    dose_NLP = NLP_result.dose;
+
+    [dvh_WP dosebins_WP] = dvhist(dose_WP,Geometry,roi_idx,nfrac);
+    [dvh_NLP dosebins_NLP] = dvhist(dose_NLP,Geometry,roi_idx,nfrac);
+    
+    figure
+    hold on
+    plot(dosebins_WP, dvh_WP,'Color', [0.9,0.2,0.2],'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+    plot(dosebins_NLP, dvh_NLP,'Color', [0.2,0.9,0.2],'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+    hold off
+end

NLP_beamlet_optimizer.m

@@ -0,0 +1,103 @@
+
+% this is a simple example for a constrained NLP solver from matlabs
+% optimization toolbox
+
+function [D_full, w_fin] = NLP_beamlet_optimizer
+fprintf('starting NLP optimization process')
+
+% load the data
+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
+targetROI_ind = 1;
+numVox  = numel(Geometry.data);
+numBeam = size(beamlet_cell_array,2);
+
+ROI_idx=Geometry.ROIS{1, 2}.ind;
+fprintf('.')
+
+%% create input for optimization
+% ideal target values
+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));
+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;
+    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));
+
+% initial 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);
+fun = @(x) eval_f(x, beamlets_pruned, target_pruned);
+
+
+% all weights*-1 should be < 0 (b)
+A = [];
+b = [];
+Aeq = [];
+beq = [];
+lb = zeros(1, numBeam);
+ub = [];
+nonlcon = [];
+
+% define optimizer options
+options = optimoptions('fmincon');
+options.MaxFunctionEvaluations = 100000;
+options.Display = 'iter';
+
+fprintf('--')
+%% Run the optimization
+tic
+x = 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;
+
+%% 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');
+end
+
+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

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -49,10 +49,9 @@ ypmax = 1.0;    % +jaw width / 2
 % y-prime points in the z-direction in the CT coordinate system
 
 % Number of beamlets in the BEV for each direction
-Mxp = 7;        % number of leaves; leaf size is 1/20cm= 0.5mm
-Nyp = 1;        % always 1 for Tomo due to binary mlc
-% 
-% Mxp = 20;        % ## Grozomah
+Mxp = 7;  % Mxp = 20;  number of leaves; leaf size is 1/20cm= 0.5mm
+Nyp = 1;  % always 1 for Tomo due to binary mlc
+
 
 % ============================================= End of user-supplied inputs
 executable_path = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\WiscPlanEXE\RyanCsphoton.x86.exe';

WiscPlanPhotonkV125/matlab_frontend/lab/XTPS.m

@@ -474,8 +474,8 @@ classdef XTPS < handle
         %-------------------------------------------------------------------------------
         % Grozomah
         function MergeBeamlets_Callback(obj, src, evt)
-            disp('Test 1 works!')
             merge_beamlets(4, obj.patient_dir);
+            disp('Beamlets merged!')
         end
         %-------------------------------------------------------------------------------
         % Grozomah