Commit before script overhaul.

NLP_beamlet_optimizer.m

@@ -1,4 +1,10 @@
 
+
+% [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer;
+% orthoslice(D_full, [0,70])
+% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
+
+
 function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
 % This function performs the beamlet optimization
 % Inputs: none. Still needs to find "Geometry" and "beamlets" from Wiscplan
@@ -11,9 +17,9 @@ function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
 % To-do:
 % - Add robusness aspect (+take worst case scenario, see REGGUI)
 
-N_fcallback1 = 1000;
-N_fcallback2 = 150000;
-patient = 'doggo';
+N_fcallback1 = 5000;
+N_fcallback2 = 200000;
+patient = 'gbm_005';
 
 switch patient
     case 'patient'
@@ -28,6 +34,9 @@ switch patient
     case 'doggo'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_dog5_3';
         blet_in_beam=5; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
+    case 'gbm_005'
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\gbm_aus_Data';
+        blet_in_beam=7; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     otherwise
         error('invalid case')
 end
@@ -98,6 +107,9 @@ switch patient
 
         optGoal = make_ROI_goals_DOG_2(Geometry, beamlets);
         optGoal_beam = make_ROI_goals_DOG_2(Geometry, beamlets_joined);
+    case 'gbm_005'
+        optGoal = make_ROI_goals_gbm_005(Geometry, beamlets);
+        optGoal_beam = make_ROI_goals_gbm_005(Geometry, beamlets_joined);
     otherwise
         error('invalid case')
 end
@@ -178,17 +190,18 @@ switch patient
     case 'doggo'
         optGoal_idx=[1,2];
         targetMinMax_idx=[1,3];
+    case 'gbm_005'
+        optGoal_idx=[1,3];
+        targetMinMax_idx=[1,2];
     otherwise
         error('invalid case')
 end
 
 plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
+colorwash(Geometry.data, D_full);
 % plot_DVH_robust(D_full, optGoal, optGoal_idx)
 end
 
-% orthoslice(D_full, [0,70])
-% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
-
 %% support functions
 % ---- PENALTY FUNCTION ----
 function penalty = get_penalty(x, optGoal)
@@ -639,6 +652,59 @@ function optGoal = make_ROI_goals_DOG_3(Geometry, beamlets)
     % -- END DEFINITION OF GOAL --
     
 end
+function optGoal = make_ROI_goals_gbm_005(Geometry, beamlets)
+    optGoal={};
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'Target_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.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 = 10 / 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, 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 = 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 = 0.8 / 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 = 'RING 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.imgDim = size(Geometry.data);
+    goal_3.D_final = 10;
+    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 = 2 / 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
 
 % ---- MAKE ROI ROBUST ----
 function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
@@ -652,15 +718,19 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % X - X>0 moves image right
     % Y - Y>0 moves image down
     % Z - in/out.
+    
+    shift_mag = 1; % vox of shift
     nrs_scene_list={[0,0,0]};
+
     sss_scene_list={[0,0,0]};
+%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0]};
+%     sss_scene_list={[0,0,0], [-shift_mag,0,0], [shift_mag,0,0], [0,-shift_mag,0], [0,shift_mag,0],...
+%         [-shift_mag*2,0,0], [shift_mag*2,0,0], [0,-shift_mag*2,0], [0,shift_mag*2,0]};
+
+
     rrs_scene_list={[0,0,0]};
-    
-%     nrs_scene_list={[0,0,0]};
-%     sss_scene_list={[0,0,0], [-2,0,0], [2,0,0], [0,-2,0], [0,2,0]};
-%     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\CDP_data\CDP5_DP_target.nrrd');
 %     [targetIn, meta] = nrrdread('C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom\Tomo_DP_target.nrrd');
     
     for i = 1:numel(optGoal)
@@ -722,68 +792,7 @@ function targetImg3=get_RO_sss(targetImg2, sss_scene_shift);
     targetImg3 = imtranslate(targetImg2,sss_scene_shift);
 end
 
-% ---- DVH PLOT FUNCTION ----
-function plot_DVH(dose, optGoal, optGoal_idx, targetMinMax_idx)
-    % this function plots the DVHs of the given dose
-    nfrac=1;
-    
-    figure
-    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
-    % get % of volume above/below threshold
-    num_vox = numel(optGoal{targetMinMax_idx(1)}.target);
-    perc_vox_min = 100* numel(find(dose(optGoal{targetMinMax_idx(1)}.ROI_idx) - optGoal{targetMinMax_idx(1)}.target*0.95 <0)) ...
-        / num_vox;
-    perc_vox_max = 100* numel(find(dose(optGoal{targetMinMax_idx(2)}.ROI_idx) - optGoal{targetMinMax_idx(1)}.target*1.07 >0)) ...
-        / num_vox;
-    title([num2str(perc_vox_min) ' % below 0.95D_0, ' num2str(perc_vox_max) ' % above D 1.07D_0'])
-    
-    
-end
-function plot_DVH_robust(dose, optGoal, optGoal_idx)
-    % this function plots the DVHs of the given dose
-    nfrac=1;
-    
-    figure
-    hold on
-    for goal_i=optGoal_idx
-        % for all the perturbations
-        for nrs_i = 1:optGoal{1}.NbrRandScenarios 
-            for sss_i = 1:optGoal{1}.NbrSystSetUpScenarios % syst. setup scenarios = ss
-                for rrs_i = 1:optGoal{1}.NbrRangeScenarios % range scenario = rs
-                    % plot the histogram
-                    ROI_idx = optGoal{goal_i}.nrs{nrs_i}.sss{sss_i}.rrs{rrs_i}.ROI_idx;
-                    
-                    [dvh dosebins] = get_dvh_data(dose,ROI_idx,nfrac);
-                    plot(dosebins, dvh,'Color', optGoal{goal_i}.dvh_col,'LineStyle', '-','DisplayName', optGoal{goal_i}.ROI_name);
-                end
-            end
-        end
-    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
+
 
 
 

nrrd2geometry.m

@@ -3,7 +3,7 @@
 function [Geometry patient_dir] = nrrd2geometry
     
     %% -----===== Get CT files =====-----
-    [IMG_in, IMG_path, ~] = uigetfile(['C:\010-work\003_localGit\WiscPlan_v2\data\', '*.nrrd']);
+    [IMG_in, IMG_path, ~] = uigetfile(['C:\010-work\003_localGit\WiscPlan_v2\data\', '*.nrrd'], 'Select CT image');
     
 %     IMG_in = 'test_ct.nrrd';
 %     IMG_path = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\';
@@ -52,7 +52,7 @@ function [Geometry patient_dir] = nrrd2geometry
     % Geometry.ROIS
     
     %% -----===== Get mask/contour files =====-----
-    [TRGT_in, TRGT_path, ~] = uigetfile(['C:\010-work\003_localGit\WiscPlan_v2\CDP_data\','*.nrrd']);
+    [TRGT_in, TRGT_path, ~] = uigetfile(['C:\010-work\003_localGit\WiscPlan_v2\CDP_data\','*.nrrd'], 'Select file with Target');
 %     TRGT_in = 'test_seg.nrrd';
 %     TRGT_path = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\';
     
@@ -66,8 +66,8 @@ function [Geometry patient_dir] = nrrd2geometry
 
     
         %% -----===== Save geometry files =====-----
-%     patient_dir = uifile('getdir', 'Save the patient data to directory');
-    patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\';
+    patient_dir = uifile('getdir', 'Save the patient data to directory');
+%     patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\';
     
     % make directories
     mkdir(patient_dir);