function count = otlwrite(X, FileName, H )
% count = otlwrite(X, FileName, H )
%
% Write data X and header H to an OTL output file.
%
% count (scalar) returns the number of bytes written, including header and data.
%
% X [NDIMxNPTS] - coordinates of the N outline points.  
%       NDIM is the number of dimensions.
% H (scalar struct) is the header, which may contain any of the fields returned by otlread.
% H is an optional argument; if you don't provide H, or if any of the fields of H is nonexistant,
% then otlwrite will either fill in reasonable defaults, or it will just not use that field.
%
%   Header fields with default values:  If any of the following is not specified,
%   a reasonable default will be used instead:
%       H.offset [1xNDIM] - the minimum value of X in each dimension.
%       H.npixels [1xNDIM] - the number of possible values of X in each dimension
%          WARNING: DATA VALUES LESS THAN H.offset, OR MORE THAN H.offset+H.npixels,
%          WILL BE CLIPPED!!!
%       H.scaling_factor [1xNDIM] - scaling factor to use in writing X.
%          If H.scaling_factor is too large to fit the number of bits, it will be reduced.
%       H.sample_n_bytes [scalar] - the number of bytes which should be used to encode each X element
%
%   Optional header fields:  If any of the following is not specified,
%   the field will not be saved to the output file.
%       H.slice_thickness [scalar] - thickness of the MRI or CT image slice which X represents
%       H.top_left_hand_corner [1x3] - absolute 3D coordinates of the top left hand corner
%             of the image slice which X represents.
%       H.top_right_hand_corner [1x3] - absolute 3D coordinates of the top right hand corner
%             of the image slice which X represents.
%       H.bottom_right_hand_corner [1x3] - absolute 3D coordinates of the bottom right hand corner
%             of the image slice which X represents.
% 
% If FileName already exists, it is backed up in a file
% with extension '.bak'.  The backup is not backed up!
%
% Example:
%  Suppose X is the data matrix for an outline, where the image starts at position [0.5,0.5]
%  and goes to position [mrhd.ih.height+0.5, mrhd.ih.width+0.5].
%  The colors of each point are given by COLORS, with values between 0 and 128.
%  The thickness and corners are given by mrhd.mri.slice_thickness and mrhd.mri.corners, in the
%  format returned by geread.m.
%  This outline can be written using
%
%  > otlwrite(FileName, [ X; double(COLORS) ], struct( ...
%      'npixels', [mrhd.ih.width, mrhd.ih.height, 128 ], ...
%      'offset', [0.5 0.5 0], ...
%      'slice_thickness', mrhd.mri.slice_thickness, ...
%      'top_left_hand_corner', mrhd.mri.corners(1,:), ...
%      'top_right_hand_corner', mrhd.mri.corners(2,:), ...
%      'bottom_right_hand_corner', mrhd.mri.corners(4,:), ...
%      'sample_n_bytes', 2));
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Required bits of information
% Set the output byte format: this must be '01', since otlwrite doesn't know any other format
if nargin<3,  H = struct('sample_byte_format', '01'); end
H.sample_byte_format = '01';

% Set the number of bytes
if ~isfield(H,'sample_n_bytes'), H.sample_n_bytes = 2; end

% Set default image offset and size
offset = min(X,[],2)';
if ~isfield(H,'offset'), 
   H.offset = offset; 
elseif length(H.offset)<length(offset),
   H.offset = [ H.offset offset((length(H.offset)+1):length(offset)) ];
end
npixels = ceil(max(X,[],2)-min(X,[],2))' + 1;
if ~isfield(H,'npixels'), 
   H.npixels = npixels;
elseif length(H.npixels)<length(npixels),
   H.npixels = [ H.npixels npixels((length(H.npixels)+1):length(npixels)) ];
end
% Arbitrary: assume that there is no such thing as an image of less than 8 pixels
H.npixels = max(8,H.npixels);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set database structure values, if not specified
if ~isfield(H,'database_id'), H.database_id= 'otlwrite'; end
if ~isfield(H,'database_version'), H.database_version = '1.0'; end
if ~isfield(H,'creation_date'), H.creation_date = date; end
if ~isfield(H,'image_id'),
  r = FileName;
  if length(findstr(filesep, FileName)) == 0, r = [pwd '/' FileName]; end
  ncomps = 0;
  while length(r)>0,
    ncomps = ncomps+1;
    [comps{ncomps}, r] = strtok(r, '/');
  end
  if ncomps > 1,
    H.image_id = fullfile(comps{max(1,ncomps-5):ncomps});
  elseif ncomps==1,
    H.image_id = comps{1};
  else,
    H.image_id = '--'
  end
end

%%%%%%%%%%%%%%%%%%%%
% If the number of points is greater than zero, then figure out that information, and add it to H
[NDIM, NPTS] = size(X);
H.channel_count = NDIM;
H.sample_count = NPTS;

if NPTS > 0,
  scaling_factor = pow2( 8 * H.sample_n_bytes - 1 ) ./ H.npixels;
  if isfield(H,'scaling_factor'), 
     if length(H.scaling_factor)<length(scaling_factor),
        H.scaling_factor = [ H.scaling_factor scaling_factor((length(H.scaling_factor)+1):length(scaling_factor)) ];
     end
     scaling_factor = min(scaling_factor, H.scaling_factor(1:length(scaling_factor)));
  end;
  H.scaling_factor = scaling_factor;
  
  for dim=1:NDIM,
    XOUT(dim,:) = H.scaling_factor(dim) * max(0, min(H.npixels(dim), ( X(dim,:) - H.offset(dim) )));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%
%   If the output filename already exists, back it up
f = dir(FileName);
if prod(size(dir(FileName))) > 0,
  ind = findstr(FileName, '.');
  BackupFile = [ FileName(1:ind(length(ind))) 'bak' ];
  if strcmp(filesep,'/'), 
    moved = unix(['mv ' FileName ' ' BackupFile ]); 
    % if moved==0, disp([ 'Backed up old ' FileName ' as ' BackupFile]), end;
 elseif strcmp(filesep,'\'), 
    moved=dos(['copy ' FileName ' ' BackupFile ]); 
    % if moved==0, disp([ 'Backed up old ' FileName ' as ' BackupFile]), end;
  end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try to open the output file
[fid, message] = fopen(FileName, 'w', 'l');
if fid <= 0,
  disp(sprintf('otlwrite: Unable to open file %s for writing, message %s',FileName, message));
  count = 0;
  return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pull out cell arrays of fields and values, for the nistwr function
fields = fieldnames(H);
for i=1:length(fields),
  values{i} = getfield(H, fields{i});
end

% Write out a NIST header
header_length = nistwr(fid,fields,values,0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Write out the data with the number of bytes specified by H.sample_n_bytes
if NPTS>0,
  data_length = H.sample_n_bytes * fwrite(fid, XOUT, sprintf('integer*%d', H.sample_n_bytes));
else
  data_length = 0;
end

% Finish up
count = header_length + data_length;
fclose(fid);
% disp(sprintf('Saved outline to %s in NIST_1A format with a %d-byte header', FileName, header_length) );
      
