function ROI_goals_prep
% this function creates and saves dose painting plans based on Geometry
% and beamlets created in WiscPlan. Manually change hardcode of whatever
% you want the plan to be!

%% ---=== INPUT PARAMS ===---
patient = 'medivation_01';
outName = 'ROI_goals_simple';

switch patient
    case 'medivation_01'
        paths.in = '\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2\';
        paths.CT_in = ['CT_crop_ 1'];
        paths.target_bin_in  = ['CTV_crop_ 1'];
        paths.target_fzy_in  = ['CTV_crop_ 1'];
        paths.body_bin_in    = ['CTV_crop_ 1'];
        paths.wiscplan = 'F:\21_WiscPlan_data\Medivation_0391811_crop1_2';
        
        
    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_64beam';
    case 'gbm_022'
        patname = 'FET_FGL022';
        paths.in = ['\\Mpufs5\data_wnx1\_Data\Glioma_aus\' patname '\B1\Processed'];
        paths.CT_in = [patname '_B1_CT2FET'];
        paths.target_bin_in  = ['\RODP_files\' patname '_B1_thr2_binPTV'];
        paths.target_fzy_in  = ['\RODP_files\' patname '_B1_thr2_fuzzyCTV'];
        paths.body_bin_in    = ['\RODP_files\' patname '_B1_head_crop_bin'];
        paths.wiscplan = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_022';
    otherwise
        error('invalid case')
end

%% ---=== LOAD DATA ===---

% load Geometry
load([paths.wiscplan '\matlab_files\Geometry.mat']);
fprintf('Loaded geometry ')

% load beamlets
[beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, paths.wiscplan);
fprintf('and beamlets.\n')

%% ---=== GET OPTGOAL ===---
switch patient
    case 'medivation_01'
        ROI_goals.optGoal = make_ROI_goals_medivation_011(Geometry, beamlets);
        ROI_goals.optGoal_beam = make_ROI_goals_medivation_011(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'
        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);
        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_022'
        ROI_goals.optGoal = make_ROI_goals_gbm_022(Geometry, beamlets);
        ROI_goals.optGoal_beam = make_ROI_goals_gbm_022(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
    otherwise
        error('invalid case')
end


%% ---=== SAVE OPTGOAL ===---
fprintf('Writing ROI_goals...')
save([paths.in '\RODP_files\' outName '.mat'], 'ROI_goals')
fprintf(' done!\n')
end

function optGoal = make_ROI_goals_medivation_011(Geometry, beamlets, minDose, maxDose)
    optGoal={};
    
    % -- 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 = 'ring_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 = 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_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={};
    
    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL015\B1\Processed';
    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL015_B1_thr2_DP_minDose.nrrd']);
    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL015_B1_thr2_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';
    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
    goal_1.target   = minDose(ROI_idx); % 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
    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   = maxDose(ROI_idx); % 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
    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) * goal_3.D_final;
    goal_3.opt_weight = 7 / 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_022(Geometry, beamlets, minDose, maxDose)
    optGoal={};
    
    DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL022\B1\Processed';
    [minDose, minDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL022_B1_thr2_DP_minDose.nrrd']);
    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL022_B1_thr2_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';
    goal_1.beamlets_pruned = beamlets(ROI_idx, :);
    goal_1.target   = minDose(ROI_idx); % 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
    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   = maxDose(ROI_idx); % 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
    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) * goal_3.D_final;
    goal_3.opt_weight = 7 / 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 beamlets = get_beamlets(beamlet_cell_array, numVox);
    wbar1 = waitbar(0, 'Creating beamlet array');
    numBeam = size(beamlet_cell_array,2);
    batchSize=100;
    beamlets = sparse(0, 0);
    for beam_i=1:numBeam
        % for each beam define how much dose it delivers on each voxel
        idx=beamlet_cell_array{1, beam_i}.non_zero_indices;

        % break the beamlets into multiple batches
        if rem(beam_i, batchSize)==1;
            beamlet_batch = sparse(numVox, batchSize);
            beam_i_temp=1;
        end

        beamlet_batch(idx, beam_i_temp) = 1000*beamlet_cell_array{1, beam_i}.non_zero_values;
        waitbar(beam_i/numBeam, wbar1, ['Adding beamlet array: #', num2str(beam_i)])

        % add the batch to full set when filled
        if rem(beam_i, batchSize)==0;
            beamlets =[beamlets, beamlet_batch];
        end
        % crop and add the batch to full set when completed
        if beam_i==numBeam;
            beamlet_batch=beamlet_batch(:, 1:beam_i_temp);
            beamlets =[beamlets, beamlet_batch];
        end
        beam_i_temp=beam_i_temp+1;

    end
    close(wbar1)

end