#include "mrilib.h"
#include "thd_niftiwrite.h"

/* prototypes */

nifti_image * populate_nifti_image(THD_3dim_dataset *dset, niftiwr_opts_t options) ;

void nifti_set_afni_extension(THD_3dim_dataset *dset,nifti_image *nim) ;

static int get_slice_timing_pattern( float * times, int len, float * delta );

/*******************************************************************/
/*!  Write an AFNI dataset as a NIfTI file.
     - dset  = AFNI dataset
     - options = structure with options to control the output

   Return value is 1 if went OK, 0 if not.
-------------------------------------------------------------------*/

int THD_write_nifti( THD_3dim_dataset *dset, niftiwr_opts_t options )
{
  nifti_image *nim ;
  nifti_brick_list nbl ;
  int ii ;
  char *fname , *cpt ;

ENTRY("THD_write_nifti") ;

  nifti_set_debug_level(options.debug_level) ;

   /*-- check inputs for goodness --*/

  fname = nifti_strdup(options.infile_name);

  if( !THD_filename_ok(fname) || fname[0] == '-' ){
    ERROR_message("Illegal filename for NIfTI output: %s\n",
                  (fname != NULL) ? fname : "(null)" ) ;
    RETURN(0) ;
  }

  if( !ISVALID_DSET(dset) ){
    ERROR_message("Illegal input dataset for NIfTI output: %s\n", fname) ;
    RETURN(0) ;
  }

  /*-- load dataset from disk, if need be --*/

  DSET_load(dset) ;
  if( !DSET_LOADED(dset) ){
    ERROR_message(
            "Can't write NIfTI file since dataset isn't loaded: %s\n", fname) ;
    RETURN(0) ;
  }

  nim = populate_nifti_image(dset,options) ;
  if( !nim ) RETURN(0) ;   /* catch failure    6 Apr 2006 [rickr] */

  /*-- construct filename --*/

  nim->fname = malloc( strlen(fname)+16 ) ;
  nim->iname = malloc( strlen(fname)+16 ) ;
  strcpy(nim->fname,fname) ;
  strcpy(nim->iname,fname) ;

  /* 11 Oct 2005: Allow for .hdr/.img file outputs -- RWCox */

  cpt = nifti_find_file_extension( nim->iname ) ;
  if( cpt != NULL && strcmp(cpt,".hdr") == 0 ){
    nim->nifti_type = 2 ;   /* indicate 2 file output    */
    memcpy(cpt,".img",4) ;  /* convert .hdr name to .img */
  }

  /*-- construct nifti_brick_list of pointers to data briks */

  if( options.debug_level > 2 ) nifti_image_infodump(nim) ;
  nbl.bricks = (void **) malloc ( DSET_NVALS(dset) * sizeof(void*) ) ;
  nbl.nbricks = DSET_NVALS(dset) ;
  nbl.bsize = DSET_BRICK_BYTES(dset,0) ;
  for (ii = 0 ; ii < DSET_NVALS(dset) ; ii++ )
    nbl.bricks[ii] = DSET_ARRAY(dset,ii) ;

  /*-- 13 Mar 2006: check disk space --*/

  { FILE *fp = fopen(nim->fname,"ab") ;
    int mm   = THD_freemegabytes(nim->fname) ;
    int rr   = (int)(dset->dblk->total_bytes/(1024*1024)) ;
    if( fp != NULL ) fclose(fp) ;
    if( mm >= 0 && mm <= rr )
      WARNING_message("Disk space: writing dataset %s (%d MB),"
                      " but only %d free MB on disk"                   ,
              nim->fname , rr , mm ) ;
  }

  /*-- use handy-dandy library function to write out data */

  nifti_set_afni_extension( dset , nim ) ;  /* 09 May 2005 - RWCox */

  nifti_image_write_bricks (nim, &nbl ) ;
  RETURN(1) ;
}

/*******************************************************************/

nifti_image * populate_nifti_image(THD_3dim_dataset *dset, niftiwr_opts_t options)
{
  int nparam, type0 , ii , jj, val;
  int nif_x_axnum, nif_y_axnum, nif_z_axnum ;
  int slast, sfirst ;
  int pattern, tlen ;
  nifti_image *nim ;
  char axcode[3], axsign[3] ;
  float axstep[3] , axstart[3] ;
  int   axnum[3] ;
  float fac0 ;
  float dumqx, dumqy, dumqz, dumdx, dumdy, dumdz ;
  float *tlist;

ENTRY("populate_nifti_image") ;
  /*-- create nifti_image structure --*/

  nim = (nifti_image *) calloc( 1 , sizeof(nifti_image) ) ;
  if( !nim ) {
    fprintf(stderr, "** ERROR: failed to allocate nifti image\n");
    RETURN(NULL) ;
  }

  /*-- calculate and set ndim and intents --*/
  /* cases to allow for:
              1. 3d+time dataset
              2. 3d func bucket
                 (not going to handle this currently, may extend later)
                 -- RWC: modified so that 'u' dimension is for buckets
              3. 3d single brick
              4. 2d and 1d spatial
                 (if 2d or 1d spatial + time, treat as #1)
              5. Don't know of any vectors in AFNI, so won't do those either
              6. Single 3d (or 2d or 1d) functional brik  */

  if (dset->dblk->nvals > 1) {
    STATUS("4D dataset") ;
    nim->ndim = (dset->taxis != NULL) ? 4 : 5 ;  /* RWC: bucket stored as 5th dimen */

    /*-- check sub-bricks for uniformity in type and scale --*/

    type0 = DSET_BRICK_TYPE(dset,0) ;
    fac0  = DSET_BRICK_FACTOR(dset,0) ;

    for( ii=1 ; ii < DSET_NVALS(dset) ; ii++ ){
      if( DSET_BRICK_TYPE(dset,ii) != type0){
        fprintf(stderr,
        "** ERROR: CANNOT WRITE NIfTI FILE; BRICK DATA TYPES NOT CONSISTENT\n") ;
        RETURN(NULL);
      } else if( DSET_BRICK_FACTOR(dset,ii) != fac0) {
        fprintf(stderr,
        "** ERROR: CANNOT WRITE NIfTI FILE; BRICK FACTORS NOT CONSISTENT\n") ;
        fprintf(stderr,
        "   (consider transforming to a float dataset before performing\n"
        "    this operation, or consider '3dAFNItoNIFTI -float')\n"
        ) ;
        RETURN(NULL);
      }
    }
  } else {  /* we only have one brick */
    STATUS("3D dataset") ;
    if( options.debug_level > 1 ) fprintf(stderr,"ONLY ONE BRICK!!!\n") ;
    type0 = DSET_BRICK_TYPE(dset,0);
    fac0  = DSET_BRICK_FACTOR(dset,0) ;
    if (ISFUNC(dset)) {
      STATUS("functional dataset") ;
      if( options.debug_level > 1 )
        fprintf(stderr,"ONLY ONE BRICK, AND IT'S FUNCTIONAL!!!\n") ;
      nim->intent_code = DSET_BRICK_STATCODE(dset,0);
      if (nim->intent_code < 0) nim->intent_code = dset->func_type ;
      if (nim->intent_code < 0) nim->intent_code = NIFTI_INTENT_NONE ;
      if( options.debug_level > 1 )
        fprintf(stderr,"ONLY ONE BRICK, AND ITS FUNCTIONAL STAT CODE IS %d !!!\n",nim->intent_code) ;
      if(PRINT_TRACING){
        char str[256]; sprintf(str,"intent_code = %d",nim->intent_code);STATUS(str);
      }
      if (nim->intent_code > -1) {
        nparam = FUNC_need_stat_aux[nim->intent_code];
        if (nparam >= 1) nim->intent_p1 = DSET_BRICK_STATPAR(dset,0,1);
        if (nparam >= 2) nim->intent_p2 = DSET_BRICK_STATPAR(dset,0,2);
        if (nparam == 3) nim->intent_p3 = DSET_BRICK_STATPAR(dset,0,3);
      }
    }
    if (dset->daxes->nzz > 1) {
      nim->ndim = 3 ;
    } else if (dset->daxes->nyy > 1) {
      nim->ndim = 2 ;
    } else {
      nim->ndim = 1;
    }
  }


  /*-- set datatype, size, etc. --*/

  STATUS("set datatype") ;
  switch(type0) {
    case MRI_byte:
      nim->datatype = DT_UNSIGNED_CHAR;
      nim->nbyper = 1 ;
      break;
    case MRI_short:
      nim->datatype = DT_SIGNED_SHORT;
      nim->nbyper = 2 ;
      break;
    case MRI_int:
      nim->datatype = DT_SIGNED_INT;
      nim->nbyper = 4 ;
      break;
    case MRI_float:
      nim->datatype = DT_FLOAT;
      nim->nbyper = 4 ;
      break;
    case MRI_double:
      nim->datatype = DT_DOUBLE;
      nim->nbyper = 8 ;
      break;
    case MRI_complex:
      nim->datatype = DT_COMPLEX;
      nim->nbyper = 8 ;
      break;
    case MRI_rgb:
      nim->datatype = DT_RGB24;
      nim->nbyper = 3 ;
      break;
    case MRI_rgba:
      fprintf(stderr,
               "** ERROR: Can't write NIfTI file since dataset is RGBA: %s\n",
               options.infile_name) ;
      RETURN(NULL) ;
      break;
    default:
      fprintf(stderr,
               "** ERROR: Can't write NIfTI file since datatype is unknown: %s\n",
               options.infile_name) ;
      RETURN(NULL) ;
      break;
  }

  /*-- scaling --*/

  nim->scl_slope = fac0 ;
  nim->scl_inter = 0 ;

  /*-- spatial transforms --*/

  STATUS("set orientation") ;

  axcode[0] = ORIENT_xyz[ dset->daxes->xxorient ] ; axnum[0] = dset->daxes->nxx ;
  axcode[1] = ORIENT_xyz[ dset->daxes->yyorient ] ; axnum[1] = dset->daxes->nyy ;
  axcode[2] = ORIENT_xyz[ dset->daxes->zzorient ] ; axnum[2] = dset->daxes->nzz ;

  axsign[0] = ORIENT_sign[ dset->daxes->xxorient ] ;
  axsign[1] = ORIENT_sign[ dset->daxes->yyorient ] ;
  axsign[2] = ORIENT_sign[ dset->daxes->zzorient ] ;

  axstep[0] = dset->daxes->xxdel ; axstart[0] = dset->daxes->xxorg ;
  axstep[1] = dset->daxes->yydel ; axstart[1] = dset->daxes->yyorg ;
  axstep[2] = dset->daxes->zzdel ; axstart[2] = dset->daxes->zzorg ;

  for (ii = 0 ; ii < 3 ; ii++ ) {
    if (axcode[ii] == 'x') {
      nif_x_axnum = ii ;
    } else if (axcode[ii] == 'y') {
      nif_y_axnum = ii ;
    } else nif_z_axnum = ii ;
  }

  nim->qto_xyz.m[0][0] = nim->qto_xyz.m[0][1] = nim->qto_xyz.m[0][2] =
  nim->qto_xyz.m[1][0] = nim->qto_xyz.m[1][1] = nim->qto_xyz.m[1][2] =
  nim->qto_xyz.m[2][0] = nim->qto_xyz.m[2][1] = nim->qto_xyz.m[2][2] = 0.0 ;

  /*-- set voxel and time deltas and units --*/

  nim->dx = nim->pixdim[1] = fabs ( axstep[0] ) ;
  nim->dy = nim->pixdim[2] = fabs ( axstep[1] ) ;
  nim->dz = nim->pixdim[3] = fabs ( axstep[2] ) ;

  nim->du = nim->pixdim[5] = 0 ;
  nim->dv = nim->pixdim[6] = 0 ;
  nim->dw = nim->pixdim[7] = 0 ;

#if 0
  val = (axsign[nif_x_axnum] == '+')  ? -1 : 1 ;
  nim->qto_xyz.m[0][nif_x_axnum] = val  * nim->pixdim[nif_x_axnum + 1];
  val = (axsign[nif_y_axnum] == '+')  ? -1 : 1 ;
  nim->qto_xyz.m[1][nif_y_axnum] = val  * nim->pixdim[nif_y_axnum + 1];
  val = (axsign[nif_x_axnum] == '-')  ? -1 : 1 ;
  nim->qto_xyz.m[2][nif_z_axnum] = val  * nim->pixdim[nif_z_axnum + 1];
#else
  nim->qto_xyz.m[0][nif_x_axnum] = - axstep[nif_x_axnum];
  nim->qto_xyz.m[1][nif_y_axnum] = - axstep[nif_y_axnum];
  nim->qto_xyz.m[2][nif_z_axnum] =   axstep[nif_z_axnum];
#endif

  /* nifti origin stuff */

#if 0
  nim->qoffset_x =  axstart[nif_x_axnum] ;
  if (axsign[nif_x_axnum] == '+') nim->qoffset_x = - nim->qoffset_x ;
  nim->qoffset_y =  axstart[nif_y_axnum];
  if (axsign[nif_y_axnum] == '+') nim->qoffset_y = - nim->qoffset_y ;
  nim->qoffset_z =  axstart[nif_z_axnum];
  if (axsign[nif_z_axnum] == '-') nim->qoffset_z = - nim->qoffset_z ;
#endif

  nim->qoffset_x =  -axstart[nif_x_axnum] ;
  nim->qoffset_y =  -axstart[nif_y_axnum];
  nim->qoffset_z =  axstart[nif_z_axnum];

#if 0
  nim->qoffset_x =  -axstart[0] ;
  nim->qoffset_y =  -axstart[1];
  nim->qoffset_z =  axstart[2];
#endif

  nim->qto_xyz.m[0][3] = nim->qoffset_x ;
  nim->qto_xyz.m[1][3] = nim->qoffset_y ;
  nim->qto_xyz.m[2][3] = nim->qoffset_z ;


  /*-- from the same above info, set the sform matrix to equal the qform --*/
  /* KRH 7/6/05 - using sform to duplicate qform for
                           interoperability with FSL                       */
  /* update with oblique transformation if available DRG 24 May 2007 */
  /* check for valid transformation matrix */
  if(!ISVALID_MAT44(dset->daxes->ijk_to_dicom_real)) {
     nim->sto_xyz = nim->qto_xyz; /* copy qform to sform */
  }
  else {
      /* fill in sform with AFNI daxes transformation matrix */
      nim->sto_xyz = dset->daxes->ijk_to_dicom_real;
     /* negate first two rows of sform for NIFTI - LPI standard versus
                                            AFNI RAI "DICOM" standard */
     for( ii = 0; ii < 2; ii++) {
	for (jj = 0 ; jj < 4; jj++) {
            nim->sto_xyz.m[ii][jj] = -(nim->sto_xyz.m[ii][jj]);
	}
     }
     /* update qform too with struct copy from sform*/
     nim->qto_xyz= nim->sto_xyz ;

  }

  /*-- from the above info, calculate the quaternion qform --*/

  STATUS("set quaternion") ;

  nifti_mat44_to_quatern( nim->qto_xyz , &nim->quatern_b, &nim->quatern_c, &nim->quatern_d,
                    &dumqx, &dumqy, &dumqz, &dumdx, &dumdy, &dumdz, &nim->qfac ) ;


  /*-- verify dummy quaternion parameters --*/

  if( options.debug_level > 2 )
    fprintf(stderr,"++ Quaternion check:\n"
          "%f , %f\n %f , %f\n %f , %f\n %f , %f\n %f , %f\n %f , %f\n; %f\n",
           nim->qoffset_x, dumqx , nim->qoffset_y, dumqy , nim->qoffset_z, dumqz ,
           nim->dx, dumdx , nim->dy, dumdy , nim->dz, dumdz, nim->qfac ) ;

  /*-- calculate inverse qform            --*/

  nim->qto_ijk = nifti_mat44_inverse( nim->qto_xyz ) ;

  /*-- set dimensions of grid array --*/

  nim->nt = nim->nu = nim->nv = nim->nw = 1 ;
  nim->nx = axnum[0] ;
  nim->ny = axnum[1] ;
  nim->nz = axnum[2] ;

  if (dset->taxis == NULL) {
    nim->nu = DSET_NVALS(dset) ;   /* RWC: bucket is 5th dimension */
  } else {
    nim->nt = DSET_NUM_TIMES(dset) ;  /* time is 4th dimension */
  }

  if ( nim->nt > 1){
    float TR = dset->taxis->ttdel ;
    if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) TR *= 0.001; /* 10 May 2005 */
    nim->dt = nim->pixdim[4] = TR ;
  }

  nim->dim[0] = nim->ndim;
  nim->dim[1] = nim->nx;
  nim->dim[2] = nim->ny;
  nim->dim[3] = nim->nz;
  nim->dim[4] = nim->nt;  /* RWC: at most one of nt and nu is > 1 */
  nim->dim[5] = nim->nu;
  nim->dim[6] = nim->nv;
  nim->dim[7] = nim->nw;

  nim->nvox = nim->nx * nim->ny * nim->nz * nim->nt
                                * nim->nu * nim->nv * nim->nw ;

  /*-- slice timing --*/

  nim->freq_dim = nim->phase_dim = 0 ;
  if (dset->taxis != NULL) {  /* if time axis exists */
    nim->slice_dim = 3 ;
    nim->slice_duration = 0 ;
    nim->slice_start = 0 ;
    nim->slice_end = nim->nz - 1;
    nim->toffset =  DSET_TIMEORIGIN(dset);

    /*-- this bit assumes that afni slice timing offsets  *
     *-- are created starting from zero and including all *
     *-- slices initially.  They may later be modified by *
     *-- zero padding at either end.  No other            *
     *-- modifications are intentionally accepted right now. */

    if (DSET_NUM_TTOFF(dset) > 0 ) { /* if time offset exists */

      /*-- Find first and last non-zero element */
#define MYEPSILON 0.00001
#define MYFPEQ(a, b) (fabs((a) - (b)) < MYEPSILON)

      tlist = dset->taxis->toff_sl;
      for (ii = 0 ; ii < nim->nz ; ii++ ) {
        if (!MYFPEQ(tlist[ii],0.0)) break ;
      }
      sfirst = ii ;
      for (ii = nim->nz - 1 ; ii >= sfirst ; ii-- ) {
        if (!MYFPEQ(tlist[ii],0.0)) break ;
      }
      slast = ii ;

      /* pattern check re-written to deal with including zeros */
      /* on either end                     14 Jun 2006 [rickr] */

      pattern = NIFTI_SLICE_UNKNOWN;

      /* do we have all zeros? */
      if( sfirst == slast && MYFPEQ(tlist[sfirst],0.0) ) {
         nim->slice_duration = 0.0;
      } else { /* see if there is a known pattern in the list */
         tlen = slast-sfirst+2;

         /* try including leading adjacent zero in the pattern, first */
         if( sfirst > 0 ) {
            pattern = get_slice_timing_pattern(tlist+sfirst-1, tlen,
                                            &nim->slice_duration);
            if( pattern != NIFTI_SLICE_UNKNOWN ) sfirst--;
         }

         /* try including trailing adjacent zero in the pattern, next */
         if( pattern == NIFTI_SLICE_UNKNOWN && slast < nim->nz-1 ) {
            pattern = get_slice_timing_pattern(tlist+sfirst, tlen,
                                               &nim->slice_duration);
            if( pattern != NIFTI_SLICE_UNKNOWN ) slast++;
         }

         /* if no pattern yet, try list without zeros */
         if( pattern == NIFTI_SLICE_UNKNOWN )
            pattern = get_slice_timing_pattern(tlist+sfirst, tlen-1,
                                               &nim->slice_duration);

         if( pattern == NIFTI_SLICE_UNKNOWN ) {
            nim->slice_code = pattern ;
            nim->slice_start = 0 ;
            nim->slice_end = 0 ;
            nim->slice_duration = 0.0 ;
         } else {
            nim->slice_start = sfirst ;
            nim->slice_end = slast ;
            nim->slice_code = pattern;
         }

         if( options.debug_level > 1)
            fprintf(stderr,"+d timing pattern '%s', slice %d to %d, stime %f\n",
               nifti_slice_string(pattern), sfirst, slast, nim->slice_duration);
      }
    }

    nim->time_units = NIFTI_UNITS_SEC ;

  } else { /* if time axis not exists */
    nim->slice_dim = 0 ;
    nim->time_units = NIFTI_UNITS_UNKNOWN ;
  }

  /*-- byte order --*/

  nim->byteorder = nifti_short_order() ;

  /* KRH 7/25/05 modified to note talairach view into NIfTI file */

  if ( dset->view_type == VIEW_TALAIRACH_TYPE ) {
    nim->qform_code = NIFTI_XFORM_TALAIRACH ;
  } else {
    nim->qform_code = NIFTI_XFORM_SCANNER_ANAT ;
  }
  nim->sform_code = nim->qform_code ; /* KRH 7/6/05 - using */
           /* sform to duplicate qform for interoperability with FSL */


  /*-- odds and ends that are constant for AFNI files --*/
  nim->cal_min = nim->cal_max = 0 ;
  nim->nifti_type = 1 ;
  nim->xyz_units = NIFTI_UNITS_MM ;
  nim->num_ext = 0;
  nim->ext_list = NULL ;
  nim->iname_offset = 352 ; /* until extensions are added */
  nim->data = NULL ;

  RETURN(nim) ;
}

/*-------------------------------------------------------------------*/
/*! List of dataset attributes NOT to save in a NIfTI-1.1 file. -----*/

static char *badlist[] = {
     "IDCODE_STRING"      ,   /* this goes in the NI_group header */
     "DATASET_RANK"       ,
     "DATASET_DIMENSIONS" ,
     "TYPESTRING"         ,
     "SCENE_DATA"         ,
     "ORIENT_SPECIFIC"    ,
     "ORIGIN"             ,
     "DELTA"              ,
     "TAXIS_NUMS"         ,
     "TAXIS_FLOATS"       ,
     "TAXIS_OFFSETS"      ,
     "BYTEORDER_STRING"   ,
     "BRICK_TYPES"        ,
     "BRICK_FLOAT_FACS"   ,
     "STAT_AUX"           ,
     "LABEL_1"            ,
     "LABEL_2"            ,
     "DATASET_NAME"       ,
 NULL } ;

/*-------------------------------------------------------------------*/
/*! Create the AFNI extension string for a NIfTI-1.1 file, and insert
    this metadata into the nifti_image struct for output to disk.
    - If something bad happens, fails silently
    - 09 May 2005 - RWCox
---------------------------------------------------------------------*/

void nifti_set_afni_extension( THD_3dim_dataset *dset , nifti_image *nim )
{
   NI_group      *ngr ;
   NI_element    *nel ;
   NI_stream      ns  ;
   char *rhs , buf[128] ;
   int ii,bb , npart,*bpart ;

   if( nim == NULL                     ) return ;  /* stupid or evil caller */
   if( AFNI_yesenv("AFNI_NIFTI_NOEXT") ) return ;  /* not allowed */

   /** write all dataset 'attributes' into a NIML group */

   ngr = THD_nimlize_dsetatr( dset ) ;
   if( ngr == NULL ) return ;            /* bad */
   NI_rename_group( ngr , "AFNI_attributes" ) ;

   /* 12 May 2005: add a signature to check the file on input to AFNI */

   sprintf(buf,"%d,%d,%d,%d,%d,%d" ,
           nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->datatype ) ;
   NI_set_attribute( ngr , "NIfTI_nums" , buf ) ;

   /** now, scan attribute elements in the group, and mark some
       of them as being useless or redundant in the NIfTI world **/

   npart = ngr->part_num ;
   bpart = (int *)calloc(sizeof(int),npart) ;
   for( ii=0 ; ii < npart ; ii++ ){
     if( ngr->part_typ[ii] != NI_ELEMENT_TYPE ) continue ;
     nel = (NI_element *) ngr->part[ii] ;
     if( strcmp(nel->name,"AFNI_atr") != 0 )    continue ;
     rhs = NI_get_attribute( nel , "AFNI_name" ) ;
     if( rhs == NULL )                          continue ;

     for( bb=0 ; badlist[bb] != NULL ; bb++ )
       if( strcmp(rhs,badlist[bb]) == 0 ){ bpart[ii] = 1; break; }
   }

   /** remove marked attributes from the NIML group **/

   for( ii=npart-1 ; ii >= 0 ; ii-- ){
     if( bpart[ii] )
       NI_remove_from_group( ngr , ngr->part[ii] ) ;
   }
   free((void *)bpart) ;  /* done with this */
   if( ngr->part_num <= 0 ){ NI_free_element(ngr); return; }

   /** format into a character string to put in the NIfTI-1.1 extension **/

   ns = NI_stream_open( "str:" , "w" ) ;
   NI_stream_writestring( ns , "<?xml version='1.0' ?>\n" ) ;
   NI_write_element( ns , ngr , NI_TEXT_MODE ) ;
   rhs = NI_stream_getbuf( ns ) ;

   /** write contents of the string into the nifti_image struct **/

   nifti_add_extension( nim , rhs , strlen(rhs)+1 , NIFTI_ECODE_AFNI ) ;

   NI_stream_close(ns) ;   /* frees the string buffer, too */
   NI_free_element(ngr) ;  /* done with this trashola */
   return ;
}


/*! given a list of floats, detect any slice timing pattern
 *                                                14 Jun 2006 [rickr]
 *
 *    - if a pattern is found and delta is set, return the time delta
 *    - return one of:
 *        NIFTI_SLICE_UNKNOWN,
 *        NIFTI_SLICE_SEQ_INC,  NIFTI_SLICE_SEQ_DEC,
 *        NIFTI_SLICE_ALT_INC,  NIFTI_SLICE_ALT_DEC,
 *        NIFTI_SLICE_ALT_INC2, NIFTI_SLICE_ALT_DEC2,
 */
static int get_slice_timing_pattern( float * times, int len, float * delta )
{
   float * flist, diff;
   int   * ilist;
   int     c, index, pattern;

ENTRY("get_slice_timing_pattern");

   if( delta ) *delta = 0.0;  /* init, in case of early return */
   if( ! times || len < 2 ) RETURN(NIFTI_SLICE_UNKNOWN);

   /* if the length is very short, deal with it separately */
   if( len == 2 ){
      if( delta ) *delta = fabs(times[1]-times[0]);
      if( times[1] > times[0] ) RETURN(NIFTI_SLICE_SEQ_INC);
      else                      RETURN(NIFTI_SLICE_SEQ_DEC);
   }

   /*** sort the list, and look for a linear pattern ***/

   /* duplicate list */
   flist = (float *)malloc(len * sizeof(float));
   ilist = (int   *)malloc(len * sizeof(int));
   if(!flist || !ilist) {   /* yeah, lazy... */
      ERROR_message(" GSTP: cannot dupe timing list\n");
      RETURN(NIFTI_SLICE_UNKNOWN);
   }
   memcpy(flist, times, len*sizeof(float));
   for(c = 0; c < len; c++) ilist[c] = c;  /* init ilist with current indices */

   /* sort flist, with ilist returning original indices */
   qsort_floatint(len, flist, ilist);

   /* and check for a fixed difference */
   diff = flist[1] - flist[0];
   pattern = 1;
   for( c = 1; c < len-1; c++ )
     if( !MYFPEQ(diff, (flist[c+1]-flist[c])) ) { pattern = 0; break; }

   /* if no pattern, just return failure */
   if( !pattern ) {
      free(flist);  free(ilist);  RETURN(NIFTI_SLICE_UNKNOWN);
   }

   /* we have linear offsets, now see if the slices match a known pattern */
   /* repeatedly: init to a pattern, and see if it fails */
   
   /* SEQ_INC  (0,1,2,3...,l-1) */
   pattern = NIFTI_SLICE_SEQ_INC;  index = 0;
   for( c = 0; c < len; c++ ) {
      if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN;  break; }
      index++;
   }

   if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (l-1,l-2,...2,1,0) */
       pattern = NIFTI_SLICE_SEQ_DEC;  index = len-1;
       for( c = 0; c < len; c++ ) {
          if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN;  break; }
          index--;
       }
   }

   if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (0,2,4,6,...,1,3,5,...) */
       pattern = NIFTI_SLICE_ALT_INC;  index = 0;
       for( c = 0; c < len; c++ ) {
          if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN;  break; }
          index += 2;  if( index >= len ) index = 1;  /* so no parity issue */
       }
   }

   if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (l-1,l-3,...1/0,l-2,l-4,...,0/1) */
       pattern = NIFTI_SLICE_ALT_DEC;  index = len-1;
       for( c = 0; c < len; c++ ) {
          if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN;  break; }
          index -= 2;  if( index < 0 ) index = len-2;
       }
   }

   if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (1,3,5,...,0,2,4...) */
       pattern = NIFTI_SLICE_ALT_INC2;  index = 1;
       for( c = 0; c < len; c++ ) {
          if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN;  break; }
          index += 2;  if( index >= len ) index = 0;
       }
   }

   if( pattern == NIFTI_SLICE_UNKNOWN ) { /* (l-2,l-4,...4,2,0,l-1,...,5,3,1) */
       pattern = NIFTI_SLICE_ALT_DEC2;  index = len-2;
       for( c = 0; c < len; c++ ) {
          if( ilist[c] != index ) { pattern = NIFTI_SLICE_UNKNOWN;  break; }
          index -= 2;  if( index < 0 ) index = len-1;
       }
   }

   if( delta && pattern != NIFTI_SLICE_UNKNOWN ) *delta = diff;

   /** done, whatever the case may be **/
   free(flist);  free(ilist);
   RETURN(pattern);
}