function [ X_cell, H_cell, fields_cell, values_cell ] = otlread(FileList)
% [ X, H, fields, values ] = otlread(FileName)
%
% Read x, y, z, etc. points from an OTL-format file.
% OTL-format file may be either a matlab '.mat'-format file,
% or a NIST-format file with little-endian integer data, 
% and the following header items:
%
%   npixels - a string of number showing max pixels in each direction, OR
%        npixels_x - number of image pixels in x direction (default 256)
%        npixels_y - number of image pixels in y direction (default 256)
%   channel_count - number of data channels (default 2: x and y)
%   sample_n_bytes - number of bytes per sample (default 2)
%
% X [NDIMxNPTS] is a list of the outline points.  NDIM is the number of dimensions,
%      NPTS is the number of points which exist in the outline.
%
% H (struct) is a header structure containing all of the same field/value pairs
%      as are given in the file header, plus default values of npixels, channel_count,
%      sample_n_bytes, and scaling_factor, if those fields were not defined
%   
% If fields and values are requested, they return the fields and values
% read from the NIST header, if it exists.
%
%
% SIMULTANEOUSLY READING MULTIPLE FILES:
%
% [ X_cell, H_cell, fields_cell, values_cell ] = otlread( FileList )
%
% FileList is a character array, e.g. FileList = strvcat('file1.otl','file2.otl');
%
% X_cell, H_cell, fields_cell, and values_cell are one-dimensional cell arrays whose elements
% are the values X, H, fields, and values returned from each individual otl file.
%

% Check to see how many files were specified
[NFILES, LEN] = size(FileList);

% Load each of the specified files
for ifile=1:NFILES,
  FileName = deblank(FileList(ifile,:));
  fid = fopen(FileName, 'r', 'l');
  if fid > 0,
    [ fields, values ] = nistread(fid);
    if length(fields) > 0,
    
      % First, create the header structure
      H = struct(fields{1}, values{1});
      for i=2:length(fields),
        H = setfield(H, fields{i}, values{i});
      end
    
      % Check to see if sample_n_bytes is specified; if not, assume short
      if ~isfield(H,'sample_n_bytes'), H.sample_n_bytes = 2; end
      fmt = sprintf('int%d', 8*H.sample_n_bytes);
    
      % Look for channel_count, sample_count, sample_n_bytes, and offset
      if ~isfield(H,'channel_count'), H.channel_count = 2; end
      if ~isfield(H,'sample_count'), H.channel_count = Inf; end
      if ~isfield(H,'sample_n_bytes'), H.sample_n_bytes = 2; end

      % Look for npixels (new format file) or npixels_x and npixels_y (old format file)
      if ~isfield(H,'npixels'), 
        if ~isfield(H,'npixels_x'), H.npixels_x = 256; end
        if ~isfield(H,'npixels_y'), H.npixels_y = 256; end
        H.npixels = [ H.npixels_x, H.npixels_y ];
      elseif ischar(H.npixels), 
        H.npixels = str2num(H.npixels); 
      end
      ls = length(H.npixels);
      if ls < H.channel_count,
        H.npixels = [ H.npixels, repmat(H.npixels(ls), [1, H.channel_count-ls]) ];
      end

      % Look for offset; if none is found, assume that it is [0.5, 0.5, ...] (old format)
      if isfield(H, 'offset'),
	if ischar(H.offset), H.offset = str2num(H.offset); end
      else
	H.offset = repmat(0.5, 1, H.channel_count);
      end
      
      % Look for scaling_factor
      if ~isfield(H, 'scaling_factor'),
        % If scaling_factor doesn't exist, fill in a default
        H.scaling_factor = repmat(pow2(8*H.sample_n_bytes-1) / max(H.npixels), 1, 2);
      elseif ischar(H.scaling_factor),
        % If scaling_factor is a string, try converting it to a vector.
        H.scaling_factor = str2num(H.scaling_factor); 
      elseif length(H.scaling_factor)==1,
        % If scaling_factor is a scalar, duplicate it
        H.scaling_factor = [ H.scaling_factor H.scaling_factor];
      end
      if length(H.scaling_factor) < H.channel_count,
        % If scaling_factor is too short, append ones.
        H.scaling_factor = [ H.scaling_factor, ones(1, H.channel_count - length(H.scaling_factor)) ];
      end
    
      % Read the data; normalize as specified in H.scaling_factor and H.offset
      X = fread(fid, [H.channel_count Inf], fmt);
      [truechan, H.sample_count] = size(X);
      for chan=1:truechan,
        X(chan,:) = X(chan,:)/H.scaling_factor(chan) + H.offset(chan);
      end
    
      % Finally, close the input file!!
      fclose(fid);
    else,
      % load it from .mat format file
      %
      fclose(fid);
      load(FileName,'-mat','x','y');
      X = [x(:)'; y(:)'; repmat(double('g'), size(x(:)')) ];
      
      % Output a default value for the header
      H = struct('channel_count', 3, 'sample_count', length(x), 'npixels', [256 256 128], ...  
  	  'scaling_factor', [1 1 1]);
    end;		% if length(fields)>0
  else,
    % If we can't open the file, assign some reasonable default values
    %
    disp(sprintf('otlread: Unable to load outline from file %', FileName))
    X = zeros(2,0);
    H = struct('channel_count',2,'sample_count',0);
    fields = {};
    values = {};
  end;		% if fid>0

  % Load the obtained values into cell arrays, if there are many files
  %
  if NFILES>1,
    X_cell{ifile} = X;
    H_cell{ifile} = H;
    fields_cell{ifile} = fields;
    values_cell{ifile} = values;
  else,
    X_cell = X;
    H_cell = H;
    fields_cell = fields;
    values_cell = values;
  end
end			% for ifile=1:NFILES


