/*****************************************************************************
   Major portions of this software are copyrighted by the Medical College
   of Wisconsin, 1994-2000, and are released under the Gnu General Public
   License, Version 2.  See the file README.Copyright for details.
******************************************************************************/
   
#include "afni.h"

#ifndef ALLOW_PLUGINS
#  error "Plugins not properly set up -- see machdep.h"
#endif

#define TEST_VOXEL 6177
#define TEST_TIME 0
#define RMB_DEBUG 0

/***********************************************************************
  Plugin that averages epochs from single trial data
************************************************************************/

/*--------------------- string to 'help' the user --------------------*/

static char helpstring[] =
  "Purpose: Averaging epochs of single trial data\n"
  "\n"
  "Input items to this plugin are:\n"
  "   Datasets:   Input  = 3D+time dataset to process\n"
  "                      = reference time-series\n"
  "               Output = Prefix for new dataset\n"
  "   Additional Parameters\n"
  "               delta     = shift timeseries by delta before splitting and averaging\n"
  "               method    = type of statistic to calculate\n"
  "               maxlength = maximum avg ts length\n"
  "               no1?      = images w/ only one img in avg ignored\n"
  "Author -- RM Birn"
;

/*------------- strings for output format -------------*/

static char * yes_no_strings[] = { "No" , "Yes" } ;
static char * method_strings[] = { "Mean" , "Sigma" } ;

#define _STAVG_NUM_METHODS (sizeof(method_strings)/sizeof(char *))
#define _STAVG_METH_MEAN 0
#define _STAVG_METH_SIGMA 1

/*--------------- prototypes for internal routines ---------------*/

char * STAVG_main( PLUGIN_interface * ) ;  /* the entry point */

float ** avg_epochs( THD_3dim_dataset * dset, float * ref, 
                    int user_maxlength, int no1, int meth,
		    PLUGIN_interface *plint );

MRI_IMARR * dset_to_mri(THD_3dim_dataset * dset);


/*---------------------------- global data ---------------------------*/

static PLUGIN_interface * global_plint = NULL ;
int M_maxlength;

/***********************************************************************
   Set up the interface to the user:
    1) Create a new interface using "PLUTO_new_interface";

    2) For each line of inputs, create the line with "PLUTO_add_option"
         (this line of inputs can be optional or mandatory);

    3) For each item on the line, create the item with
        "PLUTO_add_dataset" for a dataset chooser,
        "PLUTO_add_string"  for a string chooser,
        "PLUTO_add_number"  for a number chooser.
************************************************************************/


DEFINE_PLUGIN_PROTOTYPE

PLUGIN_interface * PLUGIN_init( int ncall )
{
   PLUGIN_interface * plint ;     /* will be the output of this routine */

   if( ncall > 0 ) return NULL ;  /* one interfaces */

   CHECK_IF_ALLOWED("SINGLETRIALAVG","SingleTrial Avg") ;  /* 30 Sep 2016 */

   /*---------------- set titles and call point ----------------*/

   plint = PLUTO_new_interface( "SingleTrial Avg" ,
                                "Averaging of epochs in Single Trial data" ,
                                helpstring ,
                                PLUGIN_CALL_VIA_MENU ,
                                (cptr_func *)STAVG_main  ) ;

   PLUTO_add_hint( plint , "Averaging of epochs in Single Trial data" ) ;

   global_plint = plint ;  /* make global copy */

   PLUTO_set_sequence( plint , "z:Birn" ) ;

   /*--------- 1st line ---------*/

   PLUTO_add_option( plint ,
                     "Datasets" ,  /* label at left of input line */
                     "Datasets" ,  /* tag to return to plugin */
                     TRUE          /* is this mandatory? */
                   ) ;

   PLUTO_add_dataset(  plint ,
                       "Input" ,          /* label next to button   */
                       ANAT_ALL_MASK ,    /* take any anat datasets */
                       FUNC_FIM_MASK ,    /* only allow fim funcs   */
                       DIMEN_4D_MASK |    /* need 3D+time datasets  */
                       BRICK_ALLREAL_MASK /* need real-valued datasets */
                    ) ;
   PLUTO_add_hint( plint , "Input 3d+t dataset" ) ;

   PLUTO_add_string( plint ,
                     "Output" ,  /* label next to textfield */
                     0,NULL ,    /* no fixed strings to choose among */
                     19          /* 19 spaces for typing in value */
                   ) ;
   PLUTO_add_hint( plint , "Name of output dataset" ) ;

   /*---------- 2nd line --------*/

   PLUTO_add_option( plint ,
                     "Timing" ,
                     "Timing" ,
                     TRUE
                   ) ;


   PLUTO_add_timeseries(plint, "Stim. Timing");
   PLUTO_add_hint( plint , "Stimulus Timing (0 = no task, 1 = task)" ) ;

   PLUTO_add_number( plint ,
                     "delta" ,   
                     -1000 ,    
                     1000 ,  
                     0 ,    
                     0 ,   
                     TRUE
                   ) ;
   PLUTO_add_hint( plint , "Shift data timecourse by delta before splitting and averaging" ) ;

   /*---------- 3rd line: computation ----------*/

   PLUTO_add_option( plint ,
                     "Compute" ,  /* label at left of input line */
                     "Compute" ,  /* tag to return to plugin */
                     TRUE         /* is this mandatory? */
                   ) ;

   PLUTO_add_string( plint ,
                     "Method" ,           /* label next to chooser button */
                     _STAVG_NUM_METHODS,  /* number of strings in list */
                     method_strings ,     /* list of strings to choose among */
                     _STAVG_METH_MEAN     /* index of default string */
                   ) ;

   PLUTO_add_hint( plint , "Choose statistic to compute" ) ;

   /*---------- 4th line --------*/

   PLUTO_add_option( plint ,
                     "Parameters" ,  /* label at left of input line */
                     "Parameters" ,  /* tag to return to plugin */
                     FALSE            /* is this mandatory? */
                   ) ;

   PLUTO_add_number( plint ,
                     "maxlength" ,    /* label next to chooser */
                     0 ,         /* smallest possible value */
                     1000 ,        /* largest possible value */
                     0 ,         /* decimal shift (none in this case) */
                     15 ,         /* default value */
                     TRUE       /* allow user to edit value? */
                   ) ;
   PLUTO_add_hint( plint , "maximum # of timepoints of output dataset" ) ;

   PLUTO_add_string( plint ,
                     "no1?" ,               /* label next to chooser button */
                     2  ,               /* number of strings to choose among */
                     yes_no_strings ,  /* list of strings to choose among */
                     1                  /* index of default string */
                   ) ;

   PLUTO_add_hint( plint , "ignore timepoints where only one image is in average" ) ;


   /*--------- done with interface setup ---------*/

   return plint ;
}

/***************************************************************************
  Main routine for this plugin (will be called from AFNI).
  If the return string is not NULL, some error transpired, and
  AFNI will popup the return string in a message box.
****************************************************************************/

/*------------------ macros to return workspace at exit -------------------*/

#undef  FREEUP
#define FREEUP(x) if((x) != NULL){free((x)); (x)=NULL;}

#define FREE_WORKSPACE                              \
  do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
      FREEUP(fout) ; \
      FREEUP(fxar) ; \
    } while(0)

/*-------------------------------------------------------------------------*/

char * STAVG_main( PLUGIN_interface * plint )
{
   MCW_idcode * idc ;                          /* input dataset idcode */
   THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
   char * new_prefix , * str , * str2;         /* strings from user */
   int   meth;                                 /* chosen computation method */
   int   new_datum ,                           /* control parameters */
         old_datum , ntime ;

   int   te, ne, tinc, kim, nia;
   int   numepochs, minlength, maxlength, lastindex, navgpts;
   int   nvox , perc , new_units, old_units ;
   int   ii, ibot,itop , kk, jj; 
   int   no1, user_maxlength, delta;
   int   *pEpochLength, *pTimeIndex;
   int   nx, ny, nz, npix;
   float *pNumAvg;
   float old_dtime;

   MRI_IMAGE * stimim;
   
   MRI_IMARR *avgimar;

   byte   ** bptr  = NULL ;  /* one of these will be the array of */
   short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
   float  ** fptr  = NULL ;  /* (depending on input datum type) */

   float   * fxar  = NULL ;  /* array loaded from input dataset */
   float   * stimar = NULL ;
   float  ** fout  = NULL ;  /* will be array of output floats */

   float   * tar   = NULL ;  /* will be array of taper coefficients */

   float   * nstimar;

   /*--------------------------------------------------------------------*/
   /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/

   /*--------- go to first input line ---------*/

   PLUTO_next_option(plint) ;

   idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
   old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
   if( old_dset == NULL )
      return "*************************\n"
             "Cannot find Input Dataset\n"
             "*************************"  ;

   ntime = DSET_NUM_TIMES(old_dset) ;
   if( ntime < 2 )
      return "*****************************\n"
             "Dataset has only 1 time point\n"
             "*****************************"  ;

   ii = DSET_NVALS_PER_TIME(old_dset) ;
   if( ii > 1 )
      return "************************************\n"
             "Dataset has > 1 value per time point\n"
             "************************************"  ;
   
   old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
   new_datum = old_datum;
   old_dtime = DSET_TIMESTEP(old_dset);
   old_units = DSET_TIMEUNITS(old_dset);
   
   nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
   npix = old_dset->daxes->nxx * old_dset->daxes->nyy;
   nx = old_dset->daxes->nxx;


   new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
   if( ! PLUTO_prefix_ok(new_prefix) )      /* check if it is OK */
      return "************************\n"
             "Output Prefix is illegal\n"
             "************************"  ;

   /*--------- go to next input line ---------*/

   PLUTO_next_option(plint);

   stimim = PLUTO_get_timeseries(plint);
   if( stimim == NULL ) return "Please specify stimulus timing";

   if( stimim->nx < ntime ){
      return "**************************************\n"
             "Not enough pts in stimulus time-series\n"
             "**************************************";
   }

   stimar = MRI_FLOAT_PTR(stimim);


   delta = PLUTO_get_number(plint);

   if( abs(delta) > ntime ){
      return "************************\n"
             "Delta shift is too large\n"
             "************************";
   }
  
   /*initialize variables if not user specified */
   user_maxlength = ntime;
   no1 = 0;

   /*--------- go to next input line ---------*/

   PLUTO_next_option(plint);

   str  = PLUTO_get_string(plint) ;      /* get string item (the method) */
   meth = PLUTO_string_index( str ,      /* find it in list it is from */
                              _STAVG_NUM_METHODS ,
                              method_strings ) ;

   /*--------- see if the 4th option line is present --------*/

   str = PLUTO_get_optiontag( plint ) ;
   if( str != NULL ){
      user_maxlength = (int) PLUTO_get_number(plint) ;
      str2  = PLUTO_get_string(plint) ;      /* get string item (the method) */
      no1   = PLUTO_string_index( str2 ,      /* find it in list it is from */
                                 2 ,
                                 yes_no_strings) ;
   }
   

   /*------------------------------------------------------*/
   /*---------- At this point, the inputs are OK ----------*/

   PLUTO_popup_meter( plint ) ;  /* popup a progress meter */

   /*________________[ Main Code ]_________________________*/
  
   fout = avg_epochs( old_dset, stimar, user_maxlength, 1, meth, plint );
   if( fout == NULL ) return " \nError in avg_epochs() function!\n " ;
   
   if( RMB_DEBUG ) fprintf(stderr, "Done with avg_epochs\n");
   maxlength = M_maxlength;
   
   
   /*______________________________________________________*/

   
   new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */

   { char * his = PLUTO_commandstring(plint) ;
     tross_Copy_History( old_dset , new_dset ) ;
     tross_Append_History( new_dset , his ) ; free( his ) ;
   }
   
   /*-- edit some of its internal parameters --*/
   ii = EDIT_dset_items(
           new_dset ,
              ADN_prefix      , new_prefix ,           /* filename prefix */
              ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
              ADN_datum_all   , new_datum ,            /* atomic datum */
              ADN_nvals       , maxlength ,            /* # sub-bricks */
              ADN_ntt         , maxlength ,            /* # time points */
           /*   ADN_ttorg       , old_dtime ,  */              /* time origin */
           /*   ADN_ttdel       , old_dtime ,  */            /* time step */
           /*   ADN_ttdur       , old_dtime ,  */            /* time duration */
           /*   ADN_nsl         , 0 ,          */        /* z-axis time slicing */
           /*   ADN_tunits      , old_units ,  */        /* time units */
           ADN_none ) ;

   if( ii != 0 ){
      THD_delete_3dim_dataset( new_dset , False ) ;
      FREE_WORKSPACE ;
      return "***********************************\n"
             "Error while creating output dataset\n"
             "***********************************"  ;
   }


   /*------------------------------------------------------------*/
   /*------- The output is now in fout[kk][ii],
             for kk=0..maxlength-1 , ii=0..nvox-1.
             We must now put this into the output dataset -------*/

   switch( new_datum ){

      /*** output is floats is the simplest:
           we just have to attach the fout bricks to the dataset ***/

      case MRI_float:
         for( kk=0 ; kk < maxlength ; kk++ )
            EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
      break ;

      /*** output is shorts:
           we have to create a scaled sub-brick from fout ***/

      case MRI_short:{
         short * bout ;
         float fac ; 

         for( kk=0 ; kk < maxlength ; kk++ ){  /* loop over sub-bricks */

            /*-- get output sub-brick --*/
            bout = (short *) malloc( sizeof(short) * nvox ) ;
            if( bout == NULL ){
               fprintf(stderr,"\nFinal malloc error in plug_stavg!\n\a") ;
               return("Final malloc error in plug_stavg!"); ;
               /*  exit(1) ;*/
            }

            /*-- find scaling and then scale --*/
            /*fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;*/
            fac = 1.0;
            EDIT_coerce_scale_type( nvox,fac ,
                                    MRI_float,fout[kk] , MRI_short,bout ) ;
            free( fout[kk] ) ;  /* don't need this anymore */

            /*-- put output brick into dataset, and store scale factor --*/
            EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
         }
      }
      break ;

      /*** output is bytes (byte = unsigned char)
           we have to create a scaled sub-brick from fout ***/

      case MRI_byte:{
         byte * bout ;
         float fac ;

         for( kk=0 ; kk < maxlength ; kk++ ){  /* loop over sub-bricks */

            /*-- get output sub-brick --*/

            bout = (byte *) malloc( sizeof(byte) * nvox ) ;
            if( bout == NULL ){
               fprintf(stderr,"\nFinal malloc error in plug_stavg!\n\a") ;
               return("Final malloc error in plug_stavg!"); ;
	       /*               exit(1) ;*/
            }

            /*-- find scaling and then scale --*/

            fac = 1.0;
            EDIT_coerce_scale_type( nvox,fac ,
                                    MRI_float,fout[kk] , MRI_byte,bout ) ;

            free( fout[kk] ) ;  /* don't need this anymore */

            /*-- put output brick into dataset, and store scale factor --*/

            EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
         }

      }
      break ;

   } /* end of switch on output data type */

   /*-------------- Cleanup and go home ----------------*/

   PLUTO_set_meter( plint , 100 ) ;  /* set progress meter to 100% */

   PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;

   FREE_WORKSPACE ;
   return NULL ;  /* null string returned means all was OK */
}

/*---------------------------------------------------------------*/
float ** avg_epochs( THD_3dim_dataset * dset, float * ref, 
                    int user_maxlength, int no1, int meth,
		    PLUGIN_interface * plint )
/*---------------------------------------------------------------*/
{

   int     numepochs, lastindex;
   int     nvox, numims, nx, ny, nz;
   int     kim, ne, te, tinc, nia;
   int     ii, kk;
   int     maxlength, minlength;
   int     datum;
   float   fac;       /* scaling factor for current subbrick */
   double  scaledval; /* temp var for scaled value to be used repeatedly */
   float ** fxar;
   float ** outar;    /* output averaged time-series */
   float * sumsq = NULL; /* sum of squared voxel values */
   float * pNumAvg;  /* array for number of pts to avg at each time*/
   int   * pTimeIndex; /* array of time markers (1st img of each epoch) */
   int   * pEpochLength; /* array of epoch lengths */
   float ** tempar;
   MRI_IMARR *inimar;
   
   
   nx = dset->daxes->nxx;
   ny = dset->daxes->nyy;
   nz = dset->daxes->nzz;
   nvox = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
   numims = DSET_NUM_TIMES(dset);
   datum = DSET_BRICK_TYPE( dset , 0 ) ; /* get old dataset datum type */
   
   PLUTO_popup_meter(plint) ;
   
   DSET_load(dset);
   
   inimar =  dset_to_mri(dset);
   if( inimar == NULL ) return NULL ;

   fxar = (float **) malloc( sizeof( float *) * numims);
   if( datum == MRI_float){
      for( kk=0; kk<numims; kk++){
         fxar[kk] = MRI_FLOAT_PTR(IMAGE_IN_IMARR(inimar,kk));
      }
   }
   else{
      for( kk=0; kk<numims; kk++){
         fxar[kk] = MRI_FLOAT_PTR(mri_to_float(IMAGE_IN_IMARR(inimar,kk)));
      }
   }

   nia = 0;    /* number of images (timepoints) averaged  where num epochs > 1*/
   
   if( RMB_DEBUG ) fprintf(stderr, "Start stavg\n");
   
   /* determine number of epochs to average */
   if( RMB_DEBUG ) fprintf(stderr, "Determining number of epochs...");
   numepochs = 1;
   for( kim=0; kim < numims; kim++ ){
      if( ref[kim] > 0) numepochs++;
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");

   /* set array of epoch lengths */
   if( RMB_DEBUG ) fprintf(stderr, "Set array of epoch lengths...");
   pEpochLength = (int *)malloc(sizeof(int) * numepochs);
   for( ne=0; ne < numepochs; ne++) pEpochLength[ne] = 0;
   ne = 0;
   for( kim=0; kim < numims; kim++ ){
      if( ref[kim] > 0) ne++;
      pEpochLength[ne]++;
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");

   /* set array of time markers (1st img of each epoch) */
   if( RMB_DEBUG ) fprintf(stderr, "Set array of time markers...");
   pTimeIndex = (int *)malloc(sizeof(int) * (numepochs - 1));
   lastindex = 0;
   minlength = numims;
   maxlength = 0;
   for( ne=0; ne < (numepochs-1); ne++){
      pTimeIndex[ne] = lastindex + pEpochLength[ne];
      lastindex = pTimeIndex[ne];
      if(pEpochLength[ne+1] > maxlength) maxlength = pEpochLength[ne+1];
      if(pEpochLength[ne+1] < minlength) minlength = pEpochLength[ne+1];
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");

   if(maxlength > user_maxlength) maxlength = user_maxlength;

   if( RMB_DEBUG ) fprintf(stderr, "init...");
   pNumAvg = (float *) malloc( sizeof(float) * maxlength);
   outar = (float **) malloc( sizeof(float *) * maxlength);
   for( te=0; te < maxlength; te++){
      outar[te] = (float *) malloc( sizeof(float) * nvox);
      for( ii=0; ii<nvox; ii++){
         outar[te][ii] = 0.0;
      }
   }
   if (meth == _STAVG_METH_SIGMA) {
       sumsq = (float *) malloc( sizeof(float) * nvox);
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");

   if( RMB_DEBUG ) fprintf(stderr, "Start averaging...");
   for( te=0; te < maxlength; te++){
      pNumAvg[te] = 0.0;   

     switch(meth) {
     default: case _STAVG_METH_MEAN:

      for( ne=0; ne < (numepochs-1); ne++){
         tinc = pTimeIndex[ne] + te;
         if( te < pEpochLength[ne+1] ){
            fac = DSET_BRICK_FACTOR(dset, tinc);
	    if (fac == 0.0 || fac == 1.0) {
		for( ii=0; ii<nvox; ii++){
		    outar[te][ii] += fxar[tinc][ii];
		}
	    } else {
		for( ii=0; ii<nvox; ii++){
		    outar[te][ii] += fxar[tinc][ii] * fac;
		}
	    }
            pNumAvg[te]++;
         }
      }
      for( ii=0; ii<nvox; ii++){
	 outar[te][ii] = outar[te][ii]/pNumAvg[te];
      }
      break;

     case _STAVG_METH_SIGMA:

      /* sigma statistic is actually the unbiased standard deviation
	 calculated with this computational formula:
	   sqrt( (sum(x^2) - sum(x)^2 / n) / (n-1) )  */
      for( ii=0; ii<nvox; ii++){
	  sumsq[ii] = 0.0;
      }
      for( ne=0; ne < (numepochs-1); ne++){
         tinc = pTimeIndex[ne] + te;
         if( te < pEpochLength[ne+1] ){
            fac = DSET_BRICK_FACTOR(dset, tinc);
	    if (fac == 0.0 || fac == 1.0) {
		for( ii=0; ii<nvox; ii++){
		    outar[te][ii] += fxar[tinc][ii];
		    sumsq[ii] += fxar[tinc][ii] * fxar[tinc][ii];
		}
	    } else {
		for( ii=0; ii<nvox; ii++){
		    scaledval = fxar[tinc][ii] * fac;
		    outar[te][ii] += scaledval;
		    sumsq[ii] += scaledval * scaledval;
		}
	    }
            pNumAvg[te]++;
         }
      }
      for( ii=0; ii<nvox; ii++){
	  outar[te][ii] = sqrt( (sumsq[ii] -
				 outar[te][ii] * outar[te][ii] /
				 pNumAvg[te]) /
				(pNumAvg[te] - 1) );
      }
      break;

     }

      if( pNumAvg[te] > 1) nia ++;
      PLUTO_set_meter(plint, (100*(te+1))/maxlength ) ;
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");
   if ( no1 ){     /* ignore images with only one average */
      if( nia < maxlength) maxlength = nia;
   }
   
   M_maxlength = maxlength;
   
   
   if( RMB_DEBUG ) fprintf(stderr, "malloc output...");
   tempar = (float **) malloc(sizeof(float *) * maxlength);
   for( te=0 ; te < maxlength ; te++ ){
      tempar[te] = (float *) malloc( sizeof(float) * nvox);
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");
   
   if( RMB_DEBUG ) fprintf(stderr, "convert to output...");
   for( te=0; te < maxlength; te++){
      for( ii=0; ii<nvox; ii++){
         tempar[te][ii] = outar[te][ii];
      }
   }
   if( RMB_DEBUG ) fprintf(stderr, "done\n");

   /* toss arrays */
   if( RMB_DEBUG ) fprintf(stderr, "free mem...");
   FREE_IMARR( inimar );
   if (meth == _STAVG_METH_SIGMA) free(sumsq);
   free( outar );
   free( pNumAvg );
   free( pTimeIndex );
   free( pEpochLength );
   if( RMB_DEBUG ) fprintf(stderr, "done\n");
   DSET_unload(dset);
   
   return(tempar);
}


/*-------------------------------------------------------*/
MRI_IMARR * dset_to_mri(THD_3dim_dataset * dset)
/*--------------------------------------------------------*/
{

   int ii, kk, ntime, datum;
   int nvox, nx, ny, nz;
   int use_fac;
   
   MRI_IMARR * ims_in;
   MRI_IMAGE * im, *temp_im;
   

   byte   ** bptr  = NULL ;  /* one of these will be the array of */
   short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
   float  ** fptr  = NULL ;  /* (depending on input datum type) */
   
   float * fac  = NULL ;  /* array of brick scaling factors */
   
   float * fout;
   

   ntime = DSET_NUM_TIMES(dset) ;
   nx = dset->daxes->nxx;
   ny = dset->daxes->nyy;
   nz = dset->daxes->nzz;
   nvox = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
   datum = DSET_BRICK_TYPE( dset , 0 ) ; /* get dataset datum type */

   switch( datum ){  /* pointer type depends on input datum type */

      default:
         return NULL  ;

      /** create array of pointers into old dataset sub-bricks **/

      /*--------- input is bytes ----------*/
      /* voxel #i at time #k is bptr[k][i] */
      /* for i=0..nvox-1 and k=0..ntime-1.  */

      case MRI_byte:
         bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
         if( bptr == NULL ) return NULL ;
         for( kk=0 ; kk < ntime ; kk++ )
            bptr[kk] = (byte *) DSET_ARRAY(dset,kk) ;
      break ;

      /*--------- input is shorts ---------*/
      /* voxel #i at time #k is sptr[k][i] */
      /* for i=0..nvox-1 and k=0..ntime-1.  */

      case MRI_short:
         sptr = (short **) malloc( sizeof(short *) * ntime ) ;
         if( sptr == NULL ) return NULL ;
         for( kk=0 ; kk < ntime; kk++ )
            sptr[kk] = (short *) DSET_ARRAY(dset,kk) ;
      break ;

      /*--------- input is floats ---------*/
      /* voxel #i at time #k is fptr[k][i] */
      /* for i=0..nvox-1 and k=0..ntime-1.  */

      case MRI_float:
         fptr = (float **) malloc( sizeof(float *) * ntime) ;
         if( fptr == NULL ) return NULL ;
         for( kk=0 ; kk < ntime; kk++ )
            fptr[kk] = (float *) DSET_ARRAY(dset,kk) ;
      break ;

   } /* end of switch on input type */
   
   INIT_IMARR(ims_in) ;
   for( kk=0 ; kk < ntime ; kk++ ){
      im = mri_new_vol_empty( nx , ny , nz , datum ) ;
      ADDTO_IMARR(ims_in,im) ;
   }
   
   for( kk=0 ; kk < ntime ; kk++ ){
      im = IMARR_SUBIMAGE(ims_in,kk) ;
      
      switch( datum ){
         case MRI_byte:  mri_fix_data_pointer( bptr[kk], im ) ; break ;
         case MRI_short: mri_fix_data_pointer( sptr[kk], im ) ; break ;
         case MRI_float: mri_fix_data_pointer( fptr[kk], im ) ; break ;
      }
   }


   
   return(ims_in);
}


