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!


%% ---=== LOAD DATA ===---
patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatData_ausGli_005';

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

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

%% ---=== GET OPTGOAL ===---
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

%% ---=== SAVE OPTGOAL ===---
fprintf('Writing ROI_goals...')
DP_dir = '\\Mpufs5\data_wnx1\_Data\Glioma_aus\FET_FGL005\B1\Processed';
save([DP_dir '\RODP_files\ROI_goals.mat'], 'ROI_goals')
fprintf(' done!\n')
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_DP_minDose.nrrd']);
    [maxDose, maxDose_meta] = nrrdread([DP_dir '\RODP_files\FET_FGL005_B1_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_sq';
    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