function Output = lyngby_read_vapet(filename, arg1, arg2, arg3, ...
    arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14);

% lyngby_read_vapet    - Read a 3D VAPET file 
%
%	function Output = lyngby_read_vapet(filename,
%	   'PropertyName', 'PropertyValue');
%
%	Input:	filename     Filename
%               Property:
%                 Output    [ {Volume} | DataType | Endian | Size |
%                           VoxelDim | Orientation ]
%                 Size      Volume dimension. 
%                 Datatype  String with precision, for example 'uint16'
%                           or 'float32'. 
%                 FilenameIndex  {1}
%                 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 a VAPET volume file.
%
%       See also: lyngby_read_vapet4d, lyngby_read_analyze,
%                 lyngby_write_ana

% cvs : $Id: lyngby_read_vapet.m,v 1.1 1999/10/01 13:45:46 fnielsen Exp $
%       $Revision: 1.1 $


    % Default properties
    output        = 'volume';
    siz           = [];
    hdrsz         = 512;
    datatype      = [];
    endian        = 'big';
    voxelMask     = 1;
    volumeIndex   = 1;
    filenameIndex = 1;
    filenamePattern = []; 
    ordering        = 'xyz';
    orientation     = '';  
    voxelDim        = [];
    
    % Parse 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 strcmp(arg, 'volume') | strcmp(arg, 'size') | ...
		strcmp(arg, 'orientation') | strcmp(arg, 'datatype') ...
		| strcmp(arg, 'voxelmask') | strcmp(arg, 'voxeldim')
	    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 ''Output'' 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, '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

    % Big endian little endian
    if strcmpi(endian, 'big')
      sMachineformat = 'ieee-be';
    elseif strcmpi(endian, 'little')
      sMachineformat = 'ieee-le';
    else
      error('Internal error');
    end

    % Handle filename pattern
    if ~isempty(findstr('%', filename))
      filenamePattern = filename;
      filename = sprintf(filename, filenameIndex(1));
    end
    
    
    % Open file
    fid = fopen(filename, 'rb', sMachineformat);
    if fid < 0, 
      error(sprintf('Could not open file: %s ', filename));
    end

    line = fgetl(fid);
    if ~strcmp(line,'vaphdr')
      if fid > 0, 
	fclose(fid);
      end
      error('File not a VAPET: Could not find ''vaphdr''');
    end
    
    fseek(fid, 0, -1);        % Rewind
    header = char(fread(fid, 512, 'char')');

    % Header size
    hdrsz = sscanf(header((findstr(header, 'hdrsz=')+6):512), '%d', 1);
    if isempty(hdrsz)
      % Probably default 512
      hdrsz = 512;
    end
    if hdrsz ~= 512
      fseek(fid, 0, -1);                           % Rewind
      header = char(fread(fid, hdrsz, 'char')');   % Reread
    end

    % Rank
    rank = sscanf(header((findstr(header, 'rank=')+5):hdrsz), '%d', 1); 
    if prod(size(rank)) ~=1 | rank < 1 | rank > 4,
      if fid > 0, 
	fclose(fid);
      end
      error(sprintf('Wrong rank in file: %d', rank));
    end
    
    % Size
    tmpSiz = sscanf(header((findstr(header, 'size=')+5):hdrsz), ...
	'%d', rank);
    if isempty(siz)
      if prod(size(tmpSiz)) == rank
	siz = tmpSiz;
      else
	error('The field ''size'' does not correspond to rank');
      end
    end
    
    % Voxel dimensions
    vdim = sscanf(header((findstr(header, 'cmpix=')+6):hdrsz), ...
	'%f', rank); 
    if isempty(voxelDim)
      if prod(size(vdim)) == rank
	voxelDim = vdim/100;
      else
	error('The field ''cmpix'' is wrong');
      end
    end

    % Orientation
    index = findstr(header, 'orient=') + 7;
    orient = header(index:index+1); 

    if isempty(orientation)
      if strcmp(orient, 'lr') | strcmp(orient, 'rl')
	orientation = [ orient 'ap' 'is'];
      else
	error('The field ''orient'' is wrong');
      end
    end

    
    % Datatype
    bytes = sscanf(header((findstr(header, 'data=')+5):hdrsz), ...
	'%c', 1); 
    prec = sscanf(header((findstr(header, 'datatype=')+9):hdrsz), ...
	'%c', 1);
    type = [bytes prec];
    if isempty(datatype)
      if prod(size(type)) == 2,
	if strcmp(type,'1u'),
	  datatype = 'uint8';
	elseif strcmp(type, '2u'),
	  datatype = 'uint16';
	elseif strcmp(type, '1i'),
	  datatype = 'int8';
	elseif strcmp(type, '2i'),
	  datatype = 'int16';
	elseif strcmp(type, '4f'),
	  datatype = 'float32';
	elseif strcmp(type, '8f'),
	  datatype = 'float64';
	else
	  error('Unknown datatype');
	end
      else
	if fid > 0, 
	  fclose(fid);
	end
	error('Wrong datatype')
      end
    end

    if ~strcmp(output, 'volume') & fid > 0 
      fclose(fid);
    end


    if strcmp(output, 'size')
      Output = siz;

    elseif strcmp(output, 'orientation')
      Output = orientation;
   
    elseif strcmp(output, 'datatype')
      Output = datatype;
    
    elseif strcmp(output, 'voxelmask')
      
      Output = 1;
      
    elseif strcmp(output, 'voxeldim')
      
      Output = voxelDim;
          
    elseif strcmp(output, 'orientation')
      
      Output = orientation;
          
   
    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')


      [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(:)');
      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
      
      Output = Volume(:)';
    
    else 
      error(sprintf('Internal error: output=%s', output))

    end









