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

% lyngby_read_sdt      - Read a Stimulate SDT/SPR file
%
%       function V = lyngby_read_sdt(filename, PropertyName, ...
%          PropertyValue) 
%
%       Input:    filename  Filename, with or without the '.spr' of
%                           '.sdt' extension
%       Property: Output    [ {Volume} | Size | DataType | Endian |
%                           Scans ]
%                 Size      Volume dimension. 
%                 Datatype  String with precision, for example 'uint16'
%                           or 'float32'. 
%                 VolumeIndex   {1}
%                 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 functions read a Stimulate SDT (data) and SPR (header)
%       file. It is not able to handle write the COMPLEX datatype. 
%
%       Ref: Stimulate, http://www.cmrr.umn.edu/stimulate/
%            http://www.cmrr.umn.edu/stimulate/stimUsersGuide/
%
%       See also LYNGBY, LYNGBY_WRITE_SDT, LYNGBY_READ_ANALYZE,
%                LYNGBY_READ_VAPET. 
%
% $Id: lyngby_read_sdt.m,v 1.3 2003/03/04 15:38:10 fnielsen Exp $



    % Default properties
    output        = 'volume';
    siz           = [];
    datatype      = [];
    endian        = 'big';
    voxelMask     = 1;
    volumeIndex   = 1;
    filenameIndex = 1;
    filenamePattern = []; 
    ordering        = 'xyz';
    orientation       = 'lrpais';  
    real2WordScale  = 1;
    
    % 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, 'datatype') | strcmp(arg, 'scans')
	    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(size(arg)) == 3
	  siz = arg;
	else
	  error('The argument with ''Output'' should be a 3 integers.'); 
	end

      elseif strcmp(arg, 'scans')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if length(arg) == 1
	  scans = arg;
	else
	  error('The argument with ''Scans'' 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, 'volumeindex')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  volumeIndex = arg;
	else
	  error(['The argument with ''VolumeIndex'' should be ' ...
		'an integers vector.']); 
	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
      
    
    if (strcmp(output, 'volume') & (isempty(siz) | isempty(datatype) ...
	  | isempty(endian))) | ~strcmp(output, 'volume')
      % Read header file
   
      % Open header file
      filename2 = [filename '.spr']; 
      fid = fopen(filename2, 'rt');
      if fid < 0 
	index = findstr(filename, '.sdt');
	if ~isempty(index) & index == (length(filename) - 3)
	  filename2 = [filename(1:index) 'spr'];
	  fid = fopen(filename2 ,'rt');
	  if fid < 0
	    error(sprintf('Files could not be opened:\n %s\n %s', ...
		filename, filename2));
	  end
	else
	  fid = fopen(filename, 'rt');
	  if fid < 0
	    error(sprintf('Files could not be opened:\n %s\n %s', ...
		[ filename '.spr' ], filename));
	  end
	end
      end

      % Read fields in SPR file
      cont = 1;
      while cont
	line = fgetl(fid);
	if ~isstr(line)
	  break 
	end
	if strcmp('numDim', line(1:6))
	  numDim = str2num(line(8:end));
	elseif strcmp('dim', line(1:3))
	  dim = str2num(line(5:end));
	elseif strcmp('origin', line(1:6))
	  origin = str2num(line(8:end));
	elseif strcmp('extent', line(1:6))
	  extent = str2num(line(8:end));
	elseif strcmp('fov', line(1:3))
	  fov = str2num(line(5:end));
	elseif strcmp('interval', line(1:8))
	  interval = str2num(line(10:end));
	elseif strcmp('dataType', line(1:8))
	  dataType = strtok(line(10:end));
	elseif strcmp('displayRange', line(1:12))
	  displayRange = str2num(line(14:end));
	elseif strcmp('fidName', line(1:7))
	  fidName = strtok(line(9:end));
	elseif strcmp('sdtOrient', line(1:9))
	  sdtOrient = strtok(line(11:end));
	elseif strcmp('Real2WordScale', line(1:14))
	  real2WordScale = strtok(line(15:end));

	else
	  warning(sprintf('Unrecognized field in SPR file: %s.', line));
	end
      end
      if fid > 0,
	fclose(fid);
      end
    end

    
    % Datatype
    if isempty(datatype)
      if strcmp('BYTE', dataType)
	datatype = 'char';
      elseif strcmp('WORD', dataType)
	datatype = 'short';
      elseif strcmp('LWORD', dataType)
	datatype = 'int';
      elseif strcmp('REAL', dataType)
	datatype = 'float';
      elseif strcmp('COMPLEX', dataType)
	datatype = 'float';
	error('dataType ''COMPLEX'' not supported');
      else
	error(sprintf('Unrecognized datatype: %s', dataType));
      end
    end
    
    % Ordering
    if isempty(ordering)
      if strcmp('ax', sdtOrient)
	ordering = 'xyz';
      elseif strcmp('sag', sdtOrient)
	ordering = 'yzx';
      else
	ordering = 'xyz';
	warning(sprintf('sdtOrient field not recognized: %s', ...
	    sdtOrient));
      end
    end
    
    
    
    if strcmp(output, 'size')
      Output = dim(1:3);
    
    elseif strcmp(output, 'datatype')
      Output = datatype;

    elseif strcmp(output, 'scans')

      Output = dim(numDim);
    
    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 '.sdt' file
	filename2 = [ filename '.sdt' ];
	fid = fopen(filename2, 'rb', 'ieee-be');
	if fid < 0 
	  index = findstr(filename, '.spr');
	  if ~isempty(index) & index == (length(filename) - 3)
	    filename2 = [filename(1:index) 'sdt'];
	    fid = fopen(filename2 ,'rb', 'ieee-be');
	    if fid < 0
	      error(sprintf('Files could not be opened:\n %s\n %s', ...
		  [ filename '.sdt' ], filename2));
	    end
	  else
	    fid = fopen(filename, 'rb', 'ieee-be');
	    if fid < 0
	      error(sprintf('Files could not be opened:\n %s\n %s', ...
		  [ filename '.sdt' ], filename));
	    end
	  end
	end

	% fread(fid, prod(siz)*(volumeIndex-1), datatype);
	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;

	if real2WordScale ~= 1
	  Volume = Volume / real2WordScale;
	end
	
	fclose(fid);

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

      end
      
    else
      error('Internal error')
    end









