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 = amiraReadleave the argument blank and use the file selector
%
%%

warning off all;

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'))';
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));
%%
%
%
%
%
%
%
%
%
%       
%% 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, not
%  sure what algorithm they use.
 
        % pre-allocate space for the output
        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;
        offset = 0;
        for i = 1:floor((byte_length-final_offset)/2) % 
        %     disp(int2str(i));
            data_counter = 2*i+offset;
            if ((data(data_counter) == 0) && (data(data_counter-1) > 0))
                net_counter = data(data_counter-1)+net_counter;
            else
                if data(data_counter-1) < 0
                    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
                    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;
        data = reshape(out_line,dimensions(1),dimensions(2),dimensions(3));
end