Preparation commit before Geometry overhaul.

Plan to: track paths, add am. loading etc.

ROI_goals_prep.m

@@ -7,17 +7,59 @@ function ROI_goals_prep
 % you want the plan to be!
 
 %% ---=== INPUT PARAMS ===---
-patient = 'gbm_015';
+patient = 'avastin_009';
+outName = 'ROI_goals_DP';
 
 switch patient
+    case 'avastin_009'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
+        paths.CT_in = ['AV009b_1ct_resized'];
+        paths.target_bin_in  = ['\RODP_files\AV009b_flt_seg_thr2_imclose_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\AV009b_flt_seg_thr2_imclose_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\AV009b_seg_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
+    case 'avastin_009_dumb'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
+        paths.CT_in = ['AV009b_1ct_resized'];
+        paths.target_bin_in  = ['\RODP_files\AV009b_flt_seg_thr2_imclose_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\AV009b_flt_seg_thr2_imclose_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\AV009b_seg_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
+    case 'avastin_009_DP'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
+        paths.CT_in = ['AV009b_1ct_resized'];
+        paths.target_bin_in  = ['\RODP_files\AV009b_flt_seg_thr2_imclose_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\AV009b_flt_seg_thr2_imclose_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\AV009b_seg_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
+    case 'gbm_005'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+        paths.CT_in = ['FET_FGL005_B1_CT2FET'];
+        paths.target_bin_in  = ['\RODP_files\FET_FGL005_B1_thr2_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\FET_FGL005_B1_thr2_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\FET_FGL005_B1_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam';
+    case 'gbm_005_dumb'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+        paths.CT_in = ['FET_FGL005_B1_CT2FET'];
+        paths.target_bin_in  = ['\RODP_files\FET_FGL005_B1_thr2_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\FET_FGL005_B1_thr2_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\FET_FGL005_B1_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam_2';
+    case 'gbm_005_DP'
+        paths.in = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+        paths.CT_in = ['FET_FGL005_B1_CT2FET'];
+        paths.target_bin_in  = ['\RODP_files\FET_FGL005_B1_thr2_binPTV'];
+        paths.target_fzy_in  = ['\RODP_files\FET_FGL005_B1_thr2_fuzzyCTV'];
+        paths.body_bin_in    = ['\RODP_files\FET_FGL005_B1_head_crop_bin'];
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam';
     case 'gbm_015'
         paths.in = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
         paths.CT_in = ['FET_FGL015_B1_CT2FET'];
         paths.target_bin_in  = ['\RODP_files\FET_FGL015_B1_thr2_binPTV'];
         paths.target_fzy_in  = ['\RODP_files\FET_FGL015_B1_thr2_fuzzyCTV'];
         paths.body_bin_in    = ['\RODP_files\FET_FGL015_B1_head_crop_bin'];
-        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
-        
+        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015_64beam';
     case 'gbm_022'
         patname = 'FET_FGL022';
         paths.in = ['\\Mpufs5\data_wnx1\_Data\Glioma_aus\' patname '\B1\Processed'];
@@ -42,6 +84,36 @@ fprintf('and beamlets.\n')
 
 %% ---=== GET OPTGOAL ===---
 switch patient
+    case 'avastin_009'
+        ROI_goals.optGoal = make_ROI_goals_avastin_009(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_avastin_009(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+    case 'avastin_009_dumb'
+        ROI_goals.optGoal = make_ROI_goals_avastin_009_dumb(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_avastin_009_dumb(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+    case 'avastin_009_DP'
+        ROI_goals.optGoal = make_ROI_goals_avastin_009(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_avastin_009(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+    case 'gbm_005'
+        ROI_goals.optGoal = make_ROI_goals_gbm_005(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_gbm_005(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+    case 'gbm_005_dumb'
+        ROI_goals.optGoal = make_ROI_goals_gbm_005_dumb(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_gbm_005_dumb(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume
+    case 'gbm_005_DP'
+        ROI_goals.optGoal = make_ROI_goals_gbm_005_DP(Geometry, beamlets);
+        ROI_goals.optGoal_beam = make_ROI_goals_gbm_005_DP(Geometry, beamlets_joined);
+        ROI_goals.optGoal_idx=[1,3]; % indeces of volumes you want on histogram
+        ROI_goals.targetMinMax_idx=[1,2]; % indeces of limits for min/max target volume    
     case 'gbm_015'
         ROI_goals.optGoal = make_ROI_goals_gbm_015(Geometry, beamlets);
         ROI_goals.optGoal_beam = make_ROI_goals_gbm_015(Geometry, beamlets_joined);
@@ -60,11 +132,375 @@ end
 
 %% ---=== SAVE OPTGOAL ===---
 fprintf('Writing ROI_goals...')
-save([paths.in '\RODP_files\ROI_goals.mat'], 'ROI_goals')
+save([paths.in '\RODP_files\' outName '.mat'], 'ROI_goals')
 fprintf(' done!\n')
 end
-
-
+function optGoal = make_ROI_goals_avastin_009_DP(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\AV009b_flt_seg_thr3_imclose_DP_maxDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\AV009b_flt_seg_thr3_imclose_DP_minDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    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 = minDose(ROI_idx);
+    goal_1.function = 'min_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1;
+    goal_1.target   = minDose(ROI_idx); % minDose(ROI_idx);
+    goal_1.opt_weight = 70 / 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 = '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 = maxDose(ROI_idx);
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target_alpha = 1;
+    goal_2.target   = maxDose(ROI_idx); % maxDose(ROI_idx);
+    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 = 'head_max';
+    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.function = 'max_sq';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target_alpha = 1;
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 5 / 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_avastin_009_dumb(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\AV009b_flt_seg_thr3_imclose_DP_maxDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\AV009b_flt_seg_thr3_imclose_DP_minDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    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.function = 'min_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+%     goal_1.target_alpha = 1;
+    goal_1.target   = ones(numel(ROI_idx), 1) * 60; % minDose(ROI_idx);
+    goal_1.opt_weight = 70 / 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 = '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 = 63;
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+%     goal_2.target_alpha = 1;
+    goal_2.target   = ones(numel(ROI_idx), 1) * 63; % maxDose(ROI_idx);
+    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 = 'head_max';
+    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.function = 'max';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target_alpha = 1;
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 5 / 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_avastin_009(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\AV009b_flt_seg_thr3_imclose_DP_maxDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\AV009b_flt_seg_thr3_imclose_DP_minDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    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 = minDose(ROI_idx);
+    goal_1.function = 'min_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1;
+    goal_1.target   = minDose(ROI_idx); % minDose(ROI_idx);
+    goal_1.opt_weight = 70 / 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 = '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 = maxDose(ROI_idx);
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target_alpha = 1;
+    goal_2.target   = maxDose(ROI_idx); % maxDose(ROI_idx);
+    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 = 'head_max';
+    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.function = 'max_sq';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target_alpha = 1;
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 5 / 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_gbm_005(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_seg_thr2.0_DP_minDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_seg_thr2.0_DP_maxDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    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 = minDose(ROI_idx);
+    goal_1.function = 'min_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1;
+    goal_1.target   = minDose(ROI_idx); % minDose(ROI_idx);
+    goal_1.opt_weight = 70 / 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 = '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 = maxDose(ROI_idx);
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target_alpha = 1;
+    goal_2.target   = maxDose(ROI_idx); % maxDose(ROI_idx);
+    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 = 'head_max';
+    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.function = 'max';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target_alpha = 1;
+    goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 5 / 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_gbm_005_dumb(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+%     [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_seg_thr2.0_DP_minDose.nrrd']);
+%     [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_seg_thr2.0_DP_maxDose.nrrd']);
+%     minDose = double(minDose);
+%     maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    goal_1.name = 'PTV_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_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target   = ones(numel(ROI_idx), 1) * 60;
+%     goal_1.target   = minDose(ROI_idx);
+    goal_1.opt_weight = 77 / 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 = 'PTV_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 = 63;
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target   = ones(numel(ROI_idx), 1) * 63;
+%     goal_2.target   = maxDose(ROI_idx);
+    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 = 'head_max';
+    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.function = 'max';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target   = ones(numel(ROI_idx), 1) * 20;
+    goal_3.opt_weight = 5 / 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_gbm_005_DP(Geometry, beamlets, minDose, maxDose)
+    optGoal={};
+    
+    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
+    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_seg_thr2.0_DP_minDose.nrrd']);
+    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_seg_thr2.0_DP_maxDose.nrrd']);
+    minDose = double(minDose);
+    maxDose = double(maxDose);
+    
+    % -- START DEFINITION OF GOAL --
+    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 = minDose(ROI_idx);
+    goal_1.function = 'min_sq';
+    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_1.target_alpha = 1;
+    goal_1.target   = minDose(ROI_idx); % minDose(ROI_idx);
+    goal_1.opt_weight = 70 / 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 = '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 = maxDose(ROI_idx);
+    goal_2.function = 'max_sq';
+    goal_2.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_2.target_alpha = 1;
+    goal_2.target   = maxDose(ROI_idx); % maxDose(ROI_idx);
+    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 = 'head_max';
+    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.function = 'max';
+    goal_3.beamlets_pruned = beamlets(ROI_idx, :);
+    goal_3.target_alpha = 1;
+    goal_2.target   = maxDose(ROI_idx); % maxDose(ROI_idx);
+%     goal_3.target   = ones(numel(ROI_idx), 1) * goal_3.D_final;
+    goal_3.opt_weight = 5 / 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_gbm_015(Geometry, beamlets, minDose, maxDose)
     optGoal={};
     

WiscPlanPhotonkV125/matlab_frontend/NLP_beamlet_optimizer.m

@@ -18,11 +18,19 @@ function [D_full, w_fin, Geometry, optGoal] = NLP_beamlet_optimizer
 
 N_fcallback1 = 5000;
 N_fcallback2 = 200000;
-patient = 'gbm_015';
+patient = 'temp';
+goalsName= 'ROI_goals_DP';
 
 switch patient
+    case 'temp'
+        patient_dir = 'C:\000-temp\DELETEME';
+        DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis'
+    case 'avastin_009'
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009';
+        DP_dir = '\\Mpufs5\data_wnx1\_Data\Avastin\AV009\PF_RODP_analysis';
     case 'gbm_005'
-        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
+%         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';
+        patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam_2';
         DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
     case 'gbm_015'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_015';
@@ -30,6 +38,7 @@ switch patient
     case 'gbm_022'
         patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_022';
         DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL022\B1\Processed';
+    case 'avastin_009'
     otherwise
         error('invalid case')
 end
@@ -46,7 +55,7 @@ load([patient_dir '\matlab_files\Geometry.mat'])
 
 
 %% -- OPTIMIZATION TARGETS --
-load([DP_dir '\RODP_files\ROI_goals.mat']);
+load([DP_dir '\RODP_files\' goalsName '.mat']);
 
 optGoal = ROI_goals.optGoal;
 optGoal_beam = ROI_goals.optGoal_beam;
@@ -219,14 +228,18 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     % Y - Y>0 moves image down
     % Z - in/out.
     
-    shift_mag = 3; % vox of shift
+    shift_mag = 2; % vox of shift
     nrs_scene_list={[0,0,0]};
 
+    
+    % ----====#### CHANGE ROBUSTNESS HERE ####====----
+    
 %     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], [0,0,-1], [0,0,1]};
 %     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]};
 
+    % ----====#### CHANGE ROBUSTNESS HERE ####====----
 
     rrs_scene_list={[0,0,0]};
 
@@ -255,7 +268,7 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
             
             for sss_i = 1 :optGoal{goal_i}.NbrSystSetUpScenarios   % syst. setup scenarios = sss
                 % modify target and beamlets
-                targetImg3=get_RO_sss(targetImg2, sss_scene_list{sss_i});
+                [targetImg3 idxValid]=get_RO_sss(targetImg2, sss_scene_list{sss_i});
                 % beamlets stay the same
                 
                 for rrs_i = 1:optGoal{goal_i}.NbrRangeScenarios   % range scenario = rrs
@@ -268,7 +281,8 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
                     ROI_idx=find(targetImg4>0);
                     
                     if isfield(optGoal{goal_i}, 'target_alpha')
-                        target = double(optGoal{goal_i}.target_alpha * targetIn(ROI_idx));
+                        target = optGoal{goal_i}.target;
+                        target = target(idxValid);
 %                         disp('exists')
                     else
                         target = ones(numel(ROI_idx), 1) .* optGoal{goal_i}.D_final;
@@ -288,9 +302,23 @@ function optGoal = make_robust_optGoal(optGoal, RO_params, beamlets);
     end
 end
 % ------ supp: RO case SSS ------
-function targetImg3=get_RO_sss(targetImg2, sss_scene_shift);
+function [targetImg3 ia]=get_RO_sss(targetImg2, sss_scene_shift);
     % translate the target image 
+    
     targetImg3 = imtranslate(targetImg2,sss_scene_shift);
+    
+    % now we need to figure out if any target voxels fell out during the
+    % shift
+    imgValid = imtranslate(targetImg3,-sss_scene_shift);
+    imgInvalid = (targetImg2-imgValid);
+    
+    idx_1 = find(targetImg2);
+    idx_2 = find(imgInvalid);
+    
+    [idxValid,ia] = setdiff(idx_1,idx_2);
+    
+    [C,ia, ib] = intersect(idx_1,idxValid);
+    
 end
 
 

WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7.m

@@ -23,7 +23,7 @@
 function helicalDosecalcSetup7(patient_dir)
 % -- INPUT:
 % patient_dir: specify where kernel and [geometry] files are located
-num_batches = 4; % number of cores you want to run the beam calculation
+num_batches = 2; % number of cores you want to run the beam calculation
 
 iso = [0 0 0]; % Point about which gantry rotations begin
 SAD = 85;  % source-axis distance for the x-ray source ##
@@ -44,7 +44,7 @@ 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 = 20;  % Mxp = 64;  number of leaves;
+Mxp = 64;  % Mxp = 64;  number of MLC leaves;
 Nyp = 1;  % always 1 for Tomo due to binary mlc
 
 

WiscPlanPhotonkV125/matlab_frontend/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'], 'Select CT image');
+    [IMG_in, IMG_path, filterIdx] = uigetfile([{'*.nrrd'; '*.am'}], 'Select CT image', 'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\');
     
 %     IMG_in = 'test_ct.nrrd';
 %     IMG_path = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\';
@@ -24,46 +24,59 @@ function [Geometry patient_dir] = nrrd2geometry
     % +Geometry.BTV
     % +Geometry.ring
     
-    if IMG_in==0 % no file selected, no change to CT data
-        warning('No file selected!')
-    else 
-        [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
+    switch filterIdx
+        case 0
+            warning('No file selected!')
+        case 1
+            [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
+%             warning('NRRD file loaded!')
+            Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
+            Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
+            [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
+            Geometry.voxel_size = meta.spacedirections;
+            Geometry.start=meta.spaceorigin;
+        case 2
+            imgIn=am2mat([IMG_path, IMG_in]);
+            Geometry.data = permute(imgIn.data, [2,1,3]);
+%             warning('NRRD file loaded!')
+            Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
+            Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
+            [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
+            
+            Geometry.voxel_size = imgIn.voxel_size;
+            Geometry.start=imgIn.start;
+        case 3
+            % something selected
+        otherwise
+%             error('You should never see this.')
     end
     
-    % Geometry.rhomw
-    Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
-    Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
-    [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
-    
-    % Geometry.voxel_size
-%     VoxSizeString = meta.spacedirections;
-%     C=textscan(VoxSizeString(2:end-1), '%f', 9,'Delimiter',{',',') ('});
-%     C= C{1};
-%     Geometry.voxel_size = [C(1), C(5), C(9)];
-    Geometry.voxel_size = meta.spacedirections;
-    
-    % Geometry.start
-%     StartString = meta.spaceorigin;
-%     C=textscan(StartString(2:end-1), '%f', 3,'Delimiter',',');
-%     Geometry.start=C{1};
-%     Geometry.start=Geometry.start';
-    Geometry.start=meta.spaceorigin;
     
     % Geometry.ROIS
     
     %% -----===== Get mask/contour files =====-----
-    [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\';
-    
-    if TRGT_in==0 % no file selected, no change to CT data
-        warning('No file selected!')
-    else 
-        [TRGT_img, meta]=nrrdread([TRGT_path, TRGT_in]);
+    [TRGT_in, TRGT_path, filterIdx2] = uigetfile([{'*.nrrd'; '*.am'}], 'Select file with Target', 'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\');
+    
+    switch filterIdx2
+        case 0
+            warning('No file selected!')
+        case 1
+            [TRGT_img, meta]=nrrdread([TRGT_path, TRGT_in]);
+            warning('NRRD file loaded!')
+            Geometry.ROIS{1}.ind = find(TRGT_img>0);
+            Geometry.ROIS{1}.name= 'Target';
+        case 2
+            TRGT_in=am2mat([TRGT_path, TRGT_in]);
+            warning('AM file loaded!')
+            TRGT_img = permute(TRGT_in.data, [2,1,3]);
+            Geometry.ROIS{1}.ind = find(TRGT_img>0);
+            Geometry.ROIS{1}.name= 'Target';
+        case 3
+            % something selected
+        otherwise
+%             error('You should never see this.')
     end
-    Geometry.ROIS{1}.ind = find(TRGT_img>0);
-    Geometry.ROIS{1}.name= 'Target';
-
+    
     
         %% -----===== Save geometry files =====-----
     patient_dir = uifile('getdir', 'Save the patient data to directory');

planEvaluation/Compare_optPlans.m

@@ -2,24 +2,57 @@
 
 function Compare_optPlans
     
-    resultPath = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1\matlab_files\';
-    
+    resultPath = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_Avastin009\matlab_files\';
     load([resultPath 'Geometry.mat'])
-    load([resultPath 'NLP_result.mat'])
-    load([resultPath 'optResults.mat'])
     
     roi_idx = 1;
     nfrac = 1;
+    
+    load([resultPath 'optResults_backup_Dumb.mat'])
     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);
+    [dvh_TP dosebins_TP] = dvhist(dose_WP,Geometry,roi_idx,nfrac);
+    [dvh_TP_head dosebins_TP_head] = dvhist(dose_WP,Geometry,2,nfrac);
+    
+    load([resultPath 'NLP_result_backup_Dumb.mat'])
+    dose_NLP_1 = NLP_result.dose;
+    [dvh_NLP_1 dosebins_NLP_1] = dvhist(dose_NLP_1,Geometry,roi_idx,nfrac);
+    [dvh_NLP_head_1 dosebins_NLP_head_1] = dvhist(dose_NLP_1,Geometry,2,nfrac);
+    
+    load([resultPath 'NLP_result_backup_prob.mat'])
+%     load([resultPath 'NLP_result.mat'])
+    dose_NLP_2 = NLP_result.dose;
+    [dvh_NLP_2 dosebins_NLP_2] = dvhist(dose_NLP_2,Geometry,roi_idx,nfrac);
+    [dvh_NLP_head_2 dosebins_NLP_head_2] = dvhist(dose_NLP_2,Geometry,2,nfrac);
+    
+    load([resultPath 'NLP_result_backup_DP.mat'])
+    dose_NLP_3 = NLP_result.dose;
+    [dvh_NLP_3 dosebins_NLP_3] = dvhist(dose_NLP_3,Geometry,roi_idx,nfrac);
+    [dvh_NLP_head_3 dosebins_NLP_head_3] = dvhist(dose_NLP_3,Geometry,2,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);
+%     plot(dosebins_TP, dvh_TP,'Color', [0.9,0.2,0.2],'LineStyle', '--','DisplayName', Geometry.ROIS{roi_idx}.name);
+%     plot(dosebins_NLP_1, dvh_NLP_1,'Color', [0.9,0.2,0.2],'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+%     plot(dosebins_NLP_2, dvh_NLP_2,'Color', [0.9,0.2,0.2],'LineStyle', '-.','DisplayName', Geometry.ROIS{roi_idx}.name);
+    plot(dosebins_NLP_3, dvh_NLP_3,'Color', [0.9,0.2,0.2],'LineStyle', '--','DisplayName', Geometry.ROIS{roi_idx}.name);
+    
+%     plot(dosebins_TP_head, dvh_TP_head,'Color', [0.2,0.9,0.2],'LineStyle', '--','DisplayName', Geometry.ROIS{roi_idx}.name);
+%     plot(dosebins_NLP_head_1, dvh_NLP_head_1,'Color', [0.2,0.9,0.2],'LineStyle', '-','DisplayName', Geometry.ROIS{roi_idx}.name);
+%     plot(dosebins_NLP_head_2, dvh_NLP_head_2,'Color', [0.5,0.9,0.2],'LineStyle', '-.','DisplayName', Geometry.ROIS{roi_idx}.name);
+    plot(dosebins_NLP_head_3, dvh_NLP_head_3,'Color', [0.5,0.9,0.2],'LineStyle', '--','DisplayName', Geometry.ROIS{roi_idx}.name);
+    
     hold off
-    legend('Classical', 'Robust')
+%     legend('DP target', 'DP head')
+%     legend('Robust target', 'Robust head')
+%     legend('Prob. target', 'DP target', 'Prob. head', 'DP head')
+%     legend('TP target', 'Robust target', 'TP head', 'Robust head')
+%     legend('Robust target', 'Prob. target', 'Robust head', 'Prob. head')
+    legend('Prob. target', 'DP target', 'Prob. head', 'DP head')
+
+    axis([0 100 0 100])
+    xlabel('Dose', 'Color','w')
+    ylabel('Volume', 'Color','w')
 end
+
+% set(gca,'Color','k') % plot background black
+% whitebg('k')

planEvaluation/compare_RO_norm_DVH_colorwash.m

@@ -7,8 +7,8 @@ function compare_RO_norm_DVH_colorwash(D_full)
 % path_in='C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\';
 % path_in='C:\010-work\003_localGit\WiscPlan_v2\data\Tomo_Phantom\';
 % path_in='C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_dog5_3\';
-path_in='C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
-
+% path_in='C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
+path_in = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005_16beam';
 
 
 
@@ -16,9 +16,10 @@ path_in='C:\010-work\003_localGit\WiscPlan_v2\data\PatientData_ozzy1';
 load([path_in '\matlab_files\optResults.mat'])
 D_norm = optResults.dose{end};
 
-orthoslice(D_full-D_norm, [-15, 15])
+orthoslice(D_full-D_norm, [-30, 30])
 
 % orthoslice(D_full, [0,70])
-% colorwash(Geometry.data, D_full, [500, 1500], [0,70])
+% colorwash(Geometry.data, D_full, [500, 1500], [0,80])
+% colorwash(Geometry.data, D_norm, [500, 1500], [0,80])
 
 end

planEvaluation/plot_DVH.m

@@ -20,7 +20,7 @@ function plot_DVH(dose, optGoal, optGoal_idx, targetMinMax_idx)
         / num_vox;
     title([num2str(perc_vox_min) ' % below 0.95D_0, ' num2str(perc_vox_max) ' % above D 1.07D_0'])
     
-    
+    axis([0 100 0 100])
 end
 function plot_DVH_robust(dose, optGoal, optGoal_idx)
     % this function plots the DVHs of the given dose

vol_prep.m

@@ -23,7 +23,7 @@ paths.in = ['\\Mpufs5\data_wnx1\_Data\Glioma_aus\' patname '\' timepoint '\Proce
 paths.CT_in = [patname '_' timepoint '_CT2FET'];
 paths.target_bin_in  = [ patname '_' timepoint '_seg_thr2.0'];
 paths.target_fzy_in  = [ patname '_' timepoint '_seg_combined'];
-paths.body_bin_in    = [ patname '_' timepoint '_head_crop_bin'];
+paths.body_bin_in    = [ patname '_' timepoint '_head'];
 
 % paths.target_bin_in  = [patname '_' timepoint '_thr2'];
 % paths.target_fzy_in  = [patname '_' timepoint '_thr2'];
@@ -31,11 +31,11 @@ paths.body_bin_in    = [ patname '_' timepoint '_head_crop_bin'];
 
 
 margins.CTV = 1.0; % GTV->CTV margin, in cm
-margins.PTV = 1.0; % CTV->PTV margin, in cm
+margins.PTV = 0.5; % CTV->PTV margin, in cm
 margins.body = 7.0; % how far from tumor do we still include body, in cm
 
-dose.min = 80;
-dose.max = 83;
+dose.min = 60;
+dose.max = 63;
 
 
 %% ---=== FUNCTION CALLS ===---
@@ -86,6 +86,12 @@ body(PTV>0) = 0;
 body(bwD>margins.body) = 0; % how far from the tumor are we still interested in the body
 
 %% save outputs
+filename = [paths.in '\RODP_files\' paths.target_bin_in '_binCTV.nrrd'];
+matrix = single(CTV);
+pixelspacing = GTV_meta.spacedirections;
+origin = GTV_meta.spaceorigin;
+nrrdWriter(filename, matrix, pixelspacing, origin, 'raw');
+
 filename = [paths.in '\RODP_files\' paths.target_bin_in '_binPTV.nrrd'];
 matrix = single(PTV);
 pixelspacing = GTV_meta.spacedirections;
@@ -256,9 +262,11 @@ end
 function dose_min = f_dose_min(CTV_fuzzy, D_min)
 % get low boundary for ideal dose plan
 dose_min = D_min* min(1, 2* max(0, (CTV_fuzzy - 0.25)));
+dose_min(CTV_fuzzy>0.98) = 1.2*D_min;
 % orthoslice(dose_min)
 end
 function dose_max = f_dose_max(CTV_fuzzy, D_max)
 % get hihg boundary for ideal dose plan
 dose_max = D_max* min(1, 2* max(0, (CTV_fuzzy - 0.0)));
+dose_max(CTV_fuzzy>0.98) = 1.2*D_max;
 end