function [beamlets, numBeamlet] = get_beamlets_fullRO(imgDim, patient_dir, optGoal)
% this function loads and returns beamlets and joined beams

% ROI_idxList = [];
% for i = 1:numel(OptGoals.data)
%     ROI_idxList = [ROI_idxList; OptGoals.data{i}.ROI_idx];
% end
% ROI_idxList = unique(ROI_idxList);

tabula = zeros(imgDim);
for i = 1:numel(optGoal)
    tabula(optGoal{i}.ROI_idx) = 1;
end
ROI_idxList = find(tabula>0);


%% add idx lists for each goal!

beamlet_batch_filename = [patient_dir '\' 'batch_dose.bin'];
% beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'ryan');
beamlet_cell_array = read_ryan_beamlets(beamlet_batch_filename, 'peter', ROI_idxList);

numVox  = prod(imgDim);
numBeamlet = size(beamlet_cell_array,2);
batchSize = 200;
beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize, ROI_idxList);
% beamlets = beamlet_cell_array;

end

% ---- GET BEAMLETS ----
function beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize, ROI_idxList);
    wbar1 = waitbar(0, 'Creating beamlet array');
    numBeam = size(beamlet_cell_array,2);
    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;
        [ROI_idxIntersect, ia, ~] = intersect (idx, ROI_idxList);
%         if ~isequal(ROI_idxIntersect,idx)
%             warning('this')
%         end
        
%         ROI_idxIntersect = beamlet_cell_array{1, beam_i}.non_zero_indices;
        
        if isempty(ROI_idxIntersect)
            warning(['Beamlet ' num2str(beam_i) ' is empty!'])
        end
        
        % 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;
        beamlet_batch(ROI_idxIntersect, beam_i_temp) = beamlet_cell_array{1, beam_i}.non_zero_values(ia);
        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
function show_joint_beamlets(beamlets, IMGsize, Beam_list)
    % this function overlays and plots multiple beamlets. The goal is to
    % check whether some beamlets are part of the same beam manually.
    
    beam=sum(beamlets(:, Beam_list), 2);
    
    beamImg=reshape(full(beam), IMGsize);
        
    orthoslice(beamImg)
    
end