1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "afni.h"
8
9 #ifndef ALLOW_PLUGINS
10 # error "Plugins not properly set up -- see machdep.h"
11 #endif
12
13 #define TEST_VOXEL 6177
14 #define TEST_TIME 0
15 #define RMB_DEBUG 0
16
17 /***********************************************************************
18 Plugin that averages epochs from single trial data
19 ************************************************************************/
20
21 /*--------------------- string to 'help' the user --------------------*/
22
23 static char helpstring[] =
24 "Purpose: Averaging epochs of single trial data\n"
25 "\n"
26 "Input items to this plugin are:\n"
27 " Datasets: Input = 3D+time dataset to process\n"
28 " = reference time-series\n"
29 " Output = Prefix for new dataset\n"
30 " Additional Parameters\n"
31 " delta = shift timeseries by delta before splitting and averaging\n"
32 " method = type of statistic to calculate\n"
33 " maxlength = maximum avg ts length\n"
34 " no1? = images w/ only one img in avg ignored\n"
35 "Author -- RM Birn"
36 ;
37
38 /*------------- strings for output format -------------*/
39
40 static char * yes_no_strings[] = { "No" , "Yes" } ;
41 static char * method_strings[] = { "Mean" , "Sigma" } ;
42
43 #define _STAVG_NUM_METHODS (sizeof(method_strings)/sizeof(char *))
44 #define _STAVG_METH_MEAN 0
45 #define _STAVG_METH_SIGMA 1
46
47 /*--------------- prototypes for internal routines ---------------*/
48
49 char * STAVG_main( PLUGIN_interface * ) ; /* the entry point */
50
51 float ** avg_epochs( THD_3dim_dataset * dset, float * ref,
52 int user_maxlength, int no1, int meth,
53 PLUGIN_interface *plint );
54
55 MRI_IMARR * dset_to_mri(THD_3dim_dataset * dset);
56
57
58 /*---------------------------- global data ---------------------------*/
59
60 static PLUGIN_interface * global_plint = NULL ;
61 int M_maxlength;
62
63 /***********************************************************************
64 Set up the interface to the user:
65 1) Create a new interface using "PLUTO_new_interface";
66
67 2) For each line of inputs, create the line with "PLUTO_add_option"
68 (this line of inputs can be optional or mandatory);
69
70 3) For each item on the line, create the item with
71 "PLUTO_add_dataset" for a dataset chooser,
72 "PLUTO_add_string" for a string chooser,
73 "PLUTO_add_number" for a number chooser.
74 ************************************************************************/
75
76
77 DEFINE_PLUGIN_PROTOTYPE
78
PLUGIN_init(int ncall)79 PLUGIN_interface * PLUGIN_init( int ncall )
80 {
81 PLUGIN_interface * plint ; /* will be the output of this routine */
82
83 if( ncall > 0 ) return NULL ; /* one interfaces */
84
85 CHECK_IF_ALLOWED("SINGLETRIALAVG","SingleTrial Avg") ; /* 30 Sep 2016 */
86
87 /*---------------- set titles and call point ----------------*/
88
89 plint = PLUTO_new_interface( "SingleTrial Avg" ,
90 "Averaging of epochs in Single Trial data" ,
91 helpstring ,
92 PLUGIN_CALL_VIA_MENU , STAVG_main ) ;
93
94 PLUTO_add_hint( plint , "Averaging of epochs in Single Trial data" ) ;
95
96 global_plint = plint ; /* make global copy */
97
98 PLUTO_set_sequence( plint , "z:Birn" ) ;
99
100 /*--------- 1st line ---------*/
101
102 PLUTO_add_option( plint ,
103 "Datasets" , /* label at left of input line */
104 "Datasets" , /* tag to return to plugin */
105 TRUE /* is this mandatory? */
106 ) ;
107
108 PLUTO_add_dataset( plint ,
109 "Input" , /* label next to button */
110 ANAT_ALL_MASK , /* take any anat datasets */
111 FUNC_FIM_MASK , /* only allow fim funcs */
112 DIMEN_4D_MASK | /* need 3D+time datasets */
113 BRICK_ALLREAL_MASK /* need real-valued datasets */
114 ) ;
115 PLUTO_add_hint( plint , "Input 3d+t dataset" ) ;
116
117 PLUTO_add_string( plint ,
118 "Output" , /* label next to textfield */
119 0,NULL , /* no fixed strings to choose among */
120 19 /* 19 spaces for typing in value */
121 ) ;
122 PLUTO_add_hint( plint , "Name of output dataset" ) ;
123
124 /*---------- 2nd line --------*/
125
126 PLUTO_add_option( plint ,
127 "Timing" ,
128 "Timing" ,
129 TRUE
130 ) ;
131
132
133 PLUTO_add_timeseries(plint, "Stim. Timing");
134 PLUTO_add_hint( plint , "Stimulus Timing (0 = no task, 1 = task)" ) ;
135
136 PLUTO_add_number( plint ,
137 "delta" ,
138 -1000 ,
139 1000 ,
140 0 ,
141 0 ,
142 TRUE
143 ) ;
144 PLUTO_add_hint( plint , "Shift data timecourse by delta before splitting and averaging" ) ;
145
146 /*---------- 3rd line: computation ----------*/
147
148 PLUTO_add_option( plint ,
149 "Compute" , /* label at left of input line */
150 "Compute" , /* tag to return to plugin */
151 TRUE /* is this mandatory? */
152 ) ;
153
154 PLUTO_add_string( plint ,
155 "Method" , /* label next to chooser button */
156 _STAVG_NUM_METHODS, /* number of strings in list */
157 method_strings , /* list of strings to choose among */
158 _STAVG_METH_MEAN /* index of default string */
159 ) ;
160
161 PLUTO_add_hint( plint , "Choose statistic to compute" ) ;
162
163 /*---------- 4th line --------*/
164
165 PLUTO_add_option( plint ,
166 "Parameters" , /* label at left of input line */
167 "Parameters" , /* tag to return to plugin */
168 FALSE /* is this mandatory? */
169 ) ;
170
171 PLUTO_add_number( plint ,
172 "maxlength" , /* label next to chooser */
173 0 , /* smallest possible value */
174 1000 , /* largest possible value */
175 0 , /* decimal shift (none in this case) */
176 15 , /* default value */
177 TRUE /* allow user to edit value? */
178 ) ;
179 PLUTO_add_hint( plint , "maximum # of timepoints of output dataset" ) ;
180
181 PLUTO_add_string( plint ,
182 "no1?" , /* label next to chooser button */
183 2 , /* number of strings to choose among */
184 yes_no_strings , /* list of strings to choose among */
185 1 /* index of default string */
186 ) ;
187
188 PLUTO_add_hint( plint , "ignore timepoints where only one image is in average" ) ;
189
190
191 /*--------- done with interface setup ---------*/
192
193 return plint ;
194 }
195
196 /***************************************************************************
197 Main routine for this plugin (will be called from AFNI).
198 If the return string is not NULL, some error transpired, and
199 AFNI will popup the return string in a message box.
200 ****************************************************************************/
201
202 /*------------------ macros to return workspace at exit -------------------*/
203
204 #undef FREEUP
205 #define FREEUP(x) if((x) != NULL){free((x)); (x)=NULL;}
206
207 #define FREE_WORKSPACE \
208 do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
209 FREEUP(fout) ; \
210 FREEUP(fxar) ; \
211 } while(0)
212
213 /*-------------------------------------------------------------------------*/
214
STAVG_main(PLUGIN_interface * plint)215 char * STAVG_main( PLUGIN_interface * plint )
216 {
217 MCW_idcode * idc ; /* input dataset idcode */
218 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
219 char * new_prefix , * str , * str2; /* strings from user */
220 int meth; /* chosen computation method */
221 int new_datum , /* control parameters */
222 old_datum , ntime ;
223
224 int te, ne, tinc, kim, nia;
225 int numepochs, minlength, maxlength, lastindex, navgpts;
226 int nvox , perc , new_units, old_units ;
227 int ii, ibot,itop , kk, jj;
228 int no1, user_maxlength, delta;
229 int *pEpochLength, *pTimeIndex;
230 int nx, ny, nz, npix;
231 float *pNumAvg;
232 float old_dtime;
233
234 MRI_IMAGE * stimim;
235
236 MRI_IMARR *avgimar;
237
238 byte ** bptr = NULL ; /* one of these will be the array of */
239 short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
240 float ** fptr = NULL ; /* (depending on input datum type) */
241
242 float * fxar = NULL ; /* array loaded from input dataset */
243 float * stimar = NULL ;
244 float ** fout = NULL ; /* will be array of output floats */
245
246 float * tar = NULL ; /* will be array of taper coefficients */
247
248 float * nstimar;
249
250 /*--------------------------------------------------------------------*/
251 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
252
253 /*--------- go to first input line ---------*/
254
255 PLUTO_next_option(plint) ;
256
257 idc = PLUTO_get_idcode(plint) ; /* get dataset item */
258 old_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
259 if( old_dset == NULL )
260 return "*************************\n"
261 "Cannot find Input Dataset\n"
262 "*************************" ;
263
264 ntime = DSET_NUM_TIMES(old_dset) ;
265 if( ntime < 2 )
266 return "*****************************\n"
267 "Dataset has only 1 time point\n"
268 "*****************************" ;
269
270 ii = DSET_NVALS_PER_TIME(old_dset) ;
271 if( ii > 1 )
272 return "************************************\n"
273 "Dataset has > 1 value per time point\n"
274 "************************************" ;
275
276 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
277 new_datum = old_datum;
278 old_dtime = DSET_TIMESTEP(old_dset);
279 old_units = DSET_TIMEUNITS(old_dset);
280
281 nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
282 npix = old_dset->daxes->nxx * old_dset->daxes->nyy;
283 nx = old_dset->daxes->nxx;
284
285
286 new_prefix = PLUTO_get_string(plint) ; /* get string item (the output prefix) */
287 if( ! PLUTO_prefix_ok(new_prefix) ) /* check if it is OK */
288 return "************************\n"
289 "Output Prefix is illegal\n"
290 "************************" ;
291
292 /*--------- go to next input line ---------*/
293
294 PLUTO_next_option(plint);
295
296 stimim = PLUTO_get_timeseries(plint);
297 if( stimim == NULL ) return "Please specify stimulus timing";
298
299 if( stimim->nx < ntime ){
300 return "**************************************\n"
301 "Not enough pts in stimulus time-series\n"
302 "**************************************";
303 }
304
305 stimar = MRI_FLOAT_PTR(stimim);
306
307
308 delta = PLUTO_get_number(plint);
309
310 if( abs(delta) > ntime ){
311 return "************************\n"
312 "Delta shift is too large\n"
313 "************************";
314 }
315
316 /*initialize variables if not user specified */
317 user_maxlength = ntime;
318 no1 = 0;
319
320 /*--------- go to next input line ---------*/
321
322 PLUTO_next_option(plint);
323
324 str = PLUTO_get_string(plint) ; /* get string item (the method) */
325 meth = PLUTO_string_index( str , /* find it in list it is from */
326 _STAVG_NUM_METHODS ,
327 method_strings ) ;
328
329 /*--------- see if the 4th option line is present --------*/
330
331 str = PLUTO_get_optiontag( plint ) ;
332 if( str != NULL ){
333 user_maxlength = (int) PLUTO_get_number(plint) ;
334 str2 = PLUTO_get_string(plint) ; /* get string item (the method) */
335 no1 = PLUTO_string_index( str2 , /* find it in list it is from */
336 2 ,
337 yes_no_strings) ;
338 }
339
340
341 /*------------------------------------------------------*/
342 /*---------- At this point, the inputs are OK ----------*/
343
344 PLUTO_popup_meter( plint ) ; /* popup a progress meter */
345
346 /*________________[ Main Code ]_________________________*/
347
348 fout = avg_epochs( old_dset, stimar, user_maxlength, 1, meth, plint );
349 if( fout == NULL ) return " \nError in avg_epochs() function!\n " ;
350
351 if( RMB_DEBUG ) fprintf(stderr, "Done with avg_epochs\n");
352 maxlength = M_maxlength;
353
354
355 /*______________________________________________________*/
356
357
358 new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
359
360 { char * his = PLUTO_commandstring(plint) ;
361 tross_Copy_History( old_dset , new_dset ) ;
362 tross_Append_History( new_dset , his ) ; free( his ) ;
363 }
364
365 /*-- edit some of its internal parameters --*/
366 ii = EDIT_dset_items(
367 new_dset ,
368 ADN_prefix , new_prefix , /* filename prefix */
369 ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
370 ADN_datum_all , new_datum , /* atomic datum */
371 ADN_nvals , maxlength , /* # sub-bricks */
372 ADN_ntt , maxlength , /* # time points */
373 /* ADN_ttorg , old_dtime , */ /* time origin */
374 /* ADN_ttdel , old_dtime , */ /* time step */
375 /* ADN_ttdur , old_dtime , */ /* time duration */
376 /* ADN_nsl , 0 , */ /* z-axis time slicing */
377 /* ADN_tunits , old_units , */ /* time units */
378 ADN_none ) ;
379
380 if( ii != 0 ){
381 THD_delete_3dim_dataset( new_dset , False ) ;
382 FREE_WORKSPACE ;
383 return "***********************************\n"
384 "Error while creating output dataset\n"
385 "***********************************" ;
386 }
387
388
389 /*------------------------------------------------------------*/
390 /*------- The output is now in fout[kk][ii],
391 for kk=0..maxlength-1 , ii=0..nvox-1.
392 We must now put this into the output dataset -------*/
393
394 switch( new_datum ){
395
396 /*** output is floats is the simplest:
397 we just have to attach the fout bricks to the dataset ***/
398
399 case MRI_float:
400 for( kk=0 ; kk < maxlength ; kk++ )
401 EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
402 break ;
403
404 /*** output is shorts:
405 we have to create a scaled sub-brick from fout ***/
406
407 case MRI_short:{
408 short * bout ;
409 float fac ;
410
411 for( kk=0 ; kk < maxlength ; kk++ ){ /* loop over sub-bricks */
412
413 /*-- get output sub-brick --*/
414 bout = (short *) malloc( sizeof(short) * nvox ) ;
415 if( bout == NULL ){
416 fprintf(stderr,"\nFinal malloc error in plug_stavg!\n\a") ;
417 return("Final malloc error in plug_stavg!"); ;
418 /* exit(1) ;*/
419 }
420
421 /*-- find scaling and then scale --*/
422 /*fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;*/
423 fac = 1.0;
424 EDIT_coerce_scale_type( nvox,fac ,
425 MRI_float,fout[kk] , MRI_short,bout ) ;
426 free( fout[kk] ) ; /* don't need this anymore */
427
428 /*-- put output brick into dataset, and store scale factor --*/
429 EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
430 }
431 }
432 break ;
433
434 /*** output is bytes (byte = unsigned char)
435 we have to create a scaled sub-brick from fout ***/
436
437 case MRI_byte:{
438 byte * bout ;
439 float fac ;
440
441 for( kk=0 ; kk < maxlength ; kk++ ){ /* loop over sub-bricks */
442
443 /*-- get output sub-brick --*/
444
445 bout = (byte *) malloc( sizeof(byte) * nvox ) ;
446 if( bout == NULL ){
447 fprintf(stderr,"\nFinal malloc error in plug_stavg!\n\a") ;
448 return("Final malloc error in plug_stavg!"); ;
449 /* exit(1) ;*/
450 }
451
452 /*-- find scaling and then scale --*/
453
454 fac = 1.0;
455 EDIT_coerce_scale_type( nvox,fac ,
456 MRI_float,fout[kk] , MRI_byte,bout ) ;
457
458 free( fout[kk] ) ; /* don't need this anymore */
459
460 /*-- put output brick into dataset, and store scale factor --*/
461
462 EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
463 }
464
465 }
466 break ;
467
468 } /* end of switch on output data type */
469
470 /*-------------- Cleanup and go home ----------------*/
471
472 PLUTO_set_meter( plint , 100 ) ; /* set progress meter to 100% */
473
474 PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
475
476 FREE_WORKSPACE ;
477 return NULL ; /* null string returned means all was OK */
478 }
479
480 /*---------------------------------------------------------------*/
avg_epochs(THD_3dim_dataset * dset,float * ref,int user_maxlength,int no1,int meth,PLUGIN_interface * plint)481 float ** avg_epochs( THD_3dim_dataset * dset, float * ref,
482 int user_maxlength, int no1, int meth,
483 PLUGIN_interface * plint )
484 /*---------------------------------------------------------------*/
485 {
486
487 int numepochs, lastindex;
488 int nvox, numims, nx, ny, nz;
489 int kim, ne, te, tinc, nia;
490 int ii, kk;
491 int maxlength, minlength;
492 int datum;
493 float fac; /* scaling factor for current subbrick */
494 double scaledval; /* temp var for scaled value to be used repeatedly */
495 float ** fxar;
496 float ** outar; /* output averaged time-series */
497 float * sumsq = NULL; /* sum of squared voxel values */
498 float * pNumAvg; /* array for number of pts to avg at each time*/
499 int * pTimeIndex; /* array of time markers (1st img of each epoch) */
500 int * pEpochLength; /* array of epoch lengths */
501 float ** tempar;
502 MRI_IMARR *inimar;
503
504
505 nx = dset->daxes->nxx;
506 ny = dset->daxes->nyy;
507 nz = dset->daxes->nzz;
508 nvox = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
509 numims = DSET_NUM_TIMES(dset);
510 datum = DSET_BRICK_TYPE( dset , 0 ) ; /* get old dataset datum type */
511
512 PLUTO_popup_meter(plint) ;
513
514 DSET_load(dset);
515
516 inimar = dset_to_mri(dset);
517 if( inimar == NULL ) return NULL ;
518
519 fxar = (float **) malloc( sizeof( float *) * numims);
520 if( datum == MRI_float){
521 for( kk=0; kk<numims; kk++){
522 fxar[kk] = MRI_FLOAT_PTR(IMAGE_IN_IMARR(inimar,kk));
523 }
524 }
525 else{
526 for( kk=0; kk<numims; kk++){
527 fxar[kk] = MRI_FLOAT_PTR(mri_to_float(IMAGE_IN_IMARR(inimar,kk)));
528 }
529 }
530
531 nia = 0; /* number of images (timepoints) averaged where num epochs > 1*/
532
533 if( RMB_DEBUG ) fprintf(stderr, "Start stavg\n");
534
535 /* determine number of epochs to average */
536 if( RMB_DEBUG ) fprintf(stderr, "Determining number of epochs...");
537 numepochs = 1;
538 for( kim=0; kim < numims; kim++ ){
539 if( ref[kim] > 0) numepochs++;
540 }
541 if( RMB_DEBUG ) fprintf(stderr, "done\n");
542
543 /* set array of epoch lengths */
544 if( RMB_DEBUG ) fprintf(stderr, "Set array of epoch lengths...");
545 pEpochLength = (int *)malloc(sizeof(int) * numepochs);
546 for( ne=0; ne < numepochs; ne++) pEpochLength[ne] = 0;
547 ne = 0;
548 for( kim=0; kim < numims; kim++ ){
549 if( ref[kim] > 0) ne++;
550 pEpochLength[ne]++;
551 }
552 if( RMB_DEBUG ) fprintf(stderr, "done\n");
553
554 /* set array of time markers (1st img of each epoch) */
555 if( RMB_DEBUG ) fprintf(stderr, "Set array of time markers...");
556 pTimeIndex = (int *)malloc(sizeof(int) * (numepochs - 1));
557 lastindex = 0;
558 minlength = numims;
559 maxlength = 0;
560 for( ne=0; ne < (numepochs-1); ne++){
561 pTimeIndex[ne] = lastindex + pEpochLength[ne];
562 lastindex = pTimeIndex[ne];
563 if(pEpochLength[ne+1] > maxlength) maxlength = pEpochLength[ne+1];
564 if(pEpochLength[ne+1] < minlength) minlength = pEpochLength[ne+1];
565 }
566 if( RMB_DEBUG ) fprintf(stderr, "done\n");
567
568 if(maxlength > user_maxlength) maxlength = user_maxlength;
569
570 if( RMB_DEBUG ) fprintf(stderr, "init...");
571 pNumAvg = (float *) malloc( sizeof(float) * maxlength);
572 outar = (float **) malloc( sizeof(float *) * maxlength);
573 for( te=0; te < maxlength; te++){
574 outar[te] = (float *) malloc( sizeof(float) * nvox);
575 for( ii=0; ii<nvox; ii++){
576 outar[te][ii] = 0.0;
577 }
578 }
579 if (meth == _STAVG_METH_SIGMA) {
580 sumsq = (float *) malloc( sizeof(float) * nvox);
581 }
582 if( RMB_DEBUG ) fprintf(stderr, "done\n");
583
584 if( RMB_DEBUG ) fprintf(stderr, "Start averaging...");
585 for( te=0; te < maxlength; te++){
586 pNumAvg[te] = 0.0;
587
588 switch(meth) {
589 default: case _STAVG_METH_MEAN:
590
591 for( ne=0; ne < (numepochs-1); ne++){
592 tinc = pTimeIndex[ne] + te;
593 if( te < pEpochLength[ne+1] ){
594 fac = DSET_BRICK_FACTOR(dset, tinc);
595 if (fac == 0.0 || fac == 1.0) {
596 for( ii=0; ii<nvox; ii++){
597 outar[te][ii] += fxar[tinc][ii];
598 }
599 } else {
600 for( ii=0; ii<nvox; ii++){
601 outar[te][ii] += fxar[tinc][ii] * fac;
602 }
603 }
604 pNumAvg[te]++;
605 }
606 }
607 for( ii=0; ii<nvox; ii++){
608 outar[te][ii] = outar[te][ii]/pNumAvg[te];
609 }
610 break;
611
612 case _STAVG_METH_SIGMA:
613
614 /* sigma statistic is actually the unbiased standard deviation
615 calculated with this computational formula:
616 sqrt( (sum(x^2) - sum(x)^2 / n) / (n-1) ) */
617 for( ii=0; ii<nvox; ii++){
618 sumsq[ii] = 0.0;
619 }
620 for( ne=0; ne < (numepochs-1); ne++){
621 tinc = pTimeIndex[ne] + te;
622 if( te < pEpochLength[ne+1] ){
623 fac = DSET_BRICK_FACTOR(dset, tinc);
624 if (fac == 0.0 || fac == 1.0) {
625 for( ii=0; ii<nvox; ii++){
626 outar[te][ii] += fxar[tinc][ii];
627 sumsq[ii] += fxar[tinc][ii] * fxar[tinc][ii];
628 }
629 } else {
630 for( ii=0; ii<nvox; ii++){
631 scaledval = fxar[tinc][ii] * fac;
632 outar[te][ii] += scaledval;
633 sumsq[ii] += scaledval * scaledval;
634 }
635 }
636 pNumAvg[te]++;
637 }
638 }
639 for( ii=0; ii<nvox; ii++){
640 outar[te][ii] = sqrt( (sumsq[ii] -
641 outar[te][ii] * outar[te][ii] /
642 pNumAvg[te]) /
643 (pNumAvg[te] - 1) );
644 }
645 break;
646
647 }
648
649 if( pNumAvg[te] > 1) nia ++;
650 PLUTO_set_meter(plint, (100*(te+1))/maxlength ) ;
651 }
652 if( RMB_DEBUG ) fprintf(stderr, "done\n");
653 if ( no1 ){ /* ignore images with only one average */
654 if( nia < maxlength) maxlength = nia;
655 }
656
657 M_maxlength = maxlength;
658
659
660 if( RMB_DEBUG ) fprintf(stderr, "malloc output...");
661 tempar = (float **) malloc(sizeof(float *) * maxlength);
662 for( te=0 ; te < maxlength ; te++ ){
663 tempar[te] = (float *) malloc( sizeof(float) * nvox);
664 }
665 if( RMB_DEBUG ) fprintf(stderr, "done\n");
666
667 if( RMB_DEBUG ) fprintf(stderr, "convert to output...");
668 for( te=0; te < maxlength; te++){
669 for( ii=0; ii<nvox; ii++){
670 tempar[te][ii] = outar[te][ii];
671 }
672 }
673 if( RMB_DEBUG ) fprintf(stderr, "done\n");
674
675 /* toss arrays */
676 if( RMB_DEBUG ) fprintf(stderr, "free mem...");
677 FREE_IMARR( inimar );
678 if (meth == _STAVG_METH_SIGMA) free(sumsq);
679 free( outar );
680 free( pNumAvg );
681 free( pTimeIndex );
682 free( pEpochLength );
683 if( RMB_DEBUG ) fprintf(stderr, "done\n");
684 DSET_unload(dset);
685
686 return(tempar);
687 }
688
689
690 /*-------------------------------------------------------*/
dset_to_mri(THD_3dim_dataset * dset)691 MRI_IMARR * dset_to_mri(THD_3dim_dataset * dset)
692 /*--------------------------------------------------------*/
693 {
694
695 int ii, kk, ntime, datum;
696 int nvox, nx, ny, nz;
697 int use_fac;
698
699 MRI_IMARR * ims_in;
700 MRI_IMAGE * im, *temp_im;
701
702
703 byte ** bptr = NULL ; /* one of these will be the array of */
704 short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
705 float ** fptr = NULL ; /* (depending on input datum type) */
706
707 float * fac = NULL ; /* array of brick scaling factors */
708
709 float * fout;
710
711
712 ntime = DSET_NUM_TIMES(dset) ;
713 nx = dset->daxes->nxx;
714 ny = dset->daxes->nyy;
715 nz = dset->daxes->nzz;
716 nvox = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
717 datum = DSET_BRICK_TYPE( dset , 0 ) ; /* get dataset datum type */
718
719 switch( datum ){ /* pointer type depends on input datum type */
720
721 default:
722 return NULL ;
723
724 /** create array of pointers into old dataset sub-bricks **/
725
726 /*--------- input is bytes ----------*/
727 /* voxel #i at time #k is bptr[k][i] */
728 /* for i=0..nvox-1 and k=0..ntime-1. */
729
730 case MRI_byte:
731 bptr = (byte **) malloc( sizeof(byte *) * ntime ) ;
732 if( bptr == NULL ) return NULL ;
733 for( kk=0 ; kk < ntime ; kk++ )
734 bptr[kk] = (byte *) DSET_ARRAY(dset,kk) ;
735 break ;
736
737 /*--------- input is shorts ---------*/
738 /* voxel #i at time #k is sptr[k][i] */
739 /* for i=0..nvox-1 and k=0..ntime-1. */
740
741 case MRI_short:
742 sptr = (short **) malloc( sizeof(short *) * ntime ) ;
743 if( sptr == NULL ) return NULL ;
744 for( kk=0 ; kk < ntime; kk++ )
745 sptr[kk] = (short *) DSET_ARRAY(dset,kk) ;
746 break ;
747
748 /*--------- input is floats ---------*/
749 /* voxel #i at time #k is fptr[k][i] */
750 /* for i=0..nvox-1 and k=0..ntime-1. */
751
752 case MRI_float:
753 fptr = (float **) malloc( sizeof(float *) * ntime) ;
754 if( fptr == NULL ) return NULL ;
755 for( kk=0 ; kk < ntime; kk++ )
756 fptr[kk] = (float *) DSET_ARRAY(dset,kk) ;
757 break ;
758
759 } /* end of switch on input type */
760
761 INIT_IMARR(ims_in) ;
762 for( kk=0 ; kk < ntime ; kk++ ){
763 im = mri_new_vol_empty( nx , ny , nz , datum ) ;
764 ADDTO_IMARR(ims_in,im) ;
765 }
766
767 for( kk=0 ; kk < ntime ; kk++ ){
768 im = IMARR_SUBIMAGE(ims_in,kk) ;
769
770 switch( datum ){
771 case MRI_byte: mri_fix_data_pointer( bptr[kk], im ) ; break ;
772 case MRI_short: mri_fix_data_pointer( sptr[kk], im ) ; break ;
773 case MRI_float: mri_fix_data_pointer( fptr[kk], im ) ; break ;
774 }
775 }
776
777
778
779 return(ims_in);
780 }
781
782
783