Commit before overhauling NLP_beamlet_optimizer to exclude optGoal generation. That will get moved into a new function.

WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

@@ -28,9 +28,9 @@ switch patient
     case 'phantom_HD'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom';
     case 'doggo'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_dog5_3';
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\archive\PatientData_dog5_3';
     case 'gbm_005'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
     otherwise
         error('invalid case')
 end
@@ -96,7 +96,7 @@ optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
 w0_beams = ones(numBeam,1);
 w0_beams = mean(optGoal_beam{1}.target ./ (optGoal_beam{1}.beamlets_pruned*w0_beams+0.1)) * w0_beams;
 w0_beams = w0_beams + (2*rand(size(w0_beams))-1) *0.1 .*w0_beams; % random perturbation
-
+w0_beams = double(w0_beams);
 
 % -- CALLBACK OPTIMIZATION FUNCTION --
 fun1 = @(x) get_penalty(x, optGoal_beam);
@@ -109,7 +109,7 @@ b = [];
 Aeq = [];
 beq = [];
 lb = zeros(1, numBeamlet);
-lb_beam = zeros(1, numBeamlet);
+lb_beam = zeros(1, numBeam);
 ub = [];
 nonlcon = [];
 
@@ -119,6 +119,8 @@ options.MaxFunctionEvaluations = N_fcallback1;
 options.Display = 'iter';
 options.PlotFcn = 'optimplotfval';
 % options.UseParallel = true;
+% options.OptimalityTolerance = 1e-9;
+
 fprintf('\n running initial optimizer:')
 
 %% Run the optimization
@@ -136,6 +138,7 @@ for beam_i = 1:numBeam % assign weights to beamlets
     w_beamlets(beam_i_list == beam_i) = w_beam(beam_i);
 end
 w_beamlets = w_beamlets + (2*rand(size(w_beamlets))-1) *0.1 .*w_beamlets; % small random perturbation
+w_beamlets = 1.1* w_beamlets; % this just kicks the beamlets up a bit to make it easier for the optimizer to start
 
 % -- GET FULL BEAMLET WEIGHTS --
 options.MaxFunctionEvaluations = N_fcallback2;
@@ -293,7 +296,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % Y - Y>0 moves image down
     % Z - in/out.
     
-    shift_mag = 1; % vox of shift
+    shift_mag = 3; % vox of shift
     nrs_scene_list={[0,0,0]};
 
 %     sss_scene_list={[0,0,0]};
@@ -304,8 +307,9 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
 
     rrs_scene_list={[0,0,0]};
 
-%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\CDP_data\CDP5_DP_target.nrrd');
+    [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
 %     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
+%     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\archive\CDP_data\CDP5_DP_target.nrrd');
     
     for i = 1:numel(optGoal)
         optGoal{i}.NbrRandScenarios     =numel(nrs_scene_list);
@@ -344,7 +348,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
                         target = double(optGoal{goal_i}.target_alpha * targetIn(ROI_idx));
 %                         disp('exists')
                     else
-                        target = ones(numel(ROI_idx), 1) * optGoal{goal_i}.D_final;
+                        target = ones(numel(ROI_idx), 1) .* optGoal{goal_i}.D_final;
 %                         disp('not exists')
                     end
                     

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -57,7 +57,7 @@ load(geometry_file);
 ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
 [target_idx, okay] = listdlg('ListString', ROI_names, ...
     'SelectionMode', 'single', 'Name', 'Target Selection', ...
-    'PromptString', 'Please select the target ROI. ');
+    'PromptString', 'Please select the target ROI for beamlet calc. ');
 if okay ~= 1
     msgbox('Plan creation aborted');
     return;

WiscPlanPhotonkV125/matlab_frontend/lab/RDXTPS_geometry_setup_wizard_ASP.m

@@ -41,63 +41,63 @@ end
 
 
 %% -----===== BTV and Ring creation =====-----
-waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
-ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
-
-[target_idx okay] = listdlg('ListString', ROI_names, ...
-    'SelectionMode', 'single', 'Name', 'Target Selection', ...
-    'PromptString', 'Please select the target ROI. ');
-if okay ~= 1
-    msgbox('Plan creation aborted');
-    delete(hWaitbar);
-    return;
-end
-
-[BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
-    'BTV margin specification', 1, {'0.6'});
-if isempty(BTV_margin_answer)
-    BTV_margin = 0.6;
-else
-    BTV_margin = str2double(BTV_margin_answer{1});
-end
-
-% target_idx and BTV_margin are set. Expand PTV to BTV
-PTVmask = false(size(Geometry.rhomw));
-% for target_idx = target_indices
-        PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
+% waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
+% ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
+% 
+% [target_idx okay] = listdlg('ListString', ROI_names, ...
+%     'SelectionMode', 'single', 'Name', 'Target Selection', ...
+%     'PromptString', 'Please select the target ROI. ');
+% if okay ~= 1
+%     msgbox('Plan creation aborted');
+%     delete(hWaitbar);
+%     return;
 % end
-if BTV_margin > 0
-    if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
-        bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-        Geometry.BTV = bwD <= BTV_margin;
-    end
-end
-
-% Create btv
-Geometry.ROIS{end+1} = Geometry.ROIS{end};
-Geometry.ROIS{end}.name = 'BTV';
-Geometry.ROIS{end}.num_curves = 0;
-Geometry.ROIS{end}.curves = {};
-Geometry.ROIS{end}.ind = find(Geometry.BTV);
-Geometry.ROIS{end}.visible = false;
-
-% Create ring
-[ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
-    'Ring margin specification', 1, {'1'});
-if isempty(ring_margin_answer)
-    ring_margin = 1;
-else
-    ring_margin = str2double(ring_margin_answer{1});
-end
-bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
-Geometry.Ring = xor(Geometry.Ring, PTVmask);
-Geometry.ROIS{end+1} = Geometry.ROIS{end};
-Geometry.ROIS{end}.name = 'Ring';
-Geometry.ROIS{end}.num_curves = 0;
-Geometry.ROIS{end}.curves = {};
-Geometry.ROIS{end}.ind = find(Geometry.Ring);
-Geometry.ROIS{end}.visible = false;
+% 
+% [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
+%     'BTV margin specification', 1, {'0.6'});
+% if isempty(BTV_margin_answer)
+%     BTV_margin = 0.6;
+% else
+%     BTV_margin = str2double(BTV_margin_answer{1});
+% end
+% 
+% % target_idx and BTV_margin are set. Expand PTV to BTV
+% PTVmask = false(size(Geometry.rhomw));
+% % for target_idx = target_indices
+%         PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
+% % end
+% if BTV_margin > 0
+%     if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
+%         bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+%         Geometry.BTV = bwD <= BTV_margin;
+%     end
+% end
+% 
+% % Create btv
+% Geometry.ROIS{end+1} = Geometry.ROIS{end};
+% Geometry.ROIS{end}.name = 'BTV';
+% Geometry.ROIS{end}.num_curves = 0;
+% Geometry.ROIS{end}.curves = {};
+% Geometry.ROIS{end}.ind = find(Geometry.BTV);
+% Geometry.ROIS{end}.visible = false;
+% 
+% % Create ring
+% [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
+%     'Ring margin specification', 1, {'1'});
+% if isempty(ring_margin_answer)
+%     ring_margin = 1;
+% else
+%     ring_margin = str2double(ring_margin_answer{1});
+% end
+% bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+% Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
+% Geometry.Ring = xor(Geometry.Ring, PTVmask);
+% Geometry.ROIS{end+1} = Geometry.ROIS{end};
+% Geometry.ROIS{end}.name = 'Ring';
+% Geometry.ROIS{end}.num_curves = 0;
+% Geometry.ROIS{end}.curves = {};
+% Geometry.ROIS{end}.ind = find(Geometry.Ring);
+% Geometry.ROIS{end}.visible = false;
 
 % % Create high comformity ring structure for adaptive step size
 % [HCRing_margin_answer] = inputdlg({sprintf('Please enter the HC ring margin (cm):\n(default 1.2 cm or 2 sigma, enter 0 to skip)')}, ...
@@ -115,9 +115,6 @@ Geometry.ROIS{end}.visible = false;
 %     Geometry.HCRing = xor(Geometry.BTV, Geometry.SBTV);
 % end
 
-%% -----===== PTV creation =====-----
-
-
 %% -----===== Save geometry files =====-----
 waitbar(2/total_num_steps, hWaitbar, 'Step 3: Save geometry files');
 
@@ -151,9 +148,9 @@ fid = fopen(fullfile(patient_dir, 'geometry_files', 'Fmw2.bin'), 'w');
 fwrite(fid, Geometry.Fmw2, 'single');
 fclose(fid);
 
-fid = fopen(fullfile(patient_dir, 'geometry_files', 'target_mask.bin'), 'w');
-fwrite(fid, Geometry.BTV, 'single');
-fclose(fid);
+% fid = fopen(fullfile(patient_dir, 'geometry_files', 'target_mask.bin'), 'w');
+% fwrite(fid, Geometry.BTV, 'single');
+% fclose(fid);
 
 % fid = fopen(fullfile(patient_dir, 'geometry_files', 'target_mask_SBTV.bin'), 'w');
 % fwrite(fid, Geometry.SBTV, 'single');

WiscPlanPhotonkV125/matlab_frontend/nrrd2geometry.m

@@ -31,7 +31,7 @@ function [Geometry patient_dir] = nrrd2geometry
     end
     
     % Geometry.rhomw
-    Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
+y    Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
     Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
     [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
     
@@ -78,66 +78,66 @@ function [Geometry patient_dir] = nrrd2geometry
     mkdir(fullfile(patient_dir, 'opt_output'));
 
     %% -----===== BTV and Ring creation =====-----
-    waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
-    ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
-
-    [target_idx okay] = listdlg('ListString', ROI_names, ...
-        'SelectionMode', 'single', 'Name', 'Target Selection', ...
-        'PromptString', 'Please select the target ROI. ');
-    if okay ~= 1
-        msgbox('Plan creation aborted');
-        delete(hWaitbar);
-        return;
-    end
-
-    [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
-        'BTV margin specification', 1, {'0.6'});
-    if isempty(BTV_margin_answer)
-        BTV_margin = 0.6;
-    else
-        BTV_margin = str2double(BTV_margin_answer{1});
-    end
-
-    % target_idx and BTV_margin are set. Expand PTV to BTV
-    PTVmask = false(size(Geometry.rhomw));
-    % for target_idx = target_indices
-            PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
-    % end
-    if BTV_margin > 0
-        if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
-            bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-            Geometry.BTV = bwD <= BTV_margin;
-        end
-    end
-
-    % Create btv
-    Geometry.ROIS{end+1} = Geometry.ROIS{end};
-    Geometry.ROIS{end}.name = 'BTV';
-    Geometry.ROIS{end}.num_curves = 0;
-    Geometry.ROIS{end}.curves = {};
-    Geometry.ROIS{end}.ind = find(Geometry.BTV);
-    Geometry.ROIS{end}.visible = false;
-
-    % Create ring
-    [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
-        'Ring margin specification', 1, {'1'});
-    if isempty(ring_margin_answer)
-        ring_margin = 1;
-    else
-        ring_margin = str2double(ring_margin_answer{1});
-    end
-    bwD = bwdistsc(PTVmask, Geometry.voxel_size);
-    Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
-    Geometry.Ring = xor(Geometry.Ring, PTVmask);
-    Geometry.ROIS{end+1} = Geometry.ROIS{end};
-    Geometry.ROIS{end}.name = 'Ring';
-    Geometry.ROIS{end}.num_curves = 0;
-    Geometry.ROIS{end}.curves = {};
-    Geometry.ROIS{end}.ind = find(Geometry.Ring);
-    Geometry.ROIS{end}.visible = false;
-
-    
-    
+%     waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
+%     ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
+% 
+%     [target_idx okay] = listdlg('ListString', ROI_names, ...
+%         'SelectionMode', 'single', 'Name', 'Target Selection', ...
+%         'PromptString', 'Please select the target ROI. ');
+%     if okay ~= 1
+%         msgbox('Plan creation aborted');
+%         delete(hWaitbar);
+%         return;
+%     end
+% 
+%     [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
+%         'BTV margin specification', 1, {'0.6'});
+%     if isempty(BTV_margin_answer)
+%         BTV_margin = 0.6;
+%     else
+%         BTV_margin = str2double(BTV_margin_answer{1});
+%     end
+% 
+%     % target_idx and BTV_margin are set. Expand PTV to BTV
+%     PTVmask = false(size(Geometry.rhomw));
+%     % for target_idx = target_indices
+%             PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
+%     % end
+%     if BTV_margin > 0
+%         if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
+%             bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+%             Geometry.BTV = bwD <= BTV_margin;
+%         end
+%     end
+% 
+%     % Create btv
+%     Geometry.ROIS{end+1} = Geometry.ROIS{end};
+%     Geometry.ROIS{end}.name = 'BTV';
+%     Geometry.ROIS{end}.num_curves = 0;
+%     Geometry.ROIS{end}.curves = {};
+%     Geometry.ROIS{end}.ind = find(Geometry.BTV);
+%     Geometry.ROIS{end}.visible = false;
+% 
+%     % Create ring
+%     [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
+%         'Ring margin specification', 1, {'1'});
+%     if isempty(ring_margin_answer)
+%         ring_margin = 1;
+%     else
+%         ring_margin = str2double(ring_margin_answer{1});
+%     end
+%     bwD = bwdistsc(PTVmask, Geometry.voxel_size);
+%     Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
+%     Geometry.Ring = xor(Geometry.Ring, PTVmask);
+%     Geometry.ROIS{end+1} = Geometry.ROIS{end};
+%     Geometry.ROIS{end}.name = 'Ring';
+%     Geometry.ROIS{end}.num_curves = 0;
+%     Geometry.ROIS{end}.curves = {};
+%     Geometry.ROIS{end}.ind = find(Geometry.Ring);
+%     Geometry.ROIS{end}.visible = false;
+
+    
+    %% -----===== Save the matlab files =====-----
     waitbar(2.33/total_num_steps, hWaitbar, 'Step 3: Save matlab geometry files');
     % Save matlab geometry file
     save(fullfile(patient_dir, 'matlab_files', 'Geometry.mat'), 'Geometry');
@@ -156,10 +156,9 @@ function [Geometry patient_dir] = nrrd2geometry
     fwrite(fid, Geometry.Fmw2, 'single');
     fclose(fid);
 
-    fid = fopen(fullfile(patient_dir, 'geometry_files', 'target_mask.bin'), 'w');
-    fwrite(fid, Geometry.BTV, 'single');
-    fclose(fid);
-
+%     fid = fopen(fullfile(patient_dir, 'geometry_files', 'target_mask.bin'), 'w');
+%     fwrite(fid, Geometry.BTV, 'single');
+%     fclose(fid);
 
     % hWaitbar = waitbar(1/total_num_steps, 'Step 4, Create optimization geometry');
 

optGoals/get_ROI_goals.m

@@ -5,7 +5,9 @@ function [optGoal, optGoal_beam, optGoal_idx, targetMinMax_idx] = get_ROI_goals(
 % this is a support function for NLP_beamlet_optimizer. It returns the
 % goals planned for specific patients
 
-load(['C:\010-work\003_localGit\WiscPlan_v2\optGoals\' patient, '.mat']);
+% load(['C:\010-work\003_localGit\WiscPlan_v2\optGoals\' patient, '.mat']);
+
+load(['\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed\RODP_files\optGoal.mat']);
 
 optGoal = ROI_goals.optGoal;
 optGoal_beam = ROI_goals.optGoal_beam;

optGoals/make_ROI_goals.m

@@ -384,16 +384,23 @@ end
 function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets)
     optGoal={};
     
+    pathIn = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed\';
+    [minDose, minDose_meta] = nrrdread([pathIn 'RODP_files\FET_FGL005_B1_DP_minDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([pathIn 'RODP_files\FET_FGL005_B1_DP_maxDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
     % -- START DEFINITION OF GOAL --
-    goal_1.name = 'BTV_min';
-    goal_1.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_1.name = 'CTV_min';
+    goal_1.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
     goal_1.ROI_idx = ROI_idx;
     goal_1.imgDim = size(Geometry.data);
-    goal_1.D_final = 60;
+    goal_1.D_final = minDose(ROI_idx);
     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.target_alpha = 1.00;
+    goal_1.target   = 60 % minDose(ROI_idx);
     goal_1.opt_weight = 40 / 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
@@ -401,15 +408,16 @@ function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets)
     % -- END DEFINITION OF GOAL --
     
     % -- START DEFINITION OF GOAL --
-    goal_2.name = 'BTV_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
+    goal_2.name = 'CTV_max';
+    goal_2.ROI_name = Geometry.ROIS{1, 1}.name;
+    ROI_idx = Geometry.ROIS{1, 1}.ind;
     goal_2.ROI_idx = ROI_idx;
     goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 62;
-    goal_2.function = 'max';
+    goal_2.D_final = maxDose(ROI_idx);
+    goal_2.function = 'max_sq';
     goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
+%     goal_1.target_alpha = 1.00;
+    goal_2.target   = 60 % maxDose(ROI_idx);
     goal_2.opt_weight = 2 / 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
@@ -418,11 +426,11 @@ function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets)
     
     % -- START DEFINITION OF GOAL --
     goal_3.name = 'head_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
+    goal_3.ROI_name = Geometry.ROIS{1, 2}.name;
+    ROI_idx = Geometry.ROIS{1, 2}.ind;
     goal_3.ROI_idx = ROI_idx;
     goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 20;
+    goal_3.D_final = 0.1;
     goal_3.function = 'max_sq';
     goal_3.beamlets_pruned = beamlets(ROI_idx, :);
     goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;