function [siz, vdim, datatype] = lyngby_readanahdr(file);

% lyngby_readanahdr    - Reads an ANALYZE header
%
%       function [siz, vdim, datatype] = lyngby_readanahdr(file);
%
%       Input:  file        Filename, with or without the '.img' of
%                          '.hdr' extension
%
%       Output: siz         Dimensions of volume (number of voxels)
%               vdim        Physical size of a voxel in meter
%               datatype    String with precision. String is
%                           compatible with fread, e.g., 'uchar'
%
%       This function returns some of the information contained in the
%       '.hdr' file of the ANALYZE volume file format. This
%       information can be used to read the actual volume data with
%       lyngby_readanavol.
%
%       The "regular" field is checked for the values 114 OR 0.
%
%       Ref: Roger Woods AIR-program "src/HEADER.h"
%
%       See also: lyngby_read_volume, lyngby_readanavol, fread

% cvs : $Id: lyngby_readanahdr.m,v 1.8 2001/07/19 08:12:15 mgl Exp $
%       $Revision: 1.8 $

% $$$   % Old file-opening code
% $$$     filename = [file '.hdr']; 
% $$$     fid = fopen(filename, 'rb', 'ieee-be');
% $$$     if fid < 0 
% $$$       index = findstr(file, '.img');
% $$$       if ~isempty(index) & index == (length(file) - 3)
% $$$ 	filename = [file(1:index) 'hdr'];
% $$$ 	fid = fopen(filename ,'rb', 'ieee-be');
% $$$ 	if fid < 0
% $$$ 	  error(sprintf('Files could not be opened: %s, %s', ...
% $$$ 	      [ file '.hdr' ], file, filename));
% $$$ 	end
% $$$       else
% $$$ 	filename = file;
% $$$ 	fid = fopen(filename, 'rb', 'ieee-be');
% $$$ 	if fid < 0
% $$$ 	  error(sprintf('Files could not be opened: %s, %s', ...
% $$$ 	      [ file '.hdr' ], file));
% $$$ 	end
% $$$       end
% $$$     end

  % New file-opening code
    % Parse filename to check for extensions
    index = findstr(file, '.img'); % Check for .img ending
    index2= findstr(file, '.hdr'); % Check for .hdr ending
    if ~isempty(index) & index == (length(file) - 3)  % Found .img at end
      filename = [file(1:index) 'hdr'];    % Strip and add .hdr
    else
      if ~isempty(index2) & index2 == (length(file) - 3)  % Found .hdr at end
	filename = file;  % Just use it as given
      else
	filename = [file '.hdr'];
      end
    end
    
    % Try to open file
    fid=fopen(filename, 'rb', 'ieee-be');  % Try big endian first
    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(filename, '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!
	error('Unable to open file - header corrupted?');
      end
      % Maybe set global endian variable here
    end
    fseek(fid, 0, 'bof');  % Ensure at start
  
  % Original code
    sizeof_hdr = fread(fid, 1, 'long');
    pad1       = fread(fid, 28, 'uchar');    % dummy read header information
    extents    = fread(fid, 1, 'int'); 
    pad2       = fread(fid, 2, 'uchar');  
    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
	vdim=[]; datatype=''; siz=[];
	if fid > 0,
	  fclose(fid);
	end
	return
      end
    end
    pad3       = fread(fid,1,'uchar'); 			
    dim        = fread(fid,4,'ushort'); 		% dimension, number of pixels
    rank       = dim(1)-1;
    if rank < 1 | rank > 4
      vdim=[]; datatype=''; siz=[];
      if fid > 0,
	fclose(fid);
      end
      return
    end
    siz = dim(2:rank+1);
    fread(fid,11,'ushort'); 
    datatype = fread(fid,1,'ushort');
    bitpix = fread(fid,1,'ushort'); 		% datatype
    fread(fid,1,'ushort');
    fread(fid,2,'ushort');        
    vdim = fread(fid,8,'float');  % size of pixels
    vdim = vdim(1:rank)/1000;
    
    type = datatype*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 !
      vdim = []; datatype = ''; siz = [];
      if fid > 0,
	fclose(fid);
      end
      return;
    end
    if fid > 0,
      fclose(fid);
    end
    
    








