00001
00002
00003
00004
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
00019
00020
00021
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
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
00048
00049 char * STAVG_main( PLUGIN_interface * ) ;
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
00059
00060 static PLUGIN_interface * global_plint = NULL ;
00061 int M_maxlength;
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077 DEFINE_PLUGIN_PROTOTYPE
00078
00079 PLUGIN_interface * PLUGIN_init( int ncall )
00080 {
00081 PLUGIN_interface * plint ;
00082
00083 if( ncall > 0 ) return NULL ;
00084
00085
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 ;
00095
00096 PLUTO_set_sequence( plint , "z:Birn" ) ;
00097
00098
00099
00100 PLUTO_add_option( plint ,
00101 "Datasets" ,
00102 "Datasets" ,
00103 TRUE
00104 ) ;
00105
00106 PLUTO_add_dataset( plint ,
00107 "Input" ,
00108 ANAT_ALL_MASK ,
00109 FUNC_FIM_MASK ,
00110 DIMEN_4D_MASK |
00111 BRICK_ALLREAL_MASK
00112 ) ;
00113 PLUTO_add_hint( plint , "Input 3d+t dataset" ) ;
00114
00115 PLUTO_add_string( plint ,
00116 "Output" ,
00117 0,NULL ,
00118 19
00119 ) ;
00120 PLUTO_add_hint( plint , "Name of output dataset" ) ;
00121
00122
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
00145
00146 PLUTO_add_option( plint ,
00147 "Compute" ,
00148 "Compute" ,
00149 TRUE
00150 ) ;
00151
00152 PLUTO_add_string( plint ,
00153 "Method" ,
00154 _STAVG_NUM_METHODS,
00155 method_strings ,
00156 _STAVG_METH_MEAN
00157 ) ;
00158
00159 PLUTO_add_hint( plint , "Choose statistic to compute" ) ;
00160
00161
00162
00163 PLUTO_add_option( plint ,
00164 "Parameters" ,
00165 "Parameters" ,
00166 FALSE
00167 ) ;
00168
00169 PLUTO_add_number( plint ,
00170 "maxlength" ,
00171 0 ,
00172 1000 ,
00173 0 ,
00174 15 ,
00175 TRUE
00176 ) ;
00177 PLUTO_add_hint( plint , "maximum # of timepoints of output dataset" ) ;
00178
00179 PLUTO_add_string( plint ,
00180 "no1?" ,
00181 2 ,
00182 yes_no_strings ,
00183 1
00184 ) ;
00185
00186 PLUTO_add_hint( plint , "ignore timepoints where only one image is in average" ) ;
00187
00188
00189
00190
00191 return plint ;
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
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 ;
00216 THD_3dim_dataset * old_dset , * new_dset ;
00217 char * new_prefix , * str , * str2;
00218 int meth;
00219 int new_datum ,
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 ;
00237 short ** sptr = NULL ;
00238 float ** fptr = NULL ;
00239
00240 float * fxar = NULL ;
00241 float * stimar = NULL ;
00242 float ** fout = NULL ;
00243
00244 float * tar = NULL ;
00245
00246 float * nstimar;
00247
00248
00249
00250
00251
00252
00253 PLUTO_next_option(plint) ;
00254
00255 idc = PLUTO_get_idcode(plint) ;
00256 old_dset = PLUTO_find_dset(idc) ;
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 ) ;
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) ;
00285 if( ! PLUTO_prefix_ok(new_prefix) )
00286 return "************************\n"
00287 "Output Prefix is illegal\n"
00288 "************************" ;
00289
00290
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
00315 user_maxlength = ntime;
00316 no1 = 0;
00317
00318
00319
00320 PLUTO_next_option(plint);
00321
00322 str = PLUTO_get_string(plint) ;
00323 meth = PLUTO_string_index( str ,
00324 _STAVG_NUM_METHODS ,
00325 method_strings ) ;
00326
00327
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) ;
00333 no1 = PLUTO_string_index( str2 ,
00334 2 ,
00335 yes_no_strings) ;
00336 }
00337
00338
00339
00340
00341
00342 PLUTO_popup_meter( plint ) ;
00343
00344
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 ) ;
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
00364 ii = EDIT_dset_items(
00365 new_dset ,
00366 ADN_prefix , new_prefix ,
00367 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
00368 ADN_datum_all , new_datum ,
00369 ADN_nvals , maxlength ,
00370 ADN_ntt , maxlength ,
00371
00372
00373
00374
00375
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
00389
00390
00391
00392 switch( new_datum ){
00393
00394
00395
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
00403
00404
00405 case MRI_short:{
00406 short * bout ;
00407 float fac ;
00408
00409 for( kk=0 ; kk < maxlength ; kk++ ){
00410
00411
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
00417 }
00418
00419
00420
00421 fac = 1.0;
00422 EDIT_coerce_scale_type( nvox,fac ,
00423 MRI_float,fout[kk] , MRI_short,bout ) ;
00424 free( fout[kk] ) ;
00425
00426
00427 EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
00428 }
00429 }
00430 break ;
00431
00432
00433
00434
00435 case MRI_byte:{
00436 byte * bout ;
00437 float fac ;
00438
00439 for( kk=0 ; kk < maxlength ; kk++ ){
00440
00441
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
00448 }
00449
00450
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] ) ;
00457
00458
00459
00460 EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
00461 }
00462
00463 }
00464 break ;
00465
00466 }
00467
00468
00469
00470 PLUTO_set_meter( plint , 100 ) ;
00471
00472 PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00473
00474 FREE_WORKSPACE ;
00475 return NULL ;
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;
00492 double scaledval;
00493 float ** fxar;
00494 float ** outar;
00495 float * sumsq = NULL;
00496 float * pNumAvg;
00497 int * pTimeIndex;
00498 int * pEpochLength;
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 ) ;
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;
00530
00531 if( RMB_DEBUG ) fprintf(stderr, "Start stavg\n");
00532
00533
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
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
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
00613
00614
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 ){
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
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 ;
00702 short ** sptr = NULL ;
00703 float ** fptr = NULL ;
00704
00705 float * fac = NULL ;
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 ) ;
00716
00717 switch( datum ){
00718
00719 default:
00720 return NULL ;
00721
00722
00723
00724
00725
00726
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
00736
00737
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
00747
00748
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 }
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