function data = amiraReader(varargin)

%  this m-file skips the header of amira files to give you the straight data
%
%  USEAGE:  either type in what you want (valuable for scripting or calling from another program) 
%           > data = amiraReader('C:\Data\Human\Test_foo.am');
%           or
%           > data = amiraReader; % leave the argument blank and use the file selector
%
%   Author: DLB - 3/20/09, help from RTF's geom2am.m @ line 38.
%
%%

warning off all;  % you'll get truncation warnings because you read inf in line 24.

switch nargin 
    case 1
        filename = varargin{1};
    case 0
        [file, pathname] = uigetfile({'*.am'},'Please locate Amira file.','Multiselect','off');
        filename = [pathname file];
end
fid = fopen(filename,'r','ieee-be');
datachar = char(fread(fid,inf,'char'))'; % this will be the char data in the header
switch isempty(findstr(datachar,'BINARY-LITTLE-ENDIAN'))
    case 1
        endianness = 'be';
    case 0
        endianness = 'le';
end
lattice = findstr(datachar,'define Lattice');
parameter = findstr(datachar,'Parameters');
start_data_position = findstr(datachar,'@1');  % note there will be two @1, so take the second one!

%% Determine if we're dealing with a labelfield or regular data
switch isempty(findstr(datachar,'HxByteRLE'))
    case 1  % regular data case
%% blatantly stolen from RTF's geom2am
        data_type = char(datachar(findstr(datachar,'Lattice { ')+10:findstr(datachar,'Data } ')-2));
        switch data_type
            case 'byte'
                out_type = 'int8';
            case 'short'
                out_type = 'int16';
            case 'ushort'
                out_type = 'uint16'; 
            case 'int32'
                out_type = 'int32';
            case 'float'
                out_type = 'float';
            case 'double' 
                out_type = 'double';
        end

%% get the dimensions from the amira header
        fseek(fid,lattice+14,'bof');
        dimensions = str2num(char(fread(fid,parameter-lattice-15,'char'))');
        
        % get the data from the amira file
        fid = fopen(filename,'r',['ieee-' endianness]);
        fseek(fid,start_data_position(2)+2,'bof');
        data = reshape(fread(fid,dimensions(1)*dimensions(2)*dimensions(3),out_type),dimensions(1),dimensions(2),dimensions(3));
%% THIS SPACE JUST SEPARATES THE LABELFIELDS AND TYPICAL AMIRA FILES
%
%   Amira's RLE compression algorithm goes as follows:
%
%   0 to 127 --> Copy the next value n times
%   -127 to -1 -->  Copy the next n values VERBATIM
%
%   The string:
%   00001111111000001010000001000
%   Saved in Amira's RLE as:
%   4 0 7 1 5 0 -125 1 0 1 6 0 -127 1 3 0
%
%   check label2short.m for more info.       
%
%% Thus begins the label2short code.  This will give you a matrix from a labelfield
    case 0 % labelfield case
        lattice = findstr(datachar,'define Lattice');
        parameter = findstr(datachar,'Parameters');
        boundingbox = findstr(datachar,'BoundingBox'); 
        CoordType = findstr(datachar,'CoordType');
        ByteRLE = findstr(datachar,'HxByteRLE,');
        start_data_position = findstr(datachar,'@1');  % note there will be two @1, so take the second one!
        DataSection = findstr(datachar,'# Data section follows');
        
%% get the dimensions from the amira header
        fseek(fid,lattice+14,'bof');
        dimensions = str2num(char(fread(fid,parameter-lattice-15,'char'))');
    
        fseek(fid,ByteRLE+9,'bof');
        byte_length = str2num(char(fread(fid,DataSection - ByteRLE-13,'char'))');

        fseek(fid,boundingbox-1,'bof');
        BB = char(fread(fid,CoordType - boundingbox-5,'char'))';

%% Switch endianess?
        switch endianness
            case 'le'
                fclose(fid);
                fid = fopen(filename,'r',['ieee-' endianness]);
                fseek(fid,start_data_position(2)+2,'bof');
                data = fread(fid,byte_length,'int8');
            case 'be'
                fseek(fid,start_data_position(2)+2,'bof');
                data = fread(fid,byte_length,'int8');
        end

%% This decodes Amira's RLE. I figured this out by trial and error
 
        % pre-allocate space for the output.
        % we're writing everything to a 1-D variable, then reshaping to 3-D.
        % note that we start with a zeros matrix b/c labels contain primarily zeros. we'll avoid writing zeros this way.
        out_line = zeros(dimensions(1)*dimensions(2)*dimensions(3),1);

        % pre-compute the final offset due to -#'s
        final_offset = sum(data(find(data<0))+127);

        net_counter = 0;  % this sifts through the out_line variable
        offset = 0;  % this adds an offset to data_counter when writing VERBATIM

        % this code will search for 0 values, and avoid explicitly writing them.
        for i = 1:floor((byte_length-final_offset)/2) % 
            data_counter = 2*i+offset;  % this will sift through the data variable
            if ((data(data_counter) == 0) && (data(data_counter-1) > 0)) % check if the value is 0 AND that we're not writing VERBATIM
                net_counter = data(data_counter-1)+net_counter;
            else
                if data(data_counter-1) < 0 % check if we're writing VERBATIM
                    for j = 1:128+data(data_counter-1)
                        out_line(j+net_counter) = data(data_counter+j-1);
        %                 disp(['i = ' int2str(i) ', length = ' int2str(data(data_counter-1)) ', written value = ' int2str(data(data_counter+j-1)) ', offset = ' int2str(offset)]);
                    end
                    net_counter = 128+data(data_counter-1)+net_counter;
                    offset = offset + 127+ data(data_counter-1);
                else  % assume that we're copying value data(data_counter), data(data_counter-1) times.
                    for j = 1:data(data_counter-1)
                        out_line(j+net_counter) = data(data_counter);
        %                 disp(['i = ' int2str(i) ', length = ' int2str(data(data_counter-1)) ', written value = ' int2str(data(data_counter))]);
                    end
                    net_counter = data(data_counter-1)+net_counter;
                end
            end
        end 
        clear data;

        % reshape the 1-D to a 3-D
        data = reshape(out_line,dimensions(1),dimensions(2),dimensions(3));
end