00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "afni.h"
00008 
00009 #ifndef ALLOW_PLUGINS
00010 #  error "Plugins not properly set up -- see machdep.h"
00011 #endif
00012 
00013 #define TEST_VOXEL 6177
00014 #define TEST_TIME 0
00015 #define RMB_DEBUG 0
00016 
00017 /***********************************************************************
00018   Plugin that averages epochs from single trial data
00019 ************************************************************************/
00020 
00021 /*--------------------- string to 'help' the user --------------------*/
00022 
00023 static char helpstring[] =
00024   "Purpose: Averaging epochs of single trial data\n"
00025   "\n"
00026   "Input items to this plugin are:\n"
00027   "   Datasets:   Input  = 3D+time dataset to process\n"
00028   "                      = reference time-series\n"
00029   "               Output = Prefix for new dataset\n"
00030   "   Additional Parameters\n"
00031   "               delta     = shift timeseries by delta before splitting and averaging\n"
00032   "               method    = type of statistic to calculate\n"
00033   "               maxlength = maximum avg ts length\n"
00034   "               no1?      = images w/ only one img in avg ignored\n"
00035   "Author -- RM Birn"
00036 ;
00037 
00038 /*------------- strings for output format -------------*/
00039 
00040 static char * yes_no_strings[] = { "No" , "Yes" } ;
00041 static char * method_strings[] = { "Mean" , "Sigma" } ;
00042 
00043 #define _STAVG_NUM_METHODS (sizeof(method_strings)/sizeof(char *))
00044 #define _STAVG_METH_MEAN 0
00045 #define _STAVG_METH_SIGMA 1
00046 
00047 /*--------------- prototypes for internal routines ---------------*/
00048 
00049 char * STAVG_main( PLUGIN_interface * ) ;  /* the entry point */
00050 
00051 float ** avg_epochs( THD_3dim_dataset * dset, float * ref, 
00052                     int user_maxlength, int no1, int meth,
00053                     PLUGIN_interface *plint );
00054 
00055 MRI_IMARR * dset_to_mri(THD_3dim_dataset * dset);
00056 
00057 
00058 /*---------------------------- global data ---------------------------*/
00059 
00060 static PLUGIN_interface * global_plint = NULL ;
00061 int M_maxlength;
00062 
00063 /***********************************************************************
00064    Set up the interface to the user:
00065     1) Create a new interface using "PLUTO_new_interface";
00066 
00067     2) For each line of inputs, create the line with "PLUTO_add_option"
00068          (this line of inputs can be optional or mandatory);
00069 
00070     3) For each item on the line, create the item with
00071         "PLUTO_add_dataset" for a dataset chooser,
00072         "PLUTO_add_string"  for a string chooser,
00073         "PLUTO_add_number"  for a number chooser.
00074 ************************************************************************/
00075 
00076 
00077 DEFINE_PLUGIN_PROTOTYPE
00078 
00079 PLUGIN_interface * PLUGIN_init( int ncall )
00080 {
00081    PLUGIN_interface * plint ;     /* will be the output of this routine */
00082 
00083    if( ncall > 0 ) return NULL ;  /* one interfaces */
00084 
00085    /*---------------- set titles and call point ----------------*/
00086 
00087    plint = PLUTO_new_interface( "SingleTrial Avg" ,
00088                                 "Averaging of epochs in Single Trial data" ,
00089                                 helpstring ,
00090                                 PLUGIN_CALL_VIA_MENU , STAVG_main  ) ;
00091 
00092    PLUTO_add_hint( plint , "Averaging of epochs in Single Trial data" ) ;
00093 
00094    global_plint = plint ;  /* make global copy */
00095 
00096    PLUTO_set_sequence( plint , "z:Birn" ) ;
00097 
00098    /*--------- 1st line ---------*/
00099 
00100    PLUTO_add_option( plint ,
00101                      "Datasets" ,  /* label at left of input line */
00102                      "Datasets" ,  /* tag to return to plugin */
00103                      TRUE          /* is this mandatory? */
00104                    ) ;
00105 
00106    PLUTO_add_dataset(  plint ,
00107                        "Input" ,          /* label next to button   */
00108                        ANAT_ALL_MASK ,    /* take any anat datasets */
00109                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00110                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00111                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00112                     ) ;
00113    PLUTO_add_hint( plint , "Input 3d+t dataset" ) ;
00114 
00115    PLUTO_add_string( plint ,
00116                      "Output" ,  /* label next to textfield */
00117                      0,NULL ,    /* no fixed strings to choose among */
00118                      19          /* 19 spaces for typing in value */
00119                    ) ;
00120    PLUTO_add_hint( plint , "Name of output dataset" ) ;
00121 
00122    /*---------- 2nd line --------*/
00123 
00124    PLUTO_add_option( plint ,
00125                      "Timing" ,
00126                      "Timing" ,
00127                      TRUE
00128                    ) ;
00129 
00130 
00131    PLUTO_add_timeseries(plint, "Stim. Timing");
00132    PLUTO_add_hint( plint , "Stimulus Timing (0 = no task, 1 = task)" ) ;
00133 
00134    PLUTO_add_number( plint ,
00135                      "delta" ,   
00136                      -1000 ,    
00137                      1000 ,  
00138                      0 ,    
00139                      0 ,   
00140                      TRUE
00141                    ) ;
00142    PLUTO_add_hint( plint , "Shift data timecourse by delta before splitting and averaging" ) ;
00143 
00144    /*---------- 3rd line: computation ----------*/
00145 
00146    PLUTO_add_option( plint ,
00147                      "Compute" ,  /* label at left of input line */
00148                      "Compute" ,  /* tag to return to plugin */
00149                      TRUE         /* is this mandatory? */
00150                    ) ;
00151 
00152    PLUTO_add_string( plint ,
00153                      "Method" ,           /* label next to chooser button */
00154                      _STAVG_NUM_METHODS,  /* number of strings in list */
00155                      method_strings ,     /* list of strings to choose among */
00156                      _STAVG_METH_MEAN     /* index of default string */
00157                    ) ;
00158 
00159    PLUTO_add_hint( plint , "Choose statistic to compute" ) ;
00160 
00161    /*---------- 4th line --------*/
00162 
00163    PLUTO_add_option( plint ,
00164                      "Parameters" ,  /* label at left of input line */
00165                      "Parameters" ,  /* tag to return to plugin */
00166                      FALSE            /* is this mandatory? */
00167                    ) ;
00168 
00169    PLUTO_add_number( plint ,
00170                      "maxlength" ,    /* label next to chooser */
00171                      0 ,         /* smallest possible value */
00172                      1000 ,        /* largest possible value */
00173                      0 ,         /* decimal shift (none in this case) */
00174                      15 ,         /* default value */
00175                      TRUE       /* allow user to edit value? */
00176                    ) ;
00177    PLUTO_add_hint( plint , "maximum # of timepoints of output dataset" ) ;
00178 
00179    PLUTO_add_string( plint ,
00180                      "no1?" ,               /* label next to chooser button */
00181                      2  ,               /* number of strings to choose among */
00182                      yes_no_strings ,  /* list of strings to choose among */
00183                      1                  /* index of default string */
00184                    ) ;
00185 
00186    PLUTO_add_hint( plint , "ignore timepoints where only one image is in average" ) ;
00187 
00188 
00189    /*--------- done with interface setup ---------*/
00190 
00191    return plint ;
00192 }
00193 
00194 /***************************************************************************
00195   Main routine for this plugin (will be called from AFNI).
00196   If the return string is not NULL, some error transpired, and
00197   AFNI will popup the return string in a message box.
00198 ****************************************************************************/
00199 
00200 /*------------------ macros to return workspace at exit -------------------*/
00201 
00202 #undef  FREEUP
00203 #define FREEUP(x) if((x) != NULL){free((x)); (x)=NULL;}
00204 
00205 #define FREE_WORKSPACE                              \
00206   do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
00207       FREEUP(fout) ; \
00208       FREEUP(fxar) ; \
00209     } while(0)
00210 
00211 /*-------------------------------------------------------------------------*/
00212 
00213 char * STAVG_main( PLUGIN_interface * plint )
00214 {
00215    MCW_idcode * idc ;                          /* input dataset idcode */
00216    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00217    char * new_prefix , * str , * str2;         /* strings from user */
00218    int   meth;                                 /* chosen computation method */
00219    int   new_datum ,                           /* control parameters */
00220          old_datum , ntime ;
00221 
00222    int   te, ne, tinc, kim, nia;
00223    int   numepochs, minlength, maxlength, lastindex, navgpts;
00224    int   nvox , perc , new_units, old_units ;
00225    int   ii, ibot,itop , kk, jj; 
00226    int   no1, user_maxlength, delta;
00227    int   *pEpochLength, *pTimeIndex;
00228    int   nx, ny, nz, npix;
00229    float *pNumAvg;
00230    float old_dtime;
00231 
00232    MRI_IMAGE * stimim;
00233    
00234    MRI_IMARR *avgimar;
00235 
00236    byte   ** bptr  = NULL ;  /* one of these will be the array of */
00237    short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
00238    float  ** fptr  = NULL ;  /* (depending on input datum type) */
00239 
00240    float   * fxar  = NULL ;  /* array loaded from input dataset */
00241    float   * stimar = NULL ;
00242    float  ** fout  = NULL ;  /* will be array of output floats */
00243 
00244    float   * tar   = NULL ;  /* will be array of taper coefficients */
00245 
00246    float   * nstimar;
00247 
00248    /*--------------------------------------------------------------------*/
00249    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00250 
00251    /*--------- go to first input line ---------*/
00252 
00253    PLUTO_next_option(plint) ;
00254 
00255    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00256    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00257    if( old_dset == NULL )
00258       return "*************************\n"
00259              "Cannot find Input Dataset\n"
00260              "*************************"  ;
00261 
00262    ntime = DSET_NUM_TIMES(old_dset) ;
00263    if( ntime < 2 )
00264       return "*****************************\n"
00265              "Dataset has only 1 time point\n"
00266              "*****************************"  ;
00267 
00268    ii = DSET_NVALS_PER_TIME(old_dset) ;
00269    if( ii > 1 )
00270       return "************************************\n"
00271              "Dataset has > 1 value per time point\n"
00272              "************************************"  ;
00273    
00274    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
00275    new_datum = old_datum;
00276    old_dtime = DSET_TIMESTEP(old_dset);
00277    old_units = DSET_TIMEUNITS(old_dset);
00278    
00279    nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
00280    npix = old_dset->daxes->nxx * old_dset->daxes->nyy;
00281    nx = old_dset->daxes->nxx;
00282 
00283 
00284    new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
00285    if( ! PLUTO_prefix_ok(new_prefix) )      /* check if it is OK */
00286       return "************************\n"
00287              "Output Prefix is illegal\n"
00288              "************************"  ;
00289 
00290    /*--------- go to next input line ---------*/
00291 
00292    PLUTO_next_option(plint);
00293 
00294    stimim = PLUTO_get_timeseries(plint);
00295    if( stimim == NULL ) return "Please specify stimulus timing";
00296 
00297    if( stimim->nx < ntime ){
00298       return "**************************************\n"
00299              "Not enough pts in stimulus time-series\n"
00300              "**************************************";
00301    }
00302 
00303    stimar = MRI_FLOAT_PTR(stimim);
00304 
00305 
00306    delta = PLUTO_get_number(plint);
00307 
00308    if( abs(delta) > ntime ){
00309       return "************************\n"
00310              "Delta shift is too large\n"
00311              "************************";
00312    }
00313   
00314    /*initialize variables if not user specified */
00315    user_maxlength = ntime;
00316    no1 = 0;
00317 
00318    /*--------- go to next input line ---------*/
00319 
00320    PLUTO_next_option(plint);
00321 
00322    str  = PLUTO_get_string(plint) ;      /* get string item (the method) */
00323    meth = PLUTO_string_index( str ,      /* find it in list it is from */
00324                               _STAVG_NUM_METHODS ,
00325                               method_strings ) ;
00326 
00327    /*--------- see if the 4th option line is present --------*/
00328 
00329    str = PLUTO_get_optiontag( plint ) ;
00330    if( str != NULL ){
00331       user_maxlength = (int) PLUTO_get_number(plint) ;
00332       str2  = PLUTO_get_string(plint) ;      /* get string item (the method) */
00333       no1   = PLUTO_string_index( str2 ,      /* find it in list it is from */
00334                                  2 ,
00335                                  yes_no_strings) ;
00336    }
00337    
00338 
00339    /*------------------------------------------------------*/
00340    /*---------- At this point, the inputs are OK ----------*/
00341 
00342    PLUTO_popup_meter( plint ) ;  /* popup a progress meter */
00343 
00344    /*________________[ Main Code ]_________________________*/
00345   
00346    fout = avg_epochs( old_dset, stimar, user_maxlength, 1, meth, plint );
00347    if( fout == NULL ) return " \nError in avg_epochs() function!\n " ;
00348    
00349    if( RMB_DEBUG ) fprintf(stderr, "Done with avg_epochs\n");
00350    maxlength = M_maxlength;
00351    
00352    
00353    /*______________________________________________________*/
00354 
00355    
00356    new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
00357 
00358    { char * his = PLUTO_commandstring(plint) ;
00359      tross_Copy_History( old_dset , new_dset ) ;
00360      tross_Append_History( new_dset , his ) ; free( his ) ;
00361    }
00362    
00363    /*-- edit some of its internal parameters --*/
00364    ii = EDIT_dset_items(
00365            new_dset ,
00366               ADN_prefix      , new_prefix ,           /* filename prefix */
00367               ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
00368               ADN_datum_all   , new_datum ,            /* atomic datum */
00369               ADN_nvals       , maxlength ,            /* # sub-bricks */
00370               ADN_ntt         , maxlength ,            /* # time points */
00371            /*   ADN_ttorg       , old_dtime ,  */              /* time origin */
00372            /*   ADN_ttdel       , old_dtime ,  */            /* time step */
00373            /*   ADN_ttdur       , old_dtime ,  */            /* time duration */
00374            /*   ADN_nsl         , 0 ,          */        /* z-axis time slicing */
00375            /*   ADN_tunits      , old_units ,  */        /* time units */
00376            ADN_none ) ;
00377 
00378    if( ii != 0 ){
00379       THD_delete_3dim_dataset( new_dset , False ) ;
00380       FREE_WORKSPACE ;
00381       return "***********************************\n"
00382              "Error while creating output dataset\n"
00383              "***********************************"  ;
00384    }
00385 
00386 
00387    /*------------------------------------------------------------*/
00388    /*------- The output is now in fout[kk][ii],
00389              for kk=0..maxlength-1 , ii=0..nvox-1.
00390              We must now put this into the output dataset -------*/
00391 
00392    switch( new_datum ){
00393 
00394       /*** output is floats is the simplest:
00395            we just have to attach the fout bricks to the dataset ***/
00396 
00397       case MRI_float:
00398          for( kk=0 ; kk < maxlength ; kk++ )
00399             EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
00400       break ;
00401 
00402       /*** output is shorts:
00403            we have to create a scaled sub-brick from fout ***/
00404 
00405       case MRI_short:{
00406          short * bout ;
00407          float fac ; 
00408 
00409          for( kk=0 ; kk < maxlength ; kk++ ){  /* loop over sub-bricks */
00410 
00411             /*-- get output sub-brick --*/
00412             bout = (short *) malloc( sizeof(short) * nvox ) ;
00413             if( bout == NULL ){
00414                fprintf(stderr,"\nFinal malloc error in plug_stavg!\n\a") ;
00415                return("Final malloc error in plug_stavg!"); ;
00416                /*  exit(1) ;*/
00417             }
00418 
00419             /*-- find scaling and then scale --*/
00420             /*fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;*/
00421             fac = 1.0;
00422             EDIT_coerce_scale_type( nvox,fac ,
00423                                     MRI_float,fout[kk] , MRI_short,bout ) ;
00424             free( fout[kk] ) ;  /* don't need this anymore */
00425 
00426             /*-- put output brick into dataset, and store scale factor --*/
00427             EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
00428          }
00429       }
00430       break ;
00431 
00432       /*** output is bytes (byte = unsigned char)
00433            we have to create a scaled sub-brick from fout ***/
00434 
00435       case MRI_byte:{
00436          byte * bout ;
00437          float fac ;
00438 
00439          for( kk=0 ; kk < maxlength ; kk++ ){  /* loop over sub-bricks */
00440 
00441             /*-- get output sub-brick --*/
00442 
00443             bout = (byte *) malloc( sizeof(byte) * nvox ) ;
00444             if( bout == NULL ){
00445                fprintf(stderr,"\nFinal malloc error in plug_stavg!\n\a") ;
00446                return("Final malloc error in plug_stavg!"); ;
00447                /*               exit(1) ;*/
00448             }
00449 
00450             /*-- find scaling and then scale --*/
00451 
00452             fac = 1.0;
00453             EDIT_coerce_scale_type( nvox,fac ,
00454                                     MRI_float,fout[kk] , MRI_byte,bout ) ;
00455 
00456             free( fout[kk] ) ;  /* don't need this anymore */
00457 
00458             /*-- put output brick into dataset, and store scale factor --*/
00459 
00460             EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
00461          }
00462 
00463       }
00464       break ;
00465 
00466    } /* end of switch on output data type */
00467 
00468    /*-------------- Cleanup and go home ----------------*/
00469 
00470    PLUTO_set_meter( plint , 100 ) ;  /* set progress meter to 100% */
00471 
00472    PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00473 
00474    FREE_WORKSPACE ;
00475    return NULL ;  /* null string returned means all was OK */
00476 }
00477 
00478 /*---------------------------------------------------------------*/
00479 float ** avg_epochs( THD_3dim_dataset * dset, float * ref, 
00480                     int user_maxlength, int no1, int meth,
00481                     PLUGIN_interface * plint )
00482 /*---------------------------------------------------------------*/
00483 {
00484 
00485    int     numepochs, lastindex;
00486    int     nvox, numims, nx, ny, nz;
00487    int     kim, ne, te, tinc, nia;
00488    int     ii, kk;
00489    int     maxlength, minlength;
00490    int     datum;
00491    float   fac;       /* scaling factor for current subbrick */
00492    double  scaledval; /* temp var for scaled value to be used repeatedly */
00493    float ** fxar;
00494    float ** outar;    /* output averaged time-series */
00495    float * sumsq = NULL; /* sum of squared voxel values */
00496    float * pNumAvg;  /* array for number of pts to avg at each time*/
00497    int   * pTimeIndex; /* array of time markers (1st img of each epoch) */
00498    int   * pEpochLength; /* array of epoch lengths */
00499    float ** tempar;
00500    MRI_IMARR *inimar;
00501    
00502    
00503    nx = dset->daxes->nxx;
00504    ny = dset->daxes->nyy;
00505    nz = dset->daxes->nzz;
00506    nvox = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
00507    numims = DSET_NUM_TIMES(dset);
00508    datum = DSET_BRICK_TYPE( dset , 0 ) ; /* get old dataset datum type */
00509    
00510    PLUTO_popup_meter(plint) ;
00511    
00512    DSET_load(dset);
00513    
00514    inimar =  dset_to_mri(dset);
00515    if( inimar == NULL ) return NULL ;
00516 
00517    fxar = (float **) malloc( sizeof( float *) * numims);
00518    if( datum == MRI_float){
00519       for( kk=0; kk<numims; kk++){
00520          fxar[kk] = MRI_FLOAT_PTR(IMAGE_IN_IMARR(inimar,kk));
00521       }
00522    }
00523    else{
00524       for( kk=0; kk<numims; kk++){
00525          fxar[kk] = MRI_FLOAT_PTR(mri_to_float(IMAGE_IN_IMARR(inimar,kk)));
00526       }
00527    }
00528 
00529    nia = 0;    /* number of images (timepoints) averaged  where num epochs > 1*/
00530    
00531    if( RMB_DEBUG ) fprintf(stderr, "Start stavg\n");
00532    
00533    /* determine number of epochs to average */
00534    if( RMB_DEBUG ) fprintf(stderr, "Determining number of epochs...");
00535    numepochs = 1;
00536    for( kim=0; kim < numims; kim++ ){
00537       if( ref[kim] > 0) numepochs++;
00538    }
00539    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00540 
00541    /* set array of epoch lengths */
00542    if( RMB_DEBUG ) fprintf(stderr, "Set array of epoch lengths...");
00543    pEpochLength = (int *)malloc(sizeof(int) * numepochs);
00544    for( ne=0; ne < numepochs; ne++) pEpochLength[ne] = 0;
00545    ne = 0;
00546    for( kim=0; kim < numims; kim++ ){
00547       if( ref[kim] > 0) ne++;
00548       pEpochLength[ne]++;
00549    }
00550    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00551 
00552    /* set array of time markers (1st img of each epoch) */
00553    if( RMB_DEBUG ) fprintf(stderr, "Set array of time markers...");
00554    pTimeIndex = (int *)malloc(sizeof(int) * (numepochs - 1));
00555    lastindex = 0;
00556    minlength = numims;
00557    maxlength = 0;
00558    for( ne=0; ne < (numepochs-1); ne++){
00559       pTimeIndex[ne] = lastindex + pEpochLength[ne];
00560       lastindex = pTimeIndex[ne];
00561       if(pEpochLength[ne+1] > maxlength) maxlength = pEpochLength[ne+1];
00562       if(pEpochLength[ne+1] < minlength) minlength = pEpochLength[ne+1];
00563    }
00564    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00565 
00566    if(maxlength > user_maxlength) maxlength = user_maxlength;
00567 
00568    if( RMB_DEBUG ) fprintf(stderr, "init...");
00569    pNumAvg = (float *) malloc( sizeof(float) * maxlength);
00570    outar = (float **) malloc( sizeof(float *) * maxlength);
00571    for( te=0; te < maxlength; te++){
00572       outar[te] = (float *) malloc( sizeof(float) * nvox);
00573       for( ii=0; ii<nvox; ii++){
00574          outar[te][ii] = 0.0;
00575       }
00576    }
00577    if (meth == _STAVG_METH_SIGMA) {
00578        sumsq = (float *) malloc( sizeof(float) * nvox);
00579    }
00580    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00581 
00582    if( RMB_DEBUG ) fprintf(stderr, "Start averaging...");
00583    for( te=0; te < maxlength; te++){
00584       pNumAvg[te] = 0.0;   
00585 
00586      switch(meth) {
00587      default: case _STAVG_METH_MEAN:
00588 
00589       for( ne=0; ne < (numepochs-1); ne++){
00590          tinc = pTimeIndex[ne] + te;
00591          if( te < pEpochLength[ne+1] ){
00592             fac = DSET_BRICK_FACTOR(dset, tinc);
00593             if (fac == 0.0 || fac == 1.0) {
00594                 for( ii=0; ii<nvox; ii++){
00595                     outar[te][ii] += fxar[tinc][ii];
00596                 }
00597             } else {
00598                 for( ii=0; ii<nvox; ii++){
00599                     outar[te][ii] += fxar[tinc][ii] * fac;
00600                 }
00601             }
00602             pNumAvg[te]++;
00603          }
00604       }
00605       for( ii=0; ii<nvox; ii++){
00606          outar[te][ii] = outar[te][ii]/pNumAvg[te];
00607       }
00608       break;
00609 
00610      case _STAVG_METH_SIGMA:
00611 
00612       /* sigma statistic is actually the unbiased standard deviation
00613          calculated with this computational formula:
00614            sqrt( (sum(x^2) - sum(x)^2 / n) / (n-1) )  */
00615       for( ii=0; ii<nvox; ii++){
00616           sumsq[ii] = 0.0;
00617       }
00618       for( ne=0; ne < (numepochs-1); ne++){
00619          tinc = pTimeIndex[ne] + te;
00620          if( te < pEpochLength[ne+1] ){
00621             fac = DSET_BRICK_FACTOR(dset, tinc);
00622             if (fac == 0.0 || fac == 1.0) {
00623                 for( ii=0; ii<nvox; ii++){
00624                     outar[te][ii] += fxar[tinc][ii];
00625                     sumsq[ii] += fxar[tinc][ii] * fxar[tinc][ii];
00626                 }
00627             } else {
00628                 for( ii=0; ii<nvox; ii++){
00629                     scaledval = fxar[tinc][ii] * fac;
00630                     outar[te][ii] += scaledval;
00631                     sumsq[ii] += scaledval * scaledval;
00632                 }
00633             }
00634             pNumAvg[te]++;
00635          }
00636       }
00637       for( ii=0; ii<nvox; ii++){
00638           outar[te][ii] = sqrt( (sumsq[ii] -
00639                                  outar[te][ii] * outar[te][ii] /
00640                                  pNumAvg[te]) /
00641                                 (pNumAvg[te] - 1) );
00642       }
00643       break;
00644 
00645      }
00646 
00647       if( pNumAvg[te] > 1) nia ++;
00648       PLUTO_set_meter(plint, (100*(te+1))/maxlength ) ;
00649    }
00650    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00651    if ( no1 ){     /* ignore images with only one average */
00652       if( nia < maxlength) maxlength = nia;
00653    }
00654    
00655    M_maxlength = maxlength;
00656    
00657    
00658    if( RMB_DEBUG ) fprintf(stderr, "malloc output...");
00659    tempar = (float **) malloc(sizeof(float *) * maxlength);
00660    for( te=0 ; te < maxlength ; te++ ){
00661       tempar[te] = (float *) malloc( sizeof(float) * nvox);
00662    }
00663    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00664    
00665    if( RMB_DEBUG ) fprintf(stderr, "convert to output...");
00666    for( te=0; te < maxlength; te++){
00667       for( ii=0; ii<nvox; ii++){
00668          tempar[te][ii] = outar[te][ii];
00669       }
00670    }
00671    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00672 
00673    /* toss arrays */
00674    if( RMB_DEBUG ) fprintf(stderr, "free mem...");
00675    FREE_IMARR( inimar );
00676    if (meth == _STAVG_METH_SIGMA) free(sumsq);
00677    free( outar );
00678    free( pNumAvg );
00679    free( pTimeIndex );
00680    free( pEpochLength );
00681    if( RMB_DEBUG ) fprintf(stderr, "done\n");
00682    DSET_unload(dset);
00683    
00684    return(tempar);
00685 }
00686 
00687 
00688 /*-------------------------------------------------------*/
00689 MRI_IMARR * dset_to_mri(THD_3dim_dataset * dset)
00690 /*--------------------------------------------------------*/
00691 {
00692 
00693    int ii, kk, ntime, datum;
00694    int nvox, nx, ny, nz;
00695    int use_fac;
00696    
00697    MRI_IMARR * ims_in;
00698    MRI_IMAGE * im, *temp_im;
00699    
00700 
00701    byte   ** bptr  = NULL ;  /* one of these will be the array of */
00702    short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
00703    float  ** fptr  = NULL ;  /* (depending on input datum type) */
00704    
00705    float * fac  = NULL ;  /* array of brick scaling factors */
00706    
00707    float * fout;
00708    
00709 
00710    ntime = DSET_NUM_TIMES(dset) ;
00711    nx = dset->daxes->nxx;
00712    ny = dset->daxes->nyy;
00713    nz = dset->daxes->nzz;
00714    nvox = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
00715    datum = DSET_BRICK_TYPE( dset , 0 ) ; /* get dataset datum type */
00716 
00717    switch( datum ){  /* pointer type depends on input datum type */
00718 
00719       default:
00720          return NULL  ;
00721 
00722       /** create array of pointers into old dataset sub-bricks **/
00723 
00724       /*--------- input is bytes ----------*/
00725       /* voxel #i at time #k is bptr[k][i] */
00726       /* for i=0..nvox-1 and k=0..ntime-1.  */
00727 
00728       case MRI_byte:
00729          bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
00730          if( bptr == NULL ) return NULL ;
00731          for( kk=0 ; kk < ntime ; kk++ )
00732             bptr[kk] = (byte *) DSET_ARRAY(dset,kk) ;
00733       break ;
00734 
00735       /*--------- input is shorts ---------*/
00736       /* voxel #i at time #k is sptr[k][i] */
00737       /* for i=0..nvox-1 and k=0..ntime-1.  */
00738 
00739       case MRI_short:
00740          sptr = (short **) malloc( sizeof(short *) * ntime ) ;
00741          if( sptr == NULL ) return NULL ;
00742          for( kk=0 ; kk < ntime; kk++ )
00743             sptr[kk] = (short *) DSET_ARRAY(dset,kk) ;
00744       break ;
00745 
00746       /*--------- input is floats ---------*/
00747       /* voxel #i at time #k is fptr[k][i] */
00748       /* for i=0..nvox-1 and k=0..ntime-1.  */
00749 
00750       case MRI_float:
00751          fptr = (float **) malloc( sizeof(float *) * ntime) ;
00752          if( fptr == NULL ) return NULL ;
00753          for( kk=0 ; kk < ntime; kk++ )
00754             fptr[kk] = (float *) DSET_ARRAY(dset,kk) ;
00755       break ;
00756 
00757    } /* end of switch on input type */
00758    
00759    INIT_IMARR(ims_in) ;
00760    for( kk=0 ; kk < ntime ; kk++ ){
00761       im = mri_new_vol_empty( nx , ny , nz , datum ) ;
00762       ADDTO_IMARR(ims_in,im) ;
00763    }
00764    
00765    for( kk=0 ; kk < ntime ; kk++ ){
00766       im = IMARR_SUBIMAGE(ims_in,kk) ;
00767       
00768       switch( datum ){
00769          case MRI_byte:  mri_fix_data_pointer( bptr[kk], im ) ; break ;
00770          case MRI_short: mri_fix_data_pointer( sptr[kk], im ) ; break ;
00771          case MRI_float: mri_fix_data_pointer( fptr[kk], im ) ; break ;
00772       }
00773    }
00774 
00775 
00776    
00777    return(ims_in);
00778 }
00779 
00780