function Output = lyngby_read_analyze(filename, arg1, arg2, arg3, ...
    arg4, arg5, arg6, arg7, arg8, arg9, arg10)

% lyngby_read_analyze  - Reads an ANALYZE file
%
%       function V = lyngby_read_analyze(filename, PropertyName, ...
%          PropertyValue) 
%
%       Input:    filename  Filename, with or without the '.img' of
%                           '.hdr' extension
%
%       Property: Output    [ {Volume} | Size | DataType | Endian |
%                           GuessFileNo | Offset | Origin | Scale |
%                           VoxelDim | structure ] 
%                 Size      Volume dimension. 
%                 Datatype  String with precision, for example 'uint16'
%                           or 'float32'. 
%                 FilenameIndex  {1}
%                 Endian    [ {big} | little ] Byteordering
%                 Ordering  [ {xyz} | xzy | yxz | yzx | zxy | zyx ]
%                 Orientation [ {lrpais} | lrapsi | lrapis | ... ]
%                 VoxelMask 
% 
%       Output:   Output    Output according to 'Output' argument in 
%                           the input
%
%       This function reads in an ANALYZE file(s). It is able to
%       return the volume from the .img-file as a vector, and it is
%       able to return some of the information from the .hdr-file.
%
%       The orientation cannot be read from the ANALYZE file. Two
%       different standards exist: The 'original' where right is
%       first in the file, and the 'spm' where left is first in the
%       file.
%
%       SPM type of scale and offset is read and returned and cannot
%       be changed from the properties.
%
%       Example:
%         % If SPM99 is in the path
%         filename = [fileparts(which('spm')) '/templates/T1'];
%
%         % Read in the volume as a vector
%         T1 = lyngby_read_analyze(filename);
%
%         % Read in the size of the volume
%         siz = lyngby_read_analyze(filename, 'output', 'size');
%
%         % Reshape and plot
%         T1 = reshape(T1, siz');
%         figure, imagesc(squeeze(T1(:,:,40)));
%
%       Ref:    http://www.fil.ion.ucl.ac.uk/spm/distrib.html#AzeFmt
%               http://www.mayo.edu/bir/analyze/AnalyzeFileInfo.html
% 
%       See also LYNGBY_WRITE_ANA, LYNGBY_READ_SDT, LYNGBY_READ_VAPET,
%                LYNGBY_READ_ANA4D.
%
% $Id: lyngby_read_analyze.m,v 1.13 2004/03/17 19:51:05 fnielsen Exp $    

    lyngby_global

    
    % Default properties
    output        = 'volume';
    siz           = [];
    datatype      = [];
    endian        = [];
    voxelMask     = 1;
    filenameIndex = 1;
    filenamePattern = []; 
    ordering        = 'xyz';
    orientation     = 'lrpais';  
    origin          = [];
    offset          = [];
    scale           = [];
    voxelDim        = [];
    
    
    % Read Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'output')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if ismember(arg, { 'volume' 'size' 'datatype' 'endian' ...
		'guessfileno' 'offset' 'origin' 'scale' 'voxeldim' ...
		'structure' })
	    output = arg;
	  else
	    error(sprintf(['Argument to ''Output'' not recognized. ' ...
		  'It was: %s'], arg))
	  end
	else
	  error('The argument with ''Output'' should be a string.'); 
	end
      
      elseif strcmp(arg, 'size')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if ~(prod(arg) == 3 )
	  siz = arg;
	else
	  error('The argument with ''Size'' should be a 3 integers.'); 
	end

      elseif strcmp(arg, 'datatype')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  datatype = arg;
	else
	  error('The argument with ''Datatype'' should be a string.'); 
	end

      elseif strcmp(arg, 'endian')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmpi(arg, 'big')
	    endian = arg;
	  elseif strcmpi(arg, 'little')
	    endian = arg;
	  end
	else
	  error('The argument with ''Endian'' should be a string.'); 
	end
	
      elseif strcmp(arg, 'filenameindex')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  filenameIndex = arg;
	else
	  error(['The argument with ''FilenameIndex'' should be ' ...
		'an integers vector.']); 
	end

      elseif strcmp(arg, 'voxelmask')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if issparse(arg) | arg == 1
	  voxelmask = arg;
	else
	  error(['The argument with ''VoxelMask'' should be ' ...
		'either the value ''1'' or a sparse matrix.']); 
	end

      elseif strcmp(arg, 'ordering')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg) & length(arg) == 3 & findstr('x', arg) & ...
	      findstr('y', arg) & findstr('z', arg)
	  ordering = arg;
	else
	  error(['The argument with ''Ordering'' should be ' ...
		'a string with length 3 and consist of x, y and z.']); 
	end

      elseif strcmp(arg, 'orientation')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg) & length(arg) == 6 & ...
	      ( length(findstr('lr', arg)) | length(findstr('rl', arg)) ) & ...
	      ( length(findstr('ap', arg)) | length(findstr('pa', arg)) ) & ...
	      ( length(findstr('is', arg)) | length(findstr('si', arg)) ) 
	  orientation = arg;
	else
	  error(['The argument with ''Orientation'' should be ' ...
		'a string with length 6 and consist of x, y and z.']); 
	end

      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end    
    
    % Handle filename pattern
    if ~isempty(findstr('%', filename))
      filenamePattern = filename;
      filename = sprintf(filename, filenameIndex(1));
    end
    
    % Construct file names
    index = findstr(filename, '.img');
    if ~isempty(index) & length(filename)-3==index
      filename = filename(1:end-4);
    end;  
    index = findstr(filename, '.hdr');
    if ~isempty(index) & length(filename)-3==index
      filename = filename(1:end-4);
    end; 
    filenameHdr = sprintf('%s.hdr', filename);
    filenameImg = sprintf('%s.img', filename);
    
    if (strcmp(output, 'volume') & (isempty(siz) | isempty(datatype) ...
	  | isempty(endian))) | ~strcmp(output, 'volume') 
      % Read header file
      
      % New code to check for endian of header file.
      % --------------------------------------------
      
      % Try to open file
      fid = fopen(filenameHdr, 'rb', 'ieee-be');  % Try big endian first
      if fid == -1 
	error(sprintf('Could not open file: %s', filenameHdr));
      end
      fseek(fid, 0, 'bof');  % Ensure at start
      sizeof_hdr = fread(fid, 1, 'int32'); % Read header length
      if (sizeof_hdr==1543569408),  % Ooops, should be 348! Probably wrong endian
				    % Close and re-open with other endian
        fclose(fid);
	fid = fopen(filenameHdr, 'rb', 'ieee-le'); % Now try little endian
	fseek(fid, 0, 'bof');  % Ensure at start
	sizeof_hdr = fread(fid, 1, 'int32'); % Read header length
	if (sizeof_hdr~=348),
	  % Err, dunno what's wrong!
	  lyngby_ui_message('Unable to open file - header maybe corrupted!');
	  endianHdr = 'unknown';
	end
	endianHdr = 'ieee-le';
      else
	endianHdr = 'ieee-be';
      end
      fseek(fid, 0, 'bof');  % Ensure at start
      sizeof_hdr = fread(fid, 1, 'int32'); % Re-read header length
                                           % to move pointer on....

      % Parse the passed endian arg - not needed for header (as have reasonably
      %  robust testing above), so use for the .img files...
      if isempty(endian)
	  % Empty, so set to same endian as we think the header is
	endianImg = endianHdr;
      elseif strcmpi(endian, 'big')          % Call wants to use big endian
	endianImg = 'ieee-be';
      elseif strcmpi(endian, 'little')       % Call wants to use little endian
	endianImg = 'ieee-le';
      else
	% Unknown endian arg!
	warning('Requested endian unknown - using that of the .hdr file as default!');
      end
      % End of endian code
      % ---------------------------------------------------

      data_type     = fread(fid, 10, 'char');
      db_name       = fread(fid, 18, 'char');
      extents       = fread(fid, 1, 'int32'); 
      session_error = fread(fid, 1, 'int16');  
      regular       = fread(fid, 1, 'uchar'); 			
      if regular ~= 114
	% There is a problem here: Not all socalled "ANALYZE" files seem
	% to contain this field! (Maybe it is only those from VAPET who
	% lacks the field)
	if regular ~= 0
	  if fid > 0,
	    fclose(fid);
	  end
	  warning('Wrong ''Regular'' field in analyze header')
	end
      end
      hkey_un0      = fread(fid,1, 'uchar'); 			
      tmpSiz        = fread(fid,8, 'ushort'); 	% Number of voxels
      rank          = tmpSiz(1)-1;

      if rank < 1 | rank > 4
 	if fid > 0,
 	  fclose(fid);
 	end
 	error(sprintf('Wrong rank field: rank=%d', rank));
      end

       % New code to test for whether file is a real 4D Analyze file
      tmp_dims=tmpSiz(1);
      tmp_x=tmpSiz(2);
      tmp_y=tmpSiz(3);
      tmp_z=tmpSiz(4);
      tmp_t=tmpSiz(5);
      if ((tmp_x>1)&(tmp_y>1)&(tmp_z>1)&(tmp_t>1)&(tmp_dims==4)) % Looks like a 4D file
	disp(['Warning - This looks like a 4D Analyze file - try' ...
			   ' using the ''Analyze 4D'' format instead!']);
      elseif (tmp_dims<2)|(tmp_dims>4)
	warning(sprintf(['Number of dimensions given in the header is' ...
			 ' %d. May not be able to load file!'], tmp_dims));
      end
      % End of new code
	
      if isempty(siz)
        siz = tmpSiz(2:4);
      end
      
      dummy       = fread(fid, 7, 'ushort'); 
      tmpDatatype = fread(fid, 1, 'ushort');
      bitpix      = fread(fid, 1, 'ushort');    % datatype
      dim_un0     = fread(fid, 1, 'ushort');
      pixdim      = fread(fid, 8, 'float32');    % size of pixels
      funused8    = fread(fid, 1, 'float32');
      funused9    = fread(fid, 1, 'float32');
      funused10   = fread(fid, 1, 'float32');
      funused11   = fread(fid, 1, 'float32');
      funused12   = fread(fid, 1, 'float32');
      funused13   = fread(fid, 1, 'float32');
      compressed  = fread(fid, 1, 'float32');
      verified    = fread(fid, 1, 'float32');
      glmax       = fread(fid, 1, 'int32');
      glmin       = fread(fid, 1, 'int32');

      % data_history
      descrip     = fread(fid, 80, 'char');
      aux_file    = fread(fid, 24, 'char');
      orient      = fread(fid, 1, 'char');
      originator  = fread(fid, 10, 'int16');
      generated   = fread(fid, 2, 'int16');
      scannum     = fread(fid, 10, 'char');
      patient_id  = fread(fid, 10, 'char');
      exp_date    = fread(fid, 10, 'char');
      exp_time    = fread(fid, 10, 'char');
      hist_un0    = fread(fid, 3, 'char');
      views       = fread(fid, 1, 'int32');
      vols_added  = fread(fid, 1, 'int32');
      start_field = fread(fid, 1, 'int32');
      field_skip  = fread(fid, 1, 'int32');
      omax        = fread(fid, 1, 'int32');
      omin        = fread(fid, 1, 'int32');
      smax        = fread(fid, 1, 'int32');
      smin        = fread(fid, 1, 'int32');
      
      if fid > 0,
	fclose(fid);
      end
    
    end

    if isempty(datatype)
      type = tmpDatatype * 256 + bitpix;
      if type == hex2dec('0208'),
	datatype = 'uchar';
      elseif type == hex2dec('0410'),
	datatype = 'int16';
      elseif  type == hex2dec('0820'),
	datatype = 'int32';
      elseif  type == hex2dec('1020'),
	datatype = 'float32';
      elseif  type == hex2dec('4040'),
	datatype = 'float64';
      else
	% Unknown type !
	if fid > 0,
	  fclose(fid);
	end
	error('Wrong datatype');
      end
    end
    
    % SPM fields
    if isempty(origin)
      origin = originator(1:3);
    end
    if isempty(offset)
      offset = funused8;
    end
    if isempty(scale)
      scale  = funused9 + ~funused9;
    end
    
    if isempty(voxelDim)
      voxelDim = pixdim(2:rank+1) / 1000;
    end

    
    if strcmp('lr', orientation(1:2))
      x = ((1:siz(1)) - origin(1)) * voxelDim(1); 
    else
      x = ((1:siz(1)) - (siz(1)+1-origin(1))) * voxelDim(1);
    end
    if strcmp('pa', orientation(3:4))
      y = ((1:siz(2)) - origin(2)) * voxelDim(2);
    else
      y = ((1:siz(2)) - (siz(2)+1-origin(2))) * voxelDim(2);
    end
    if strcmp('is', orientation(5:6))
      z = ((1:siz(3)) - origin(3)) * voxelDim(3);
    else
      z = ((1:siz(3)) - (siz(3)+1-origin(3))) * voxelDim(3);
    end

    

    % Return the appropriate field
    if strcmp(output, 'size')
      Output = siz;
    
    elseif strcmp(output, 'datatype')
      Output = datatype;

    elseif strcmp(output, 'endian') % Returns the estimated endian format of
                                    % the hdr file (not img file)
      if strcmpi(endianHdr, 'ieee-be')
	Output = 'big';
      else
	Output = 'little';
      end
    
    elseif strcmp(output, 'guessfileno')

      [path, name, ext, ver] = fileparts(filename);
      filepattern = name(1:min(length(name),5));
      D = dir(path);
      Output = round(length(strmatch(filepattern, strvcat(D.name))) / 2);

    elseif strcmp(output, 'offset')
      Output = offset;
    
    elseif strcmp(output, 'origin')
      Output = origin;
      
    elseif strcmp(output, 'scale')
      Output = scale;
      

    elseif strcmp(output, 'structure')

      % Open '.img' file
      fid = fopen(filenameImg, 'rb', endianImg);
      if fid < 0
	fid = fopen(filenameImgGz, 'rb', endianImg);
	if fid < 0 
	  error(sprintf('File could not be opened: %s', filenameImg));
	end
	fclose(fid);
	error('File is compressed')
      end

      [Volume, count] = fread(fid, prod(siz), datatype);
      fclose(fid);
      if count ~= prod(siz),
	error(sprintf([ ...
	      'Could not read all bytes from: %s\n' ...
	      'Read: %d, Requested: %dx%dx%d = %d'], ...
	    filename, count, siz, prod(siz))); 
      end  

      Volume = reshape(Volume, siz(:)');
      ox = findstr('x', ordering);
      oy = findstr('y', ordering);
      oz = findstr('z', ordering);
      if ~strcmp('xyz', ordering)
	Volume = permute(Volume, [ox oy oz]);
      end
      if strcmp('rl', orientation(1:2))
	Volume = flipdim(Volume, 1);
        origin(1) = length(x) + 1 - origin(1);
      end
      if strcmp('ap', orientation(3:4))
	Volume = flipdim(Volume, 2);
        origin(2) = length(y) + 1 - origin(2);
      end
      if strcmp('si', orientation(5:6))
	Volume = flipdim(Volume, 3);
	origin(3) = length(z) + 1 - origin(3);
      end
      
      if strcmp(output, 'volume')
      
	Output = scale * Volume + offset;
	
      else
	
	Output.volume = scale * Volume + offset;
	xyz = { x y z };
	Output.x = xyz{ox};
	Output.y = xyz{oy};
	Output.z = xyz{oz};
	Output.filename  = filename;
	Output.origin    = origin([ox oy oz]);
	Output.voxelSize = voxelDim([ox oy oz]);
	Output.type      = 'vol';
	
      end
	

    
    elseif strcmp(output, 'volume')

      if isempty(filenamePattern)
	rOutput = length(filenameIndex);
      else 
	rOutput = 1;
      end
      if issparse(voxelMask)
	cOutput = size(voxelMask, 2);
	if size(voxelMask,1) ~= prod(siz)
	  error('Size of volume and ''VoxelMask'' are not equal');
	end
      else
	cOutput = prod(siz);
      end
      Output = zeros(rOutput, cOutput);
      
      
      % Iterate over files
      for n = 1:length(filenameIndex)

	if ~isempty(filenamePattern)
	  filename = sprintf(filenamePattern, filenameIndex(n));
	end

	% Open '.img' file
	fid = fopen(filenameImg, 'rb', endianImg);
	if fid < 0
	  error(sprintf('File could not be opened: %s', filenameImg));
	end
	
	[Volume, count] = fread(fid, prod(siz), datatype);
	if count ~= prod(siz),
	  error(sprintf([ ...
		'Could not read all bytes from: %s\n' ...
		'Read: %d, Requested: %dx%dx%d = %d'], ...
	      filename, count, siz, prod(siz))); 
	end  
	Volume = reshape(Volume, siz(:)');  % (:)' siz must be a row vector 
	if ~strcmp('xyz', ordering)
	  x = findstr('x', ordering);
	  y = findstr('y', ordering);
	  z = findstr('z', ordering);
	  Volume = permute(Volume, [x y z]);
	end
	if strcmp('rl', orientation(1:2))
	  Volume = flipdim(Volume, 1);
	end
	if strcmp('ap', orientation(3:4))
	  Volume = flipdim(Volume, 2);
	end
	if strcmp('si', orientation(5:6))
	  Volume = flipdim(Volume, 3);
	end
	
	Volume = scale * Volume + offset;
	Volume = Volume(:)' * voxelMask;
	
	fclose(fid);
	
	if isempty(filenamePattern)
	  Output = Volume;
	else
	  Output(n, :) = Volume;
	end
	
      end

    elseif strcmp(output, 'voxeldim')
      
      Output = voxelDim;
      
    else
      error('Internal error')
    end



