function Output = lyngby_read_ana4d(filename, varargin)

% lyngby_read_ana4d    - Reads a 4D dimensional ANALYZE file
%
%       function V = lyngby_read_ana4d(filename, PropertyName, ...
%          PropertyValue) 
%
%       Input:  filename      Filename, with or without the '.img' of
%                             '.hdr' extension
%       Property:
%         Datatype    String with precision, for example 'uint16'
%                     or 'float32'. 
%         Endian      [ {big} | little ]
%         FilenameIndex  {1}
%         Ordering    [ {xyz} | xzy | yxz | yzx | zxy | zyx ]
%         Orientation [ {lrpais} | lrapsi | lrapis | ... ]
%         Output      [ {Volume} | Size | DataType | Endian |
%                     Scans ]
%         Scans       Number of scans
%         Size        Volume dimension. 
%         VolumeIndex {1}
%         VoxelMask 
% 
%       Output: Output	    Output according to 'Output' argument in 
%                           the input
% 
%       Reads a four dimensional ANALYZE file. 
% 
%       See also LYNGBY, LYNGBY_READ_ANALYZE.
%
% $Id: lyngby_read_ana4d.m,v 1.4 2002/03/21 19:01:00 fnielsen Exp $


    % Default properties
    output        = 'volume';
    siz           = [];
    scans         = [];
    datatype      = [];
    endian        = [];
    endianHdr     = [];
    voxelMask     = 1;
    volumeIndex   = 1;
    filenameIndex = 92;
    filenamePattern = []; 
    ordering        = 'xyz';
    orientation     = 'lrpais';  
    
    % Parse Properties
    n = 1;
    while n <= nargin-2
      arg = lower(varargin{n});

      if strcmp(arg, 'datatype')
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg)
	  datatype = arg;
	else
	  error('The argument with ''Datatype'' should be a string.'); 
	end

      elseif strcmp(arg, 'filenameindex')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg)
	  filenameIndex = arg;
	else
	  error(['The argument with ''FilenameIndex'' should be ' ...
		'an integers vector.']); 
	end

      elseif strcmp(arg, 'endian')
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg) &  ( strcmp(arg, 'big') | strcmp(arg, 'little'))
	  endian = arg;
	else
	  error('The argument with ''Endian'' should be a string.'); 
	end

      elseif strcmp(arg, 'endianhdr')
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg)
	  endianHdr = arg;
	else
	  error('The argument with ''EndianHdr'' should be a string.'); 
	end

      elseif strcmp(arg, 'ordering')
	n = n + 1;
	arg = lower(varargin{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;
	arg = lower(varargin{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
      
      elseif strcmp(arg, 'output')
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg)
	  if strcmp(arg, 'volume') | strcmp(arg, 'size') | ...
		strcmp(arg, 'datatype') | strcmp(arg, 'scans') | ...
		strcmp(arg, 'endian')
	    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, 'scans')
	n = n + 1;
	arg = varargin{n};
	if ~(prod(arg) == 1 )
	  scans = arg;
	else
	  error('The argument with ''Scans'' should be one integer.'); 
	end
      
      elseif strcmp(arg, 'size')
	n = n + 1;
	arg = varargin{n};

	if ~(prod(arg) == 3 )
	  siz = arg;
	else
	  error('The argument with ''Output'' should be a 3 integers.'); 
	end

      elseif strcmp(arg, 'volumeindex')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg)
	  volumeIndex = arg;
	else
	  error(['The argument with ''VolumeIndex'' should be ' ...
		'an integers vector.']); 
	end

      elseif strcmp(arg, 'voxelmask')
	n = n + 1;
	arg = varargin{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
	
      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

      % 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 = 'little';  
      else
	endianHdr = 'big';
      end
      fseek(fid, 0, 'bof');  % Ensure at start
      % End of endian code
      % ------------------
      
      sizeof_hdr = fread(fid, 1, 'int32');     % Size of header = 348
      pad1       = fread(fid, 28, 'uchar');    % dummy read header information
      extents    = fread(fid, 1, 'int32'); 
      pad2       = fread(fid, 2, 'uchar');  
      regular    = fread(fid, 1, 'uchar');     % Regular = 114
      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
      pad3       = fread(fid, 1, 'uchar'); 			
      tmpSiz     = fread(fid, 8, 'uint16'); 	% Number of voxels
      rank       = tmpSiz(1)-1;
      if rank < 1 | rank > 4
	if fid > 0,
	  fclose(fid);
	end
	warning(sprintf('Wrong rank field: %d. Setting rank to 4', rank));
	rank = 4;
      end

      % Test for whether file is a real 4D Analyze file
      tmp_dims = tmpSiz(1);  % Num of dimensions
      tmp_x    = tmpSiz(2);     % Num x voxels
      tmp_y    = tmpSiz(3);     % Num y voxels
      tmp_z    = tmpSiz(4);     % Num z voxels
      tmp_t    = tmpSiz(5);     % Num time points
      if tmp_dims~=4
	warning(sprintf([ ...
	      'Wrong dimension field: %d. Assuming dimension is 4' ], ...
	    tmp_dims))
      end
      if ~( tmp_x>1 & tmp_y>1 & tmp_z>1 & tmp_t>1 )
	warning(['This doesn''t look like a 4D Analyze file - try' ...
			   ' using the ''Analyze'' format instead!']);
      end
      if (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
      
      
      if isempty(siz)
	siz = tmpSiz(2:4); % Size of volume, from header (if not 
			   %  specified from function-call 
      end
      if isempty(scans)
	scans = tmpSiz(5); % Number of time-points, from header (if
                           %  not specified from function-call
      end
      
      dummy       = fread(fid, 7, 'uint16'); 
      tmpDatatype = fread(fid, 1, 'uint16');
      bitpix      = fread(fid, 1, 'uint16');    % datatype
      dummy       = fread(fid, 1, 'uint16');
      dummy       = fread(fid, 2, 'uint16');        
      vdim        = fread(fid, 8, 'float32');     % size of pixels
      vdim        = vdim(1:rank)/1000;

% $$$       % Parse tmpDatatype - needs changes to parent functions
% $$$       % before name changes are made....
% $$$       if isempty(datatype),
% $$$ 	switch tmpDatatype
% $$$ 	 case{0}
% $$$ 	  datatype = 'unknown';
% $$$ 	 case{1}
% $$$ 	  datatype = 'binary';  % 1 bit per voxel
% $$$ 	 case{2}
% $$$ 	  datatype = 'uchar';   % Unsigned char, 8 bits
% $$$ 	 case{4}
% $$$ 	  datatype = 'int16';   % 16-bit signed integer
% $$$ 	 case{8}
% $$$ 	  datatype = 'int32';   % 32-bit signed integer
% $$$ 	 case{16}
% $$$ 	  datatype = 'float32'; % 32-bit float
% $$$ 	 case{32}
% $$$ 	  datatype = 'complex'; % 64 bits (2 x 32-bit floats)
% $$$ 	 case{64}
% $$$ 	  datatype = 'float64'; % Double precision (64-bit float)
% $$$ 	 case{128}
% $$$ 	  datatype = 'rgb';
% $$$ 	 case{255}
% $$$ 	  datatype = 'all';
% $$$ 	 otherwise
% $$$ 	  % Unsupported datatype!
% $$$ 	  datatype = 'undefined';
% $$$ 	  error('Unsupported datatype!');
% $$$ 	  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='short';
  	elseif  type == hex2dec('0820'),
  	  datatype='long';
  	elseif  type == hex2dec('1020'),
  	  datatype='float';
  	elseif  type == hex2dec('4040'),
  	  datatype='double';
  	else
  	  % Unknown type !
  	  if fid > 0,
  	    fclose(fid);
  	  end
  	  error('Wrong datatype');
  	end
      end
      
      if fid > 0,
	fclose(fid);
      end
    
    end

    % Set endian of img files
    if isempty(endian)
      % Don't know what endian the data files are stored in - so
      %  use the same as discovered for the header files
      if strcmp(endianHdr, 'little') | strcmp(HDR_ENDIAN, 'big')
	endian = endianHdr;
      else
	% No endian discovered from header - default to big endian
	endian = 'big';
      end	  
    end
    
    
    if strcmp(output, 'size')
      Output = siz;
    
    elseif strcmp(output, 'datatype')
      Output = datatype;
    
    elseif strcmp(output, 'endian')
      Output = endian;
    
    elseif strcmp(output, 'scans')
      Output = scans;
    
    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, '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

	% 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(endian, 'little') 
	  endianFile = 'ieee-le';
	elseif strcmp(endian, 'big')
	  endianFile = 'ieee-be';
	end
	  
	% Open '.img' file
	fid = fopen(filenameImg, 'rb', endianFile);
	if fid < 0 
	  fid = fopen(filename, 'rb', endianFile);
	  if fid < 0
	    error(sprintf('Files could not be opened:\n %s\n %s', ...
		[ filename '.img' ], filename));
	  end
	end

	status = fseek(fid, prod(siz)*(volumeIndex-1) * ...
	    lyngby_datatype2bit(datatype)/8, -1); 
	if status == -1
	  error('Could not read file')
	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);
	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 = Volume(:)' * voxelMask;
	
	fclose(fid);

	if isempty(filenamePattern)
	  Output = Volume;
	else
	  Output(n, :) = Volume;
	end
	
      end
      
    else
      error('Internal error')
    end


