function [ results ] = read_ryan_beamlets ( filename, mode, weights )
%READ_RYAN_BEAMLETS Read beamlets in Ryan's format into cell array
%
% Example:
%   beamlet_cell_array = read_ryan_beamlets (beamlet_batch_filename)
%   beamlet_cell_array = read_ryan_beamlets (beamlet_batch_filename, 'ryan')
%   full_dose_grid = read_ryan_beamlets (beamlet_batch_filename, 'ryan sum', weights)
%
% See also open_beamlet_batch.
%
% Copyleft 1999-, XMO.

if nargin < 2
    mode = 'Ryan';
end

fid = fopen(filename,'r');

if ( fid == -1 )
    error(['cannot open "' filename '"']);
end

Nbeamlets = fread(fid,1,'int');
if Nbeamlets < 0 || Nbeamlets > 327680
    error([num2str(Nbeamlets) ' beamlets designated']);
elseif Nbeamlets == 0
    if strcmpi(mode, 'ryan')
        results = cell(0);
        return;
    elseif strcmpi(mode, 'ryan sum')
        results = 0;
        return;
    end
end

switch lower(mode)
    case 'info'
        results = Nbeamlets;
        % This returns an empty cell array with all meta info
        fprintf('Number of beamlets = %8d\n', Nbeamlets);
        beamlets = cell(1,Nbeamlets);
        for k=1:Nbeamlets
            beamlet.num = fread(fid,1,'int');
            beamlet.x_count = fread(fid,1,'int');
            beamlet.y_count = fread(fid,1,'int');
            beamlet.z_count = fread(fid,1,'int');
            Nind = fread(fid,1,'int');  % number of non-zero indices & values
            % skip beamlet data (one int32 and one single array)
            fseek(fid, Nind * 2 * 4, 'cof');
            
            beamlets{k} = beamlet;
        end
    case {'ryan', 'ryan sum'}
        beamlets = cell(1,Nbeamlets);
        for k=1:Nbeamlets
            beamlet.num = fread(fid,1,'int');
            beamlet.x_count = fread(fid,1,'int');
            beamlet.y_count = fread(fid,1,'int');
            beamlet.z_count = fread(fid,1,'int');
            Nind = fread(fid,1,'int');  % number of non-zero values
            beamlet.non_zero_indices = int32(fread(fid,Nind,'int'));
            beamlet.non_zero_values = single(fread(fid,Nind,'float'));
            % add unity to the indices to convert them from C to Matlab style indexing:
            beamlet.non_zero_indices = beamlet.non_zero_indices + 1;
            beamlet.num = beamlet.num + 1;
            
            if any(beamlet.non_zero_indices < 0)
                error('index negative');
            end
            if any(beamlet.non_zero_indices > beamlet.x_count*beamlet.y_count*beamlet.z_count)
                error('index out of range');
            end
            
            if strcmpi(mode, 'ryan') % save all beamlets in cell array
                beamlets{k} = beamlet;
            elseif strcmpi(mode, 'ryan sum')
                if k == 1 % allocate dose matrix in the first loop
                    dose = zeros(beamlet.x_count, beamlet.y_count, beamlet.z_count);
                end
                dose(beamlet.non_zero_indices) = dose(beamlet.non_zero_indices) + ...
                    beamlet.non_zero_values * weights(beamlet.num);
%                 dose(beamlet.non_zero_indices) = dose(beamlet.non_zero_indices) + ...
%                     beamlet.non_zero_values * weights(k);
            end
            
            results = beamlets;
        end
        if strcmpi(mode, 'ryan sum')
            % rename for return
            % TODO ambiguous output value
            if Nbeamlets > 0
                results = dose;
            else
                results = 0;
            end
        end
    case 'peter'
        beamlets = cell(1,Nbeamlets);
        ROI_idxList = weights;
        wbar1 = waitbar(0, 'Creating beamlet array');
        for k=1:Nbeamlets
           
            waitbar(k/Nbeamlets, wbar1, ['Extracting beamlet array: #', num2str(k)])
            beamlet.num = fread(fid,1,'int');
            beamlet.x_count = fread(fid,1,'int');
            beamlet.y_count = fread(fid,1,'int');
            beamlet.z_count = fread(fid,1,'int');
            Nind = fread(fid,1,'int');  % number of non-zero values
            non_zero_indices = int32(fread(fid,Nind,'int'));
            non_zero_values = single(fread(fid,Nind,'float'));
            % add unity to the indices to convert them from C to Matlab style indexing:
            non_zero_indices = non_zero_indices + 1;
            beamlet.num = beamlet.num + 1;
            if any(non_zero_indices < 0)
                error('index negative');
            end
            if any(non_zero_indices > beamlet.x_count*beamlet.y_count*beamlet.z_count)
                error('index out of range');
            end
            
            % load only the voxels within provided ROI
            [ROI_idxIntersect, ia, ~] = intersect (non_zero_indices, ROI_idxList);
            beamlet.non_zero_indices = ROI_idxIntersect;
            beamlet.non_zero_values = non_zero_values(ia);
                
            beamlets{k} = beamlet;
            
        end
        results = beamlets;
        close(wbar1)
        
    case 'sparse'
        nzmax = 0;
        nBmax = 500;
        for k=1:nBmax%Nbeamlets
            % skip 4 integers
            fseek(fid, 4 * 4, 'cof');
            % get number of non-zero values
            Nind = fread(fid,1,'int');
            nzmax = nzmax + Nind;
            % skip beamlet data (one int32 and one single array)
            fseek(fid, Nind * 2 * 4, 'cof');
        end
        
        % close and open again
        % TODO this doesn't fit the normal workflow
        fclose(fid);
        fid = fopen(filename,'r');
        if ( fid == -1 )
            error(['cannot open "' filename '"']);
        end
        Nbeamlets = fread(fid,1,'int');
        
        % Return is num_voxels x num_beamlets matrix
        % i === row === 1 to num_voxels
        % j === col === 1 to num_beamlets
        S = sparse([], [], [], 256*256*158, nBmax, nzmax);
        
        idx = 1;
        xpg = xprogressbar();
        for k=1:nBmax%Nbeamlets
            xpg.update(k/Nbeamlets);
            beamlet.num = fread(fid,1,'int')
            beamlet.x_count = fread(fid,1,'int');
            beamlet.y_count = fread(fid,1,'int');
            beamlet.z_count = fread(fid,1,'int');
            Nind = fread(fid,1,'int');  % number of non-zero values
            beamlet.non_zero_indices = int32(fread(fid,Nind,'int'));
            beamlet.non_zero_values = single(fread(fid,Nind,'float'));
            % add unity to the indices to convert them from C to Matlab style indexing:
            beamlet.non_zero_indices = beamlet.non_zero_indices + 1;
            
            if any(beamlet.non_zero_indices < 0)
                error('index negative');
            end
            if any(beamlet.non_zero_indices > beamlet.x_count*beamlet.y_count*beamlet.z_count)
                error('index out of range');
            end
            
            % TODO CRITICAL This will break if beamlet file is not in sequence
            S(beamlet.non_zero_indices, beamlet.num+1) = beamlet.non_zero_values;
%             iVector(idx:idx+Nind-1) = beamlet.non_zero_indices;
%             jVector(idx:idx+Nind-1) = (beamlet.num+1) * ones(Nind, 1);
%             sVector(idx:idx+Nind-1) = beamlet.non_zero_values;
%             idx = idx + Nind;
        end
        results = S;
    otherwise
        error('unknown mode');
end

fclose(fid);

end