Cleanup prep commit

Compare_optPlans.m

@@ -2,11 +2,11 @@
 
 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')
+    resultPath = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1\matlab_files\';
+    
+    load([resultPath 'Geometry.mat'])
+    load([resultPath 'NLP_result.mat'])
+    load([resultPath 'optResults.mat'])
     
     roi_idx = 1;
     nfrac = 1;
@@ -21,4 +21,5 @@ function Compare_optPlans
     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
+    legend('Classical', 'Robust')
 end

NLP_beamlet_optimizer.m

@@ -11,12 +11,9 @@ function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
 % Outputs: full dose image dose: D_full, optimal beamlet weights: w_fin
 %
 % Made by Peter Ferjancic 1. May 2018
-% Last updated: 1. May 2018
+% Last updated: 14. August 2018
 % Inspired by Ana Barrigan's REGGUI optimization procedures
 
-% To-do:
-% - Add robusness aspect (+take worst case scenario, see REGGUI)
-
 N_fcallback1 = 5000;
 N_fcallback2 = 200000;
 patient = 'gbm_005';
@@ -24,19 +21,14 @@ patient = 'gbm_005';
 switch patient
     case 'patient'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
-        blet_in_beam=5; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     case 'tomoPhantom'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
-        blet_in_beam=7; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';        
     case 'phantom_HD'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PD_HD_dicomPhantom';
-        blet_in_beam=20; % this is the number of beamlets in a beam. Called "Mxp" in helicalDosecalcSetup
     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
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
     otherwise
         error('invalid case')
 end
@@ -89,30 +81,10 @@ close(wbar2)
 
 
 %% -- OPTIMIZATION TARGETS --
-switch patient
-    case 'patient'
-%         optGoal = make_ROI_goals(Geometry, beamlets);
-%         optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
-%         N_beamlets_in_beam = 10;
-    case 'tomoPhantom'   
-        optGoal = make_ROI_goals_2(Geometry, beamlets);
-        optGoal_beam = make_ROI_goals_2(Geometry, beamlets_joined);
-        N_beamlets_in_beam = 7;
-    case 'phantom_HD'
-        optGoal = make_ROI_goals(Geometry, beamlets);
-        optGoal_beam = make_ROI_goals(Geometry, beamlets_joined);
-    case 'doggo'
-%         optGoal = make_ROI_goals_DOG(Geometry, beamlets);
-%         optGoal_beam = make_ROI_goals_DOG(Geometry, beamlets_joined);
+make_ROI_goals(Geometry, beamlets, beamlets_joined, patient);
+
+[optGoal, optGoal_beam, optGoal_idx, targetMinMax_idx] = get_ROI_goals(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
 % -- make them robust --
 RO_params=0;
 optGoal_beam = make_robust_optGoal(optGoal_beam, RO_params, beamlets_joined);
@@ -176,26 +148,7 @@ D_full = reshape(beamlets * w_fin, size(Geometry.data));
 %% save outputs
 NLP_result.dose = D_full;
 NLP_result.weights = w_fin;
-% save('C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\matlab_files\NLP_result.mat', 'NLP_result');
-
-switch patient
-    case 'patient'
-        optGoal_idx=[1,2];
-        targetMinMax_idx=[1,3];
-    case 'tomoPhantom'   
-        optGoal_idx=[1,3];
-        targetMinMax_idx=[1,2];
-    case 'phantom_HD'
-
-    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
+save([patient_dir '\matlab_files\NLP_result.mat'], 'NLP_result');
 
 plot_DVH(D_full, optGoal, optGoal_idx, targetMinMax_idx)
 colorwash(Geometry.data, D_full);
@@ -229,6 +182,7 @@ function penalty = eval_f(x, optGoal, nrs_i, sss_i, rrs_i)
     % for each condition
     for goal_i = 1:numel(optGoal)
         switch optGoal{goal_i}.function
+            % min, max, min_sq, max_sq, LeastSquare, min_perc_Volume, max_perc_Volume
             case 'min'
                 % penalize if achieved dose is lower than target dose
                 d_penalty = 1.0e0 * sum(max(0, ...
@@ -320,391 +274,6 @@ function show_joint_beamlets(beamlets, IMGsize, Beam_list)
     
 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.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, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.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 = '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.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
-function optGoal = make_ROI_goals_2(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.imgDim = size(Geometry.data);
-    goal_1.D_final = 62;
-    goal_1.function = 'min';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.00;
-    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, 2}.name;
-    ROI_idx = Geometry.ROIS{1, 2}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 63;
-    goal_2.function = 'max';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.05;
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    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
-    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.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 --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_4.name = 'Ring_max';
-    goal_4.ROI_name = Geometry.ROIS{1, 6}.name;
-    ROI_idx = Geometry.ROIS{1, 6}.ind;
-    goal_4.ROI_idx = ROI_idx;
-    goal_4.imgDim = size(Geometry.data);
-    goal_4.D_final = 56;
-    goal_4.function = 'max_sq';
-    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
-    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_4.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_4;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_DOG(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 = 62;
-    goal_1.function = 'min_sq';
-    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 = 'Doggo_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 20;
-    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_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.2, 0.9, 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 = 'Target_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 63;
-    goal_3.function = 'max_sq';
-    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.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_3.name = 'Doggo_max2';
-    goal_3.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 50;
-    goal_3.function = 'max_sq';
-    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
-function optGoal = make_ROI_goals_DOG_2(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 = 62;
-    goal_1.function = 'min_sq';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.00;
-    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 = 'Doggo_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 20;
-    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_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.2, 0.9, 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 = 'Target_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 63;
-    goal_3.function = 'max_sq';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target_alpha = 1.05;
-    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.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_4.name = 'Doggo_max2';
-    goal_4.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_4.ROI_idx = ROI_idx;
-    goal_4.imgDim = size(Geometry.data);
-    goal_4.D_final = 50;
-    goal_4.function = 'max_sq';
-    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
-    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_4.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_4;
-    % -- END DEFINITION OF GOAL --
-    
-end
-function optGoal = make_ROI_goals_DOG_3(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 = 62;
-    goal_1.function = 'min_perc_Volume';
-    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_1.target_alpha = 1.00;
-    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 = 'Doggo_max';
-    goal_2.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_2.ROI_idx = ROI_idx;
-    goal_2.imgDim = size(Geometry.data);
-    goal_2.D_final = 20;
-    goal_2.function = 'max_sq';
-    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_2.targetImg   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.target   = ones(numel(ROI_idx), 1) * goal_2.D_final;
-    goal_2.opt_weight = 0.2 / numel(ROI_idx); % normalize to volume of target area
-    goal_2.dvh_col = [0.2, 0.9, 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 = 'Target_max';
-    goal_3.ROI_name = Geometry.ROIS{1, 1}.name;
-    ROI_idx = Geometry.ROIS{1, 1}.ind;
-    goal_3.ROI_idx = ROI_idx;
-    goal_3.imgDim = size(Geometry.data);
-    goal_3.D_final = 63;
-    goal_3.function = 'max_perc_Volume';
-    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_3.target_alpha = 1.05;
-    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
-    goal_3.opt_weight = 0.5 / numel(ROI_idx); % normalize to volume of target area
-    goal_3.dvh_col = [0.9, 0.2, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_3;
-    % -- END DEFINITION OF GOAL --
-    
-    % -- START DEFINITION OF GOAL --
-    goal_4.name = 'Doggo_max2';
-    goal_4.ROI_name = Geometry.ROIS{1, 4}.name;
-    ROI_idx = Geometry.ROIS{1, 4}.ind;
-    goal_4.ROI_idx = ROI_idx;
-    goal_4.imgDim = size(Geometry.data);
-    goal_4.D_final = 50;
-    goal_4.function = 'max_sq';
-    goal_4.beamlets_pruned = beamlets(ROI_idx, :);
-    goal_4.target   = ones(numel(ROI_idx), 1) * goal_4.D_final;
-    goal_4.opt_weight = 1 / numel(ROI_idx); % normalize to volume of target area
-    goal_4.dvh_col = [0.2, 0.9, 0.2]; % color of the final DVH plot
-    % assign target
-    optGoal{end+1}=goal_4;
-    % -- 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);
@@ -722,8 +291,8 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     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]};
+    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]};
 
@@ -795,6 +364,3 @@ end
 
 
 
-
-
-

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -1,7 +1,8 @@
 %%
 
 % This version was modified by Peter Ferjancic to accept .nrrd file format
-% for target.
+% for target. It was also rollbacked on some parameters back to the
+% clinical system settings.
 
 % This file has been modified by Surendra Prajapati especially to run
 % WiscPlan for KV beams. Works for running locally as well as in Server.
@@ -18,38 +19,32 @@
 % Surendra edits: num_batches, SAD, pitch, xpmin/max, ypmin/max, Mxp, Nphi 
 % Also edit: 
 % Kernel should always be named Kernels.mat 
-% SAD < 20 was edited from SAD < 50 for clinical system (error message)
-% 
 
 function helicalDosecalcSetup7(patient_dir)
-
-
-% specify where kernel and [geometry] files are located
-% patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\PatientData';
+% -- INPUT:
+% patient_dir: specify where kernel and [geometry] files are located
 num_batches = 4; % number of cores you want to run the beam calculation
 
 iso = [0 0 0]; % Point about which gantry rotations begin
-% SAD = 35;  % source-axis distance for the x-ray source
-SAD = 60;  % source-axis distance for the x-ray source ##
+SAD = 85;  % source-axis distance for the x-ray source ##
 pitch = 0.86; % fraction of beam with couch translates per rotation
 
 % define the overall beam field size for each beam angle
-xpmin = -0.5;     % -field width / 2
-xpmax = 0.5;      % +field width / 2
-ypmin = -0.5;   % -jaw width / 2
-ypmax = 0.5;    % +jaw width / 2
-
+% xpmin = -0.5;     % -field width / 2
+% xpmax = 0.5;      % +field width / 2
+% ypmin = -0.5;   % -jaw width / 2
+% ypmax = 0.5;    % +jaw width / 2
 
 % ##
-xpmin = -5.0;     % -field width / 2
-xpmax = 5.0;      % +field width / 2
-ypmin = -1.0;   % -jaw width / 2
-ypmax = 1.0;    % +jaw width / 2
+xpmin = -20.0;     % -field width / 2
+xpmax = 20.0;      % +field width / 2
+ypmin = -0.3125;  % total jaw width is 0.625 cm
+ypmax = 0.3125;
 
 % y-prime points in the z-direction in the CT coordinate system
 
 % Number of beamlets in the BEV for each direction
-Mxp = 7;  % Mxp = 20;  number of leaves; leaf size is 1/20cm= 0.5mm
+Mxp = 20;  % Mxp = 20;  number of leaves;
 Nyp = 1;  % always 1 for Tomo due to binary mlc
 
 
@@ -95,13 +90,10 @@ Nrot = ceil(abs(zBow - zStern)/(pitch*fieldWidth));
 % Nphi = Nrot*51;  % number of angles used in the calculation
 Nphi = Nrot*21;  % Grozomah
 
-% load(geometry_file);
-
 % define the limits of the angles that will be used for the calculation
 % ##phimin = 0;  % starting angle in radians
 % ##phimax = 2*pi*Nphi;
 
-
 phi = [0:Nphi-1]/Nphi *2*pi*Nrot;
 
 condor_folder = patient_dir;
@@ -151,8 +143,8 @@ if Mxp <= 0 || Nyp <= 0 || Nphi <= 0
     error('Mxp, Nyp, and Nphi must be greater than zero.');
 end
 
-if SAD < 20
-    error('It is recommended that the SAD be greater than 20 cm.');
+if SAD < 50
+    error('It is recommended that the SAD be greater than 50 cm.');
 end
 
 % the xy plane is perpendicular to the isocenter axis of the linac gantry

plot_DVH.m

@@ -0,0 +1,62 @@
+% ---- 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(2)}.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

test_ReadBeamlets.m

@@ -2,6 +2,7 @@
 
 
 patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData';
+patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
 
 
 target_filename = fullfile(patient_dir, 'beamlet_batch_files', 'beamletbatch0.bin');
@@ -14,7 +15,7 @@ beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
 
 % close all
 
-N=95;
+N=59;
 data=beamlet_cell_array{N};
 tabula=zeros(data.x_count, data.y_count, data.z_count);
 tabula(data.non_zero_indices)=data.non_zero_values;