function [beamlets, numBeamlet] = get_beamlets(Geometry, patient_dir, OptGoals)
% 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(size(Geometry.data));
for i = 1:numel(OptGoals.data)
    tabula(OptGoals.data{i}.ROI_idx) = 1;
end
ROI_idxList = find(tabula>0);

for i_sss = OptGoals.sss_scene_list
    tabula2 = imtranslate(tabula, i_sss{1});
    ROI_idxList2 = find(tabula2>0);
    ROI_idxList = [ROI_idxList; ROI_idxList2];
end
ROI_idxList = unique(ROI_idxList);


%% 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  = numel(Geometry.data);
numBeamlet = size(beamlet_cell_array,2);
batchSize = 200;
beamlets = get_beamlets2(beamlet_cell_array, numVox, batchSize, ROI_idxList);
% beamlets = beamlet_cell_array;

% --- join beamlets into beams
% load([patient_dir '\all_beams.mat'])
% beamletOrigin=[0 0 0];
% beam_i=0;
% beam_i_list=[];
% for beamlet_i = 1:numel(all_beams)
%     newBeamletOrigin = all_beams{beamlet_i}.ip;
%     if any(newBeamletOrigin ~= beamletOrigin)
%         beam_i = beam_i+1;
%         beamletOrigin = newBeamletOrigin;
%     end
%     beam_i_list=[beam_i_list, beam_i];
% end
% 
% 
% wbar2 = waitbar(0, 'merging beamlets into beams');
% numBeam=numel(unique(beam_i_list));
% beamlets_joined=sparse(size(beamlets,1), numBeam);
% for beam_i = 1:numel(unique(beam_i_list))
%     beam_full = sum(beamlets(:,beam_i_list == beam_i), 2); 
%     beamlets_joined(:,beam_i) = beam_full;
%     waitbar(beam_i/numBeam, wbar2)
% end
% 
% close(wbar2)
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