function [beamlets, beamlets_joined, numBeamlet, numBeam, beam_i_list] = get_beam_lets(Geometry, patient_dir)
% this function loads and returns beamlets and joined beams

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

numVox  = numel(Geometry.data);
numBeamlet = size(beamlet_cell_array,2);

if size(Geometry.data, 1)<129
    batchSize = 400;
else
    batchSize = 200;
end

beamlets = get_beamlets(beamlet_cell_array, numVox, batchSize);
% 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_beamlets(beamlet_cell_array, numVox, batchSize);
    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;

        % 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
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