function lyngby_write_vapet(filename, V, siz, arg1, ...
    arg2, arg3, arg4, arg5, arg6, arg7);

% lyngby_write_vapet   - Write a volume to VAPET file 
%
%        function lyngby_write_vapet(filename, V, siz, ...
%                'PropertyName', 'PropertyValue');
%                
%        Input:  filename  Filename with no extension
%                V         Volume
%                siz       x,y,z, No of pixels in each direction
%                PropertyName:
%                   'VoxelDim'   {[0.001 0.0001 0.0001]} Voxel size
%                                in meters
%                   'Precision'  [ 8 | 16 | 32 | {64} ] Precision
%                                (bit per voxel).
%                   'Datatype'   [ {float} | signed | unsigned ] 
%                                If float is off then the default
%                                precision is 16
%                   'Scale'      [ {off} | on ] Scaling of pixel
%                                values. Does not apply in connection
%                                with 'float'.
%
%       Ref: structHeader.h from MVAMC library

% cvs : $Id: lyngby_write_vapet.m,v 1.2 1998/09/18 12:43:44 fnielsen Exp $
%       $Revision: 1.2 $

    if prod(size(V)) ~= prod(siz)
      error(sprintf('V and siz do not correspond. V: %d, siz: %d', ...
	  prod(size(V)), prod(siz)));
    end

    % Default properties
    voxelDim = [ 0.001 0.001 0.001 ];
    precision = 32;
    bScale = 0;
    dataType = 1;
    name = 'unknown';
    scandate = '00/00/00';
    scanstart = '00:00:00';
    
    % Parse Properties
    n = 1;
    while n <= nargin-3
      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;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg) 
	  precision = arg;
	else
	  error('The argument with ''Precision'' should be a number.'); 
	end
      elseif strcmp(arg, 'datatype')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if strcmp(lower(arg), 'float') 
	  dataType = 1;
	else
	  if strcmp(lower(arg), 'signed') 
	    dataType = 2;
	  elseif strcmp(lower(arg), 'unsigned') 
	    dataType = 3;
	  else
	    error(['The argument with ''Datatype'' should be a string: '...
		  '''float'', ''signed'' or ''unsigned''']); 
	  end
	  if ~(precision == 8 | precision == 16)
	    precision = 16;
	  end
	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
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    if dataType == 1
      if precision ~= 32 & precision ~= 64
	error('When using float the precision should be 32 or 64');
      end
    else
      if precision ~= 8 &  precision ~= 16
	error('When not using float the precision should be 8 or 16');
      end
    end
    
    % Scale
    if dataType ~= 1 & bScale
      if precision == 8
	scale = 256 / (max(max(V)) - min(min(V)));
	V = V * scale;
	if dataType == 2
	  scale1 = max(max(V)) / 127;
	  scale2 = min(min(V)) / (-128);
	  scale = max([scale1 scale2]);
	  V = V / scale;
	else
	  offset = - min(min(V));
	  V = V + offset - 128;
	  offset = offset - 0.5;
	  index = find(V == 128);
	  for i=1:length(index)
	    V(index(i)) = 127;
	  end;
	end
      else 
	if dataType == 2
	  scale1 = max(max(V)) / 32767;
	  scale2 = min(min(V)) / (-32768);
	  scale = max([scale1 scale2]);
	  V = V / scale;
	else   
	  scale = 65536 / (max(max(V)) - min(min(V)));
	  V = V * scale;
	  offset = - min(min(V));
	  V = V + offset;
	  offset = offset - 0.5;
	  index = find(V == 65536);
	  for i=1:length(index)
	    V(index(i)) = 65535;
	  end
	end
      end 
    end
  


    % From VAPET structHeader.h ...
    % char    type;           /* type of image, p=pet, m=mri, c=cat */
    % char    name[50];       /* patient id/name */
    % short   rank;           /* number of dimensions (2=image, 3=volume) */
    % short   size[3];        /* size of image */
    % float   min;            /* minimum pixel/voxel value */
    % float   max;            /* minimum pixel/voxel value */
    % float   cmpix[3];       /* pixel/voxel size in cm's */
    % char    datatype;       /* data type, u=unsigned, i=integer, f=float */
    % short   data;           /* size of data in bytes */
    % char    orient[3];      /* image orienation */
    % char    scandate[9];    /* scan date */
    % char    scanstart[9];   /* scan start time */
    
    rank = 3;
    if dataType == 1
      datatype = 'f';
    elseif dataType == 2
      datatype = 'i';
    else
      datatype = 'u';
    end
    data = precision/8;
    orient = 'lr';
    echo = -1;
    tr = -1;
    
    % Header
    s = sprintf('vaphdr\ntype=m\n');       
    s = sprintf('%sname=%s\n', s, name);     
    s = sprintf('%sscandate=%s\n', s, scandate);     
    s = sprintf('%sscanstart=%s\n', s, scanstart);     
    s = sprintf('%srank=%d\n', s, rank);     
    s = sprintf('%ssize=%d %d %d\n', s, siz);     
    s = sprintf('%sdatatype=%s\n', s, datatype);     
    s = sprintf('%sdata=%d\n', s, data);     
    s = sprintf('%scmpix=%f %f %f\n', s, voxelDim*100); 
    s = sprintf('%sorient=%s\n', s, orient); 
    s = sprintf('%smin=%f\n', s, min(min(V))); 
    s = sprintf('%smax=%f\n', s, max(max(V))); 
    s = sprintf('%secho=%f\n', s, echo); 
    s = sprintf('%str=%f\n', s, tr); 

    if length(s) > 512
      error(['The length of the header is too long. ', ...
	    'Try to make for example the name shorter.']);
    else
      s = sprintf('%s%s', s, char(32*ones(512-length(s),1)));
    end
    
    % Filename
    if isempty(findstr(filename, '.'))
      Filename = sprintf('%s.vap',filename);
    else
      Filename = filename;
    end
    
    % Write
    fid = fopen(Filename,'wb','ieee-be');
   
    if (fid ~= -1),
      fwrite(fid, s, 'char');
      if dataType == 1
	if precision == 32
	  count = fwrite(fid, V, 'float32');
	else
	  count = fwrite(fid, V, 'float64');
	end
      elseif dataType == 2;
	if precision == 8
	  count = fwrite(fid, V, 'int8');
	else
	  count = fwrite(fid, V, 'int16');
	end
      else
	if precision == 8
	  count = fwrite(fid, V, 'uint8');
	else
	  count = fwrite(fid, V, 'uint16');
	end
      end
      fclose(fid);
      if prod(size(V)) ~= count
	error(sprintf('Could not write all data, only: %d elements', ...
	    count));
      end
    else
      error('Not possible to open VAPET for writing'); 
    end
      
      























