%  Make avw structure specified by img, dims, voxel_size, and datatype.
%  The header can also include origin, description, and fileprefix.
%
%  Usage: avw = rri_make_avw(img, dims, voxel_size, datatype, ...
%	[origin], [description], [fileprefix])
%
%  Where:
%	img:		volume data, 3D matrix of image data
%	dims:		volume dimension, size of image data
%	voxel_size:	size of each voxel size in millimeter
%	datatype:	storage data type of image data:
%		2 - uint8, 4 - int16, 8 - int32, 16 - float, 64 - double
%	origin (optional):	the AC origin [default: (0,0,0)]
%	description (optional):	description about the image data
%	fileprefix(optional):	a filename string without .img/.hdr extension
%
%  e.g.:
%     dims = [65 87 26]; v_size = [2 2 4]; datatype = 16; origin = [33 44 13];
%     avw = rri_make_avw(img,dims,v_size,datatype,origin);
%
%
%  Part (or all) of this file is copied and modified under GNU license
%  from MRI_TOOLBOX developed by CNSP in Flinders University, Australia
%
function avw = rri_make_avw(varargin)

   avw = [];
   img = double(varargin{1});
   dims = varargin{2};
   voxel_size = varargin{3};
   datatype = varargin{4};
   origin = [0 0 0];
   descrip = '';
   fileprefix = '';

   if nargin > 4
      origin = varargin{5};
   end

   if nargin > 5
      descrip = varargin{6};
   end

   if nargin > 6
      fileprefix = varargin{7};
   end

   maxval = round(max(img(:)));
   minval = round(min(img(:)));

   hdr = make_header(fileprefix, dims, voxel_size, datatype, ...
	origin, descrip, maxval, minval);

   avw.hdr = hdr;
   avw.fileprefix = fileprefix;
   avw.img = img;


   switch double(avw.hdr.dime.bitpix)
   case 1
      avw.img = uint8(avw.img);			% binary
   case 8
      avw.img = uint8(avw.img);			% unsigned char
   case 16
      avw.img = int16(avw.img);			% signed short
   case 32
      if isequal(avw.hdr.dime.datatype, 8), 
         avw.img = int32(avw.img);		% signed integer
      else
         avw.img = single(avw.img);		% float
      end
   case 64
      avw.img = avw.img;
   otherwise
      msgbox('Incorrect image header');
      return;
   end


   return;


%---------------------------------------------------------------------
function hdr = make_header(fileprefix, dims, voxel_size, datatype, ...
	origin, descrip, maxval, minval)

   hdr.hk   = header_key(fileprefix);
   hdr.dime = image_dimension(dims, voxel_size, datatype, maxval, minval);
   hdr.hist = data_history(origin, descrip);
    
   return


%---------------------------------------------------------------------
function hk = header_key(fileprefix)

    hk.sizeof_hdr       = int32(348);   % must be 348!
    hk.data_type(1:10)  = sprintf('%-10s',['dsr      ' 0]);

    db_name = [[fileprefix '.hdr'] '                 '];
    hk.db_name(1:18)    = sprintf('%-18s',[db_name(1:17) 0]);

    hk.extents          = int32(0);
    hk.session_error    = int16(0);
    hk.regular          = sprintf('%1s','r');  % might be uint8
    hk.hkey_un0         = uint8(0);
    
return


%---------------------------------------------------------------------
function dime = image_dimension(dims, voxel_size, datatype, maxval, minval)

   if (length(dims) < 4);
      dims = [dims(:)' 1];
   else
      error('dims should be in form of [x y z]');
      return;
   end;

   if (length(voxel_size) < 4);
      voxel_size = [voxel_size(:)' 0];
   else
      error('voxel_size should be in form of [x y z]');
      return;
   end;
   
   dims = [4 dims 0 0 0];
   pixdim = [0 voxel_size 0 0 0];
   glmin = 0; 

   switch (datatype) 
      case  {1},  bits=1;  glmax = 1;
      case  {2},  bits=8;  glmax = 255;
      case  {4},  bits=16; glmax = 32767;
      case  {8},  bits=32; glmax = (2^31-1);
      case  {16}, bits=32; glmax = maxval; glmin = minval; % glmax = 1;
      case  {64}, bits=64; glmax = maxval; glmin = minval; % glmax = 1;
      otherwise,  
	 error('specified datatype is not supported');
	 return;
   end
    
   dime.dim(1:8)       = int16(dims);
   dime.vox_units(1:4) = sprintf('%-4s',['mm']);
   dime.cal_units(1:8) = sprintf('%8s','');
   dime.unused1        = int16(0);
   dime.datatype       = int16(datatype);
   dime.bitpix         = int16(bits);
   dime.dim_un0        = int16(0);
   dime.pixdim(1:8)    = single(pixdim);
   dime.vox_offset     = single(0);
   % Set default 8bit intensity scale (from MRIcro), otherwise funused3
   dime.roi_scale      = single(1);
   dime.funused1       = single(0);
   dime.funused2       = single(0);
   dime.cal_max        = single(0);
   dime.cal_min        = single(0);
   dime.compressed     = int32(0);
   dime.verified       = int32(0);
   dime.glmax          = int32(glmax);
   dime.glmin          = int32(glmin);

return


%---------------------------------------------------------------------
function hist = data_history(origin, descrip)

   if isempty(origin)
      origin = [0 0 0 0 0];
   else
      origin = [origin(:)' 0 0];
   end;

   descrip = [descrip sprintf('%79s','')];
   hist.descrip(1:80)      = sprintf('%-80s',descrip(1:79));

   hist.aux_file(1:24)     = sprintf('%-24s','');
   hist.orient             = uint8(0); % sprintf( '%1s',''); % see notes in avw_hdr_read
   hist.originator         = int16(origin);
   hist.generated(1:10)    = sprintf('%10s','');
   hist.scannum(1:10)      = sprintf('%10s','');
   hist.patient_id(1:10)   = sprintf('%10s','');
   hist.exp_date(1:10)     = sprintf('%10s','');
   hist.exp_time(1:10)     = sprintf('%10s','');
   hist.hist_un0(1:3)      = sprintf('%3s','');
   hist.views              = int32(0);
   hist.vols_added         = int32(0);
   hist.start_field        = int32(0);
   hist.field_skip         = int32(0);
   hist.omax               = int32(0);
   hist.omin               = int32(0);
   hist.smax               = int32(0);
   hist.smin               = int32(0);
	
return

