function lyngby_write_ana(filename, V, siz, arg1, ...
    arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, ...
    arg12, arg13, arg14);

% lyngby_write_ana     - Write a volume to ANALYZE files
%
%        function lyngby_write_ana(filename, V, siz, ...
%           'PropertyName', 'PropertyValue');
%                
%        Input:    filename  Filename with no extension
%                  V         Volume as a row vector
%                  siz       x,y,z, No of pixels in each dimension
%
%        Property: Datatype   [ uchar | uint8 | int16 | {float32} |
%                             float64 ] 
%                  Endian     [ {big} | little ]
%                  Origin     [ {round(siz/2)} ] Center voxel (SPM
%                             field) 
%                  Scale      [ {off} | on ] Scaling of pixel
%                             values. Only applicable for uchar, uint8
%                             and int16  
%                  VoxelDim   {[0.002 0.002 0.002]} Voxel size in
%                             meters 
%
%        SPM type of ANALYZE header is not fully supported.
%        
%        The volume x-coordinate is not flip, ie. it expects a
%        left-to-right volume and writes and left-to-right volume (in
%        accordance with the SPM99 ANALYZE definition NOT the original
%        ANALYZE definition.
%
%        Ref: http://www.fil.ion.ucl.ac.uk/spm/distrib.html#AzeFmt
%             http://www.mayo.edu/bir/analyze/AnalyzeFileInfo.html
%
%        See also LYNGBY_READ_ANALYZE, LYNGBY_READANAVOL,
%                 LYNGBY_READANAHDR, LYNGBY_WRITE_SDT,
%                 LYNGBY_WRITE_VAPET. 
%
% $Id: lyngby_write_ana.m,v 1.11 2003/05/21 09:21:28 fnielsen Exp $


    % Sizes
    [rV, cV] = size(V);
    
    % Check 
    if ~isstr(filename)
      error('filename should be a string');
    end
    if prod(size(siz)) ~= 3
      error('siz should be of size 3');
    end
    if cV ~= prod(siz)
      error(['The number of elements in V should correspond ' ...
	    'to that specified in siz']);
    end

    % Default values
    voxelDim    = [ 0.002 0.002 0.002 ];
    bScale      = 0;
    datatype    = 'float32';
    endian      = 'big';
    description = 'generated with lyngby';
    origin      = [];
    offsetSPM   = 0;
    scaleSPM   = 1;
    
    % Parse Properties
    n = 1;
    while n <= nargin-4
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'voxeldim')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if length(arg) == 3
	  voxelDim = arg;
	else
	  error('Argument with ''VoxelDim'' should be a 3 element vector');
	end
      elseif strcmp(arg, 'precision')
	n = n + 1;
	error('''Precision'' does not apply any longer.'); 
      
      elseif strcmp(arg, 'datatype')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'uchar') | strcmp(arg, 'uint8') | ...
		strcmp(arg, 'int16') | strcmp(arg, 'int32') | ...
		strcmp(arg, 'float32') | strcmp(arg, 'float64')
	    datatype = arg;
	  else
	    error(sprintf('Unrecognized datatype: %s', arg));
	  end
	    
	else
	  error('The argument with ''Datatype'' should be a string.'); 
	end

      elseif strcmp(arg, 'endian')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmpi(arg, 'big') 
	    endian = 'big';
	  elseif strcmpi(arg, 'little')
	    endian = 'little';
	  else
	    error(sprintf('Wrong argument to ''Endian'': %s', arg));
	  end
	else
	  error('The argument with ''Endian'' should be a string.'); 
	end

      elseif strcmp(arg, 'scale')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isstr(arg)
	  if strcmp(lower(arg), 'on') 
	    bScale = 1;
	  end
	else
	  error([ 'The argument with ''Scale'' should ', ...
		'be a string: ''on''  or ''off''']); 
	end
	
      elseif strcmp(arg, 'origin')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg) & length(arg) == 3
	  origin = arg;
	else
	  error([ 'The argument with ''Origin'' should ', ...
		'be a vector with length 3']); 
	end

      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    if strcmpi(endian, 'big')
      sMachineformat = 'ieee-be';
    else
      sMachineformat = 'ieee-le';
    end

    if isempty(origin)
      origin = round(siz/2);
    end
    
    % Flip the X-coordinate (The reference states the ANALYZE files is
    % from right to left.)
    % However, SPM-analyze
    % V = reshape(flipdim(reshape(V, siz(:)'), 1), prod(siz), 1);
    % origin(1) = siz(1) - origin(1) + 1;

    % Scale
    if (bScale)
      if strcmp(datatype, 'uint8') | strcmp(datatype, 'uchar')
	scale = (max(max(V)) - min(min(V))) / 256;
	V = V / scale;
	offset = - min(min(V));
	V = V + offset;
	offset = offset - 0.5;
	index = find(V == 256);
	for i=1:length(index)
	  V(index(i)) = 255;
	end;
      elseif strcmpi(datatype, 'int16')

	offset = 0;
	scale1 = max(max(V)) / 32767;
	scale2 = min(min(V)) / (-32768);
	scale = max([scale1 scale2]);
	V = V / scale;

      else
	warning('scaling ignored');
      end; 

      offsetSPM = offset;
      scaleSPM = scale;
      
    end;
    
    % Discard .img and .hdr file extension
    if strcmp(filename(max(1,end-3):end), '.img')
      filename = filename(1:end-4);
    elseif strcmp(filename(max(1,end-3):end), '.hdr')
      filename = filename(1:end-4);
    end

    % Write the volumefile
    result = 1;
    FileName = sprintf('%s.img', filename);
    fid = fopen(FileName, 'wb', sMachineformat);
      
    if (fid ~= -1)

      fwrite(fid, V, datatype);
	
      if strcmp(datatype, 'uint8')
	anaDatatype = 2;
	bitpix      = 8;
	glmin       = 0;
	glmax       = 255;
	
      elseif strcmp(datatype, 'int16')
	anaDatatype = 4;
	bitpix      = 16;
	glmin       = -32768;
	glmax       = 32767;
	
      elseif strcmp(datatype, 'int32')
	anaDatatype = 8;
	bitpix      = 32;
	glmin       = -2147483648;
	glmax       = 2147483647;
	
      elseif strcmp(datatype, 'float32')
	anaDatatype = 16;
	bitpix      = 32;
	glmin       = -2147483648;
	glmax       = 2147483647;
	
      elseif strcmp(datatype, 'float64')
	anaDatatype = 64;
	bitpix      = 64;
	glmin       = -2147483648;
	glmax       = 2147483647;
	
      else
	error(sprintf('Internal error: datatype=%s', datatype));
      end
	
      
      fclose(fid);
    else
      result = 0;
      error('Not possible to open ANALYZE .img file'); 
    end
    
    % Write the headerfile
    FileName = sprintf('%s.hdr',filename);
    fid = fopen(FileName, 'wb', sMachineformat);

    if fid ~= -1
      fwrite(fid, 348,             'int32');      % sizeof_hdr
      fwrite(fid, zeros(28,1),     'char');       % data_type, db_name
      fwrite(fid, 16384,           'int32');      % extents (16384?)
      fwrite(fid, zeros(2,1),      'char');
      fwrite(fid, 'r',             'char');       % regular
      fwrite(fid, zeros(1,1),      'char');

      fwrite(fid, 4,               'int16');      % dims
      fwrite(fid, siz,             'int16');      % x_dim y_dim z_dim
      fwrite(fid, 1,               'int16');      % t_dim
      fwrite(fid, zeros(3,1),      'int16');      % Rest of dim

      fwrite(fid, zeros(1,7),      'int16');      % funused 
      fwrite(fid, anaDatatype,     'int16');      % datatype
      fwrite(fid, bitpix,          'int16');      % bitpix
      fwrite(fid, 0,               'int16');      % dim_un0

      fwrite(fid, 0,               'float32');    % First pixdim
      fwrite(fid, voxelDim * 1000, 'float32');    % x_size y_size z_size
      fwrite(fid, zeros(4,1),      'float32');    % Rest of pixdim

      fwrite(fid, offsetSPM,       'float32');    % funused 9
      fwrite(fid, scaleSPM,        'float32');    % funused 10
      fwrite(fid, 0,               'float32');    % funused 10 
      fwrite(fid, 0,               'float32');    % funused 11
      fwrite(fid, 0,               'float32');    % funused 12
      fwrite(fid, 0,               'float32');    % funused 13
      fwrite(fid, 0,               'float32');    % compressed
      fwrite(fid, 0,               'float32');    % verified
      fwrite(fid, [glmax glmin],   'int32');      % glmax glmin

      descrip = zeros(80,1);
      indices = 1:min(length(description), 80);
      descrip(indices) = description(indices);
      fwrite(fid, descrip, 'char');               % description

      fwrite(fid, zeros(24,1), 'char');           % aux_file
      fwrite(fid, 0, 'char');                     % orient
      fwrite(fid, [ origin(:)' 0 0 ], 'int16');   % originator, SPM: origin
      fwrite(fid, zeros(85,1), 'char'); 
      
      fclose(fid);
    else
      result = 0;
      error('Not possible to open ANALYZE .hdr file');
    end
    





