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 #ifndef __MAIN_PLUGIN_REORDER_C__
8 #define __MAIN_PLUGIN_REORDER_C__
9
10 /***************************************************
11 * File: plug_reorder.c
12 *
13 * Author: Jay Brian Kummer, Medical College of WI
14 * [kummer@its.mcw.edu]
15 *
16 * Date: May 8, 1997
17 *
18 * Purpose:
19 *
20 * Plugin to reorder the epochs (i.e., stimulus
21 * presentation regions) in a 3D+Time data set.
22 ****************************************************/
23
24 /***********************************************************************/
25 /* C O N D I T I O N A L C O M P I L A T I O N D I R E C T I V E S */
26 /***********************************************************************/
27
28 #define MAIN_PLUGIN_REORDER
29 /* Uncomment to make this plugin verbose:
30 #define DEBUG_PLUGIN_REORDER
31 */
32
33 /*********************************/
34 /* D A T A S T R U C T U R E S */
35 /*********************************/
36
37 typedef struct { /* returned from 'plug_reorder_parseMap' */
38 char classKRH; /* the class (i.e., [a-zA-Z]) of a time-course segment */
39 int length; /* the length of the class in time points */
40 } ClassInfo;
41
42 static char helpstring[] = /* Help string for plugin interface */
43 "Purpose: Reorders voxel time-courses within a 3D+time BRIK.\n"
44 "Input : Name of the 3D+time dataset to reorder.\n"
45 "\n"
46 "Output : Prefix of filename for the new dataset.\n"
47 "\n"
48 "Map file : Text file with instructions a reordering the input dataset.\n"
49 " Sample map file:\n"
50 " # Experiment: Bogus\n"
51 " # Date: April 1, 1997\n"
52 " -\n"
53 " -\n"
54 " D1\n"
55 " D2\n"
56 " D3\n"
57 " B1\n"
58 " B2\n"
59 " B3\n"
60 " B4\n"
61 " -\n"
62 " a1\n"
63 " a2\n"
64 " a3\n"
65 " A1\n"
66 " A2\n"
67 " A3\n"
68 " A1\n"
69 " A2\n"
70 " A3\n"
71 " -\n"
72 "\n"
73 " # : Comment line (ignored)\n"
74 " - : Point to be omitted from the output time-series. In the sample map,\n"
75 " time points 1, 2, 10, and 20 will be omitted.\n"
76 " Each non-comment line in a map file corresponds to a point in the time-series of\n"
77 " an input 3D+time dataset. Codes for each time point define the plug-in action for\n"
78 " those points. For example, the code in the fifth non-comment line of the sample\n"
79 " above, 'D3', means that the fifth point in the time-series is the third point in\n"
80 " the presentation of a stimulus named 'D'. \n"
81 " Stimulus codes are single alphabetic characters and are case-sensitive, with\n"
82 " uppercase preceding lowercase in the reordering (i.e., A-Z then a-z). Stimulus\n"
83 " regions are numbered starting from 1 (i.e., [ X0 X1 X2 ] and [ X1 X3 X5 ] will\n"
84 " cause the plug-in to abort), and, as shown in the sample with the code 'A', this\n"
85 " delimits duplicate instances of the same stimulus code. For the sample map, images\n"
86 " would be reordered so that the output time-series is:\n"
87 " [ A1 A2 A3 A1 A2 A3 B1 B2 B3 B4 C1 C2 C3 a1 a2 a3 ]\n"
88 "\n"
89 "Duplicates : If there are duplicates of a stimulus presentation within a map (e.g.,\n"
90 "[ A1 A2 A3 ] in the sample), they can be:\n"
91 " 1) Collated: Duplicate instances of a stimulus will be placed in the order in\n"
92 " which they occur in the map. In the case of the sample, this results in the\n"
93 " ordering [ A1 A2 A3 A1 A2 A3 ... ].\n"
94 " 2) Averaged: Duplicate instances of a stimulus will be averaged. In the sample,\n"
95 " this would result in the ordering [ A1' A2' A3' ... ], where Ai' is the mean\n"
96 " of the i-th point of all of the stimulus A's present in the map file.\n"
97 " **Duplicate stimulus regions must have the same length! For example,\n"
98 " [ D1 D2 X1 D1 D2 D3 ] will cause the plug-in to abort.\n"
99 ;
100
101 static char *dupaction_strings[] = /* strings for duplicate action in plugin interface */
102 { "Collate", "Average" };
103
104 #define NUM_DUPACTION_STRINGS (sizeof(dupaction_strings)/sizeof(char *))
105
106 /***************************/
107 /* H E A D E R F I L E S */
108 /***************************/
109
110 #include "afni.h"
111
112 #include <unistd.h>
113 #include "plug_reorder_parseMap.c"
114
115 /***********************************************************/
116 /* I N T E R N A L F U N C T I O N P R O T O T Y P E S */
117 /***********************************************************/
118
119 char *REORDER_main(PLUGIN_interface *);/* the entry point */
120
121 /*************************************/
122 /* S U P P O R T F U N C T I O N S */
123 /*************************************/
124
compare_class(const void * i,const void * j)125 int compare_class(const void *i, const void *j)
126 /* qsort compare routine for
127 sorting 'ClassInfo' array
128 */
129 {
130 ClassInfo *iTmp = (ClassInfo *)i;
131 ClassInfo *jTmp = (ClassInfo *)j;
132
133 if(iTmp->classKRH == jTmp->classKRH) {
134 return(0);
135 }
136 else if(iTmp->classKRH < jTmp->classKRH) {
137 return(-1);
138 }
139 else {
140 return(1);
141 }
142 }
143
144 /***********************************************************************
145 * Set up the interface to the user:
146 * 1) Create a new interface using "PLUTO_new_interface";
147 *
148 * 2) For each line of inputs, create the line with "PLUTO_add_option"
149 * (this line of inputs can be optional or mandatory);
150 *
151 * 3) For each item on the line, create the item with
152 * "PLUTO_add_dataset" for a dataset chooser,
153 * "PLUTO_add_string"for a string chooser,
154 * "PLUTO_add_number"for a number chooser.
155 ************************************************************************/
156
157
158 DEFINE_PLUGIN_PROTOTYPE
159
PLUGIN_init(int ncall)160 PLUGIN_interface *PLUGIN_init(int ncall)
161 {
162 PLUGIN_interface *plint; /* will be the output of this routine */
163
164 if(ncall >= 1) { /* only one interface */
165 return NULL;
166 }
167
168 CHECK_IF_ALLOWED("REORDER","Reorder") ; /* 30 Sep 2016 */
169
170 /*---------------- set titles and call point ----------------*/
171
172 plint = PLUTO_new_interface("Reorder",
173 "Reorders voxel time-courses within a 3D+Time BRIK",
174 helpstring,
175 PLUGIN_CALL_VIA_MENU,
176 REORDER_main);
177
178 PLUTO_set_sequence( plint , "z:Kummer" ) ;
179
180 /*--------- 1st line: Input dataset ---------*/
181
182 PLUTO_add_option(plint,
183 "Input",/* label at left of input line */
184 "Input",/* tag to return to plugin */
185 TRUE /* is this mandatory? */
186 );
187
188 PLUTO_add_dataset(plint,
189 "---->>", /* label next to button */
190 ANAT_ALL_MASK, /* take any anat datasets */
191 0, /* don't allow fim funcs*/
192 DIMEN_4D_MASK | /* need 3D+time datasets*/
193 BRICK_ALLTYPE_MASK /* handle any type */
194 );
195
196 /*---------- 2nd line: Output dataset prefix ----------*/
197
198 PLUTO_add_option(plint,
199 "Output", /* label at left of input line */
200 "Output", /* tag to return to plugin */
201 TRUE /* is this mandatory? */
202 );
203
204 PLUTO_add_string(plint,
205 "Prefix", /* label next to textfield */
206 0,NULL, /* no fixed strings to choose among */
207 19 /* 19 spaces for typing in value */
208 );
209
210 /*---------- 3rd line: Map file name ----------*/
211
212 PLUTO_add_option(plint,
213 "Map file", /* label at left of input line */
214 "MapFile", /* tag to return to plugin */
215 TRUE /* is this mandatory? */
216 );
217
218 PLUTO_add_string(plint,
219 "FileName", /* label next to textfield */
220 0,NULL, /* no fixed strings to choose among */
221 19 /* 19 spaces for typing in value */
222 );
223
224 PLUTO_add_string(plint,
225 "Duplicates", /* label next to textfield */
226 NUM_DUPACTION_STRINGS, /* number of strings to choose among */
227 dupaction_strings, /* fixed strings to choose among */
228 0 /* index of default (collate) */
229 );
230
231 /*--------- done with interface setup ---------*/
232
233 return plint;
234 }
235
236 /*************************************/
237 /* M A C R O D E F I N I T I O N S */
238 /*************************************/
239
240 /* Free vectors */
241 #undef FREEUP
242 #define FREEUP(thing) if((thing) != NULL) {free((thing)); (thing) = NULL;}
243
244 #define FREE_WORKSPACE do{ FREEUP(bPtr); FREEUP(sPtr); FREEUP(fPtr); \
245 FREEUP(map); FREEUP(classKRH); } while(0)
246
247 /*****************************************************************
248 * Main routine for this plugin (will be called from AFNI).
249 * If the return string is not NULL, some error transpired, and
250 * AFNI will popup the return string in a message box.
251 ******************************************************************/
252
REORDER_main(PLUGIN_interface * plint)253 char *REORDER_main(PLUGIN_interface *plint)
254 {
255 MCW_idcode *idc; /* input dataset idcode */
256 THD_3dim_dataset *old_dset; /* input dataset */
257 THD_3dim_dataset *new_dset; /* output dataset */
258
259 ClassInfo *classKRH; /* array of ClassInfo structures, returned by reference from 'parseMap' */
260
261 char *new_prefix; /* pointer to the output dataset prefix */
262 char *mapFile; /* pointer to the name of map file to open */
263 char *dupstr; /* pointer to duplicate action label string */
264 char currentClass; /* character sentinel */
265
266 register int ii; /* register counter variables */
267 register int jj;
268 register int kk;
269 register int mm;
270
271 int perc; /* percentage value for progress meter */
272 int ninp; /* number of time points */
273 int nvox; /* number of voxels */
274 int old_datum; /* data type of input dataset */
275 int error; /* general status variable */
276 int startIndex; /* index sentinel */
277 int currentLength; /* length of the current epoch */
278 int newLength; /* length of the time-course after averaging duplicates */
279 int instances; /* number of instances of the current epoch class */
280 int timePoint; /* current time point being processed */
281 int *map; /* array of indices for remapping time points, returned from 'parseMap' */
282 int mapLength; /* length of 'map' array, returned by reference from 'parseMap' */
283 int dupaction; /* 0: collate duplicate epochs, !0: average duplicate epochs */
284 int classCount; /* length of 'ClassInfo' array, returned by reference from 'parseMap' */
285
286 double sum; /* sum variable for averaging */
287 double divisor; /* divisor for calculating the mean, this will equal 'instances' */
288
289 byte **bPtr = NULL; /* for byte data, will point to the 3D+time data to be reordered */
290 short **sPtr = NULL; /* for short data, will point to the 3D+time data to be reordered */
291 float **fPtr = NULL; /* for float data, will point to the 3D+time data to be reordered */
292
293 byte **bData = NULL; /* for byte data, will point to the reordered 3D+time data */
294 short **sData = NULL; /* for short data, will point to the reordered 3D+time data */
295 float **fData = NULL; /* for float data, will point to the reordered 3D+time data */
296
297 Bool dupsToAverage = False; /* flag to indicate whether there are segments to average */
298
299 /**********************************************************************/
300 /****** Check inputs from AFNI to see if they are reasonable-ish ******/
301
302 /*--------- go to first input line ---------*/
303
304 PLUTO_next_option(plint);
305
306 idc = PLUTO_get_idcode(plint); /* get dataset item */
307
308 old_dset = PLUTO_find_dset(idc); /* get Ptr to dataset */
309 if(NULL == old_dset) {
310 return "*************************\n"
311 "Cannot find Input Dataset\n"
312 "*************************";
313 }
314
315 /*--------- go to second input line ---------*/
316
317 PLUTO_next_option(plint);
318
319 new_prefix = PLUTO_get_string(plint); /* get string item (the output prefix) */
320 if(! PLUTO_prefix_ok(new_prefix)) { /* check if it is OK */
321 return "************************\n"
322 "Output Prefix is illegal\n"
323 "************************";
324 }
325
326 /*--------- go to next input lines ---------*/
327
328 PLUTO_next_option(plint); /* skip to next line */
329
330 mapFile = PLUTO_get_string(plint); /* get the map filename */
331 if(NULL == mapFile || access(mapFile, R_OK) < 0) { /* check if it is OK and readable */
332 return "****************************\n"
333 "Epoch map file is unreadable\n"
334 "****************************";
335 }
336
337 dupstr = PLUTO_get_string(plint); /* get the string item */
338 dupaction = PLUTO_string_index(dupstr, /* find it in the list it came from */
339 NUM_DUPACTION_STRINGS,
340 dupaction_strings);
341
342 /********************************************************/
343 /*********** At this point, the inputs are OK ***********/
344
345 PLUTO_popup_meter(plint); /* popup a progress meter */
346
347 /*--------- set up pointers to each sub-brick in the input dataset ---------*/
348
349 #ifdef DEBUG_PLUGIN_REORDER
350 printf("[Reorder] Loading dataset into memory...\n");
351 #endif
352 DSET_load(old_dset); /* must be in memory before we get pointers to it */
353
354 nvox = /* number of voxels per sub-brick */
355 old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
356
357 #ifdef DEBUG_PLUGIN_REORDER
358 printf("[Reorder] %d voxels per subBRIK [%dx%dx%d].\n"
359 , nvox, old_dset->daxes->nxx, old_dset->daxes->nyy, old_dset->daxes->nzz);
360 #endif
361
362 ninp = DSET_NUM_TIMES(old_dset); /* get number of time points in old dataset */
363 #ifdef DEBUG_PLUGIN_REORDER
364 printf("[Reorder] %d time points [subBRIK].\n", ninp);
365 #endif
366
367 /********************************************************/
368 /*************** Parse the epoch map file ***************/
369
370 mapLength = ninp; /* initialize with the expected length for mapFile entries */
371 if(NULL == (map = REORDER_parseMap(mapFile
372 , &mapLength
373 , &classKRH
374 , &classCount))) {
375 return "*******************************\n"
376 "Critical error parsing map file\n"
377 "*******************************";
378 }
379
380 #ifdef DEBUG_PLUGIN_REORDER0
381 printf("\n[Reorder] Indices for epoch remapping:\n");
382 for(ii = 0; ii < mapLength; ii++) {
383 printf("%d\n", map[ii]);
384 }
385
386 printf("\n[Reorder] Meta-sequence of epoch classes is:\n");
387 for(ii = 0; ii < classCount; ii++) {
388 printf("[%d] Class %c [Width %d TRs]\n", ii, classKRH[ii].classKRH, classKRH[ii].length);
389 }
390 #endif
391
392 /* Sort 'class' to reflect the final order
393 to aid in any epoch averaging required */
394 qsort((void *)classKRH, classCount, sizeof(ClassInfo), compare_class);
395
396 #ifdef DEBUG_PLUGIN_REORDER
397 printf("\n[Reorder] Sorted meta-sequence of epoch classes is:\n");
398 for(ii = 0; ii < classCount; ii++) {
399 printf("[%d] Class %c [Width %d TRs]\n", ii, classKRH[ii].classKRH, classKRH[ii].length);
400 }
401 #endif
402
403 old_datum /* get old dataset datum type */
404 = DSET_BRICK_TYPE(old_dset, 0);
405
406 switch(old_datum) { /* pointer type depends on input datum type */
407
408 default:
409 return "*******************************\n"
410 "Illegal datum in Input Dataset\n"
411 "******************************";
412
413 /**create array of pointers into old dataset sub-bricks **/
414 /**Note that we skip the first 'ignore' sub-bricks here **/
415
416 /*---------- input is bytes -----------*/
417 /* voxel #i at time #k is bPtr[k][i] */
418 /* for i = 0..nvox-1 and k = 0..ninp-1.*/
419
420 case MRI_byte:
421 /* Workspace for 3D+time data */
422 bPtr = (byte **)calloc(sizeof(byte *), ninp);
423 if(NULL == bPtr) {
424 return "************************\n"
425 "!! Allocation Failure !!\n"
426 "************************";
427 }
428
429 for(kk = 0; kk < ninp; kk++) {
430 bPtr[kk] = (byte *)DSET_ARRAY(old_dset,kk);
431 }
432
433 /* Array of pointers for reordering subBRIKS */
434 bData = (byte **)calloc(sizeof(byte *), mapLength);
435 if(NULL == bData) {
436 FREEUP(bPtr);
437 return "************************\n"
438 "!! Allocation Failure !!\n"
439 "************************";
440 }
441 /* Allocate memory for subBRIKS */
442 for(kk = 0; kk < mapLength; kk++) {
443 bData[kk] = (byte *)calloc(sizeof(byte), nvox);
444 if(NULL == bData[kk]) {
445 FREEUP(bPtr);
446 FREEUP(bData);
447 return "************************\n"
448 "!! Allocation Failure !!\n"
449 "************************";
450 }
451 }
452 break;
453
454 /*---------- input is shorts ----------*/
455 /* voxel #i at time #k is sPtr[k][i] */
456 /* for i = 0..nvox-1 and k = 0..ninp-1.*/
457
458 case MRI_short:
459 /* Workspace for 3D+time data */
460 sPtr = (short **)calloc(sizeof(short *), ninp);
461 if(NULL == sPtr) {
462 return "************************\n"
463 "!! Allocation Failure !!\n"
464 "************************";
465 }
466
467 for(kk = 0; kk < ninp; kk++) {
468 sPtr[kk] = (short *)DSET_ARRAY(old_dset,kk);
469 }
470
471 /* Array of pointers for reordering subBRIKS */
472 sData = (short **)calloc(sizeof(short *), mapLength);
473 if(NULL == sData) {
474 FREEUP(sPtr);
475 return "************************\n"
476 "!! Allocation Failure !!\n"
477 "************************";
478 }
479 /* Allocate memory for subBRIKS */
480 for(kk = 0; kk < mapLength; kk++) {
481 sData[kk] = (short *)calloc(sizeof(short), nvox);
482 if(NULL == sData[kk]) {
483 FREEUP(sPtr);
484 FREEUP(sData);
485 return "************************\n"
486 "!! Allocation Failure !!\n"
487 "************************";
488 }
489 }
490 break;
491
492 /*---------- input is floats ----------*/
493 /* voxel #i at time #k is fPtr[k][i] */
494 /* for i = 0..nvox-1 and k = 0..ninp-1.*/
495
496 case MRI_float:
497 /* Workspace for 3D+time data */
498 fPtr = (float **)calloc(sizeof(float *), ninp);
499 if(NULL == fPtr) {
500 return "************************\n"
501 "!! Allocation Failure !!\n"
502 "************************";
503 }
504
505 for(kk = 0; kk < ninp; kk++) {
506 fPtr[kk] = (float *)DSET_ARRAY(old_dset,kk);
507 }
508
509 /* Array of pointers for reordering subBRIKS */
510 fData = (float **)calloc(sizeof(float *), mapLength);
511 if(NULL == fData) {
512 FREEUP(fPtr);
513 return "************************\n"
514 "!! Allocation Failure !!\n"
515 "************************";
516 }
517 /* Allocate memory for subBRIKS */
518 for(kk = 0; kk < mapLength; kk++) {
519 fData[kk] = (float *)calloc(sizeof(float), nvox);
520 if(NULL == fData[kk]) {
521 FREEUP(fPtr);
522 FREEUP(fData);
523 return "************************\n"
524 "!! Allocation Failure !!\n"
525 "************************";
526 }
527 }
528 break;
529
530 } /* end of switch on input type */
531
532 /*********************** Make a new dataset ***********************/
533
534 new_dset /* start with copy of old one */
535 = EDIT_empty_copy(old_dset);
536
537 { char * his = PLUTO_commandstring(plint) ;
538 tross_Copy_History(old_dset,new_dset) ;
539 tross_Append_History(new_dset,his) ; free(his);
540 }
541
542 /*-- edit some of its internal parameters --*/
543
544 ii = EDIT_dset_items(
545 new_dset,
546 ADN_prefix , new_prefix, /* filename prefix */
547 ADN_malloc_type, DATABLOCK_MEM_MALLOC, /* store in memory */
548 ADN_nvals , mapLength, /* # time points */
549 ADN_ntt , mapLength, /* # time points */
550 ADN_none);
551 if(0 != ii) {
552 THD_delete_3dim_dataset(new_dset, False);
553 switch(old_datum) {
554 case MRI_byte:
555 FREEUP(bPtr);
556 FREEUP(bData);
557 break;
558 case MRI_short:
559 FREEUP(sPtr);
560 FREEUP(sData);
561 break;
562 case MRI_float:
563 FREEUP(fPtr);
564 FREEUP(fData);
565 break;
566 }
567 FREE_WORKSPACE;
568 return "***********************************\n"
569 "Error while creating output dataset\n"
570 "***********************************";
571 }
572
573 /*****************************************************/
574 /****** Setup has ended. Now do some real work. ******/
575
576 switch(old_datum) {
577
578 /*****************************************/
579 /* voxel #i at time #k is array[kk][ii] */
580 /* for ii = 0..nvox-1 and kk = 0..ninp-1.*/
581
582 case MRI_byte:
583 #ifdef DEBUG_PLUGIN_REORDER
584 printf("\n[Reorder] Collating/byte data...\n");
585 #endif
586
587 /* Reorder the time-course of each voxel */
588 for(kk = 0; kk < mapLength; kk++) {
589 if(NULL == memcpy((void *)bData[kk]
590 , (void *)bPtr[map[kk]]
591 , sizeof(byte)*nvox)) {
592 THD_delete_3dim_dataset(new_dset, False);
593 FREEUP(bPtr);
594 FREEUP(bData);
595 FREE_WORKSPACE;
596 return "***********************************\n"
597 "!! Error performing memory copy !!\n"
598 "**********************************";
599 }
600 perc = (100 * kk) / mapLength; /* display percentage done */
601 PLUTO_set_meter(plint, perc); /* on the progress meter */
602 }
603
604 PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
605
606 if(dupaction) { /* new_dset may change in length again */
607 /* Verify that averaging is possible:
608 All duplicates must have the same length.
609 */
610 error = 0; /* start with no error */
611 for(ii = 0; ii < classCount; ii++) {
612 for(kk = 1; kk < (classCount-1); kk++) {
613 if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
614 if(classKRH[ii].length != classKRH[kk].length) {
615 error = 1; /* lengths differ! */
616 break;
617 }
618 dupsToAverage = True; /* there is at least one duplicate to average */
619 }
620 }
621 if(error) {
622 PLUTO_popup_message(plint
623 ,"*********************************************************\n"
624 "* Duplicate stimulus regions in map file must have the *\n"
625 "* same length for averaging. Returning collated regions.*\n"
626 "*********************************************************");
627 break;
628 }
629 }
630
631 if(!error) { /* average the crap */
632 #ifdef DEBUG_PLUGIN_REORDER
633 printf("[Reorder] All duplicate epochs have the same length...\n");
634 #endif
635
636 if(dupsToAverage) {
637 #ifdef DEBUG_PLUGIN_REORDER
638 printf("[Reorder] Averaging duplicate epochs...\n");
639 printf("[Reorder] Initial time-length is %d...\n", mapLength);
640 #endif
641
642 newLength = mapLength; /* length may change after averaging, track new length */
643
644 /* Perform averaging */
645 /* For each class, do: */
646 for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
647 #ifdef DEBUG_PLUGIN_REORDER
648 printf("[Reorder] Processing instances of class %c...\n"
649 , classKRH[ii].classKRH);
650 #endif
651
652 /* Get the number of duplicates of the current class */
653 for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
654 if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
655 break;
656 }
657 }
658
659 #ifdef DEBUG_PLUGIN_REORDER
660 printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
661 #endif
662
663 if(instances > 1) {
664 #ifdef DEBUG_PLUGIN_REORDER
665 printf(" Averaging %d contiguous epochs...\n", instances);
666 #endif
667
668 divisor = (double)instances;
669
670 /* For each time point in the duplicated classes: */
671 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
672 /* For each voxel: */
673 for(kk = 0; kk < nvox; kk++) {
674 /* Average the duplicate time-points for the current voxel */
675 for(mm = 0, sum = 0.0; mm < instances; mm++) {
676 sum = sum + (double)bData[jj+(mm*classKRH[ii].length)][kk];
677 }
678 bData[timePoint][kk] = (byte)((sum / divisor) + 0.5);
679 }
680 }
681
682 /* Update final length of time-course */
683 newLength = newLength - (classKRH[ii].length * (instances - 1));
684
685 /* Jump over the indices of all instances of the current class */
686 startIndex = startIndex + (instances *classKRH[ii].length);
687
688 /* Go to the next class to process */
689 ii += instances;
690 }
691 else { /* just copy the data points */
692 #ifdef DEBUG_PLUGIN_REORDER
693 printf(" Copying single epoch...\n");
694 #endif
695
696 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
697 for(kk = 0; kk < nvox; kk++) {
698 bData[timePoint][kk] = bData[jj][kk];
699 }
700 }
701
702 /* Jump to the index of the next class in the time-course */
703 startIndex = startIndex + classKRH[ii].length;
704
705 /* eat up another class */
706 ++ii;
707 }
708
709 /* Update progress meter with the percentage done */
710 perc = (100 * timePoint) / mapLength;
711 PLUTO_set_meter(plint, perc);
712 }
713
714 /* Set progress meter to 100% */
715 PLUTO_set_meter(plint, 100);
716
717 #ifdef DEBUG_PLUGIN_REORDER
718 printf("[Reorder] Final time-length is %d...\n", newLength);
719 #endif
720 ii = EDIT_dset_items(
721 new_dset,
722 ADN_nvals, newLength, /* # time points */
723 ADN_ntt, newLength, /* # time points */
724 ADN_none);
725 if(0 != ii) {
726 THD_delete_3dim_dataset(new_dset, False);
727 FREEUP(bPtr);
728 FREEUP(bData);
729 FREE_WORKSPACE;
730 return "******************************************\n"
731 "!! Error while creating output dataset !!\n"
732 "*****************************************";
733 }
734
735 /* Free unused subBRIKS */
736 if(newLength < mapLength) {
737 #ifdef DEBUG_PLUGIN_REORDER
738 printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
739 #endif
740 for(ii = newLength; ii < mapLength; ii++) {
741 FREEUP(bData[ii]);
742 }
743 bData = (byte **)realloc((void *)bData, sizeof(byte *)*newLength);
744 if(NULL == bData) {
745 THD_delete_3dim_dataset(new_dset, False);
746 FREEUP(bPtr);
747 FREE_WORKSPACE;
748 return "************************\n"
749 "!! Reallocation error !!\n"
750 "************************";
751 }
752 }
753 mapLength = newLength;
754 }
755 }
756 }
757
758 /* Write results into new dataset */
759 #ifdef DEBUG_PLUGIN_REORDER
760 printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
761 #endif
762 for(kk = 0; kk < mapLength; kk++) {
763 EDIT_substitute_brick(new_dset, kk, MRI_byte, bData[kk]);
764 }
765 break;
766
767 case MRI_short:
768 #ifdef DEBUG_PLUGIN_REORDER
769 printf("\n[Reorder] Collating/short data...\n");
770 #endif
771
772 /* Reorder the time-course of each voxel */
773 for(kk = 0; kk < mapLength; kk++) {
774 if(NULL == memcpy((void *)sData[kk]
775 , (void *)sPtr[map[kk]]
776 , sizeof(short)*nvox)) {
777 THD_delete_3dim_dataset(new_dset, False);
778 FREEUP(sPtr);
779 FREEUP(sData);
780 FREE_WORKSPACE;
781 return "***********************************\n"
782 "!! Error performing memory copy !!\n"
783 "**********************************";
784 }
785 perc = (100 * kk) / mapLength; /* display percentage done */
786 PLUTO_set_meter(plint, perc); /* on the progress meter */
787 }
788
789 PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
790
791 if(dupaction) { /* new_dset may change in length again */
792 /* Verify that averaging is possible:
793 All duplicates must have the same length.
794 */
795 error = 0; /* start with no error */
796 for(ii = 0; ii < classCount; ii++) {
797 for(kk = 1; kk < (classCount-1); kk++) {
798 if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
799 if(classKRH[ii].length != classKRH[kk].length) {
800 error = 1; /* lengths differ! */
801 break;
802 }
803 dupsToAverage = True; /* there is at least one duplicate to average */
804 }
805 }
806 if(error) {
807 PLUTO_popup_message(plint
808 ,"*********************************************************\n"
809 "* Duplicate stimulus regions in map file must have the *\n"
810 "* same length for averaging. Returning collated regions.*\n"
811 "*********************************************************");
812 break;
813 }
814 }
815
816 if(!error) { /* average the crap */
817 #ifdef DEBUG_PLUGIN_REORDER
818 printf("[Reorder] All duplicate epochs have the same length...\n");
819 #endif
820
821 if(dupsToAverage) {
822 #ifdef DEBUG_PLUGIN_REORDER
823 printf("[Reorder] Averaging duplicate epochs...\n");
824 printf("[Reorder] Initial time-length is %d...\n", mapLength);
825 #endif
826
827 newLength = mapLength; /* length may change after averaging, track new length */
828
829 /* Perform averaging */
830 /* For each class, do: */
831 for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
832 #ifdef DEBUG_PLUGIN_REORDER
833 printf("[Reorder] Processing instances of class %c...\n"
834 , classKRH[ii].classKRH);
835 #endif
836
837 /* Get the number of duplicates of the current class */
838 for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
839 if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
840 break;
841 }
842 }
843
844 #ifdef DEBUG_PLUGIN_REORDER
845 printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
846 #endif
847
848 if(instances > 1) {
849 #ifdef DEBUG_PLUGIN_REORDER
850 printf(" Averaging %d contiguous epochs...\n", instances);
851 #endif
852
853 divisor = (double)instances;
854
855 /* For each time point in the duplicated classes: */
856 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
857 /* For each voxel: */
858 for(kk = 0; kk < nvox; kk++) {
859 /* Average the duplicate time-points for the current voxel */
860 for(mm = 0, sum = 0.0; mm < instances; mm++) {
861 sum = sum + (double)sData[jj+(mm*classKRH[ii].length)][kk];
862 }
863 sData[timePoint][kk] = (short)((sum / divisor) + 0.5);
864 }
865 }
866
867 /* Update final length of time-course */
868 newLength = newLength - (classKRH[ii].length * (instances - 1));
869
870 /* Jump over the indices of all instances of the current class */
871 startIndex = startIndex + (instances *classKRH[ii].length);
872
873 /* Go to the next class to process */
874 ii += instances;
875 }
876 else { /* just copy the data points */
877 #ifdef DEBUG_PLUGIN_REORDER
878 printf(" Copying single epoch...\n");
879 #endif
880
881 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
882 for(kk = 0; kk < nvox; kk++) {
883 sData[timePoint][kk] = sData[jj][kk];
884 }
885 }
886
887 /* Jump to the index of the next class in the time-course */
888 startIndex = startIndex + classKRH[ii].length;
889
890 /* eat up another class */
891 ++ii;
892 }
893
894 /* Update progress meter with the percentage done */
895 perc = (100 * timePoint) / mapLength;
896 PLUTO_set_meter(plint, perc);
897 }
898
899 /* Set progress meter to 100% */
900 PLUTO_set_meter(plint, 100);
901
902 #ifdef DEBUG_PLUGIN_REORDER
903 printf("[Reorder] Final time-length is %d...\n", newLength);
904 #endif
905 ii = EDIT_dset_items(
906 new_dset,
907 ADN_nvals, newLength, /* # time points */
908 ADN_ntt, newLength, /* # time points */
909 ADN_none);
910 if(0 != ii) {
911 THD_delete_3dim_dataset(new_dset, False);
912 FREEUP(sPtr);
913 FREEUP(sData);
914 FREE_WORKSPACE;
915 return "******************************************\n"
916 "!! Error while creating output dataset !!\n"
917 "*****************************************";
918 }
919
920 /* Free unused subBRIKS */
921 if(newLength < mapLength) {
922 #ifdef DEBUG_PLUGIN_REORDER
923 printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
924 #endif
925 for(ii = newLength; ii < mapLength; ii++) {
926 FREEUP(sData[ii]);
927 }
928 sData = (short **)realloc((void *)sData, sizeof(short *)*newLength);
929 if(NULL == sData) {
930 THD_delete_3dim_dataset(new_dset, False);
931 FREEUP(sPtr);
932 FREE_WORKSPACE;
933 return "*************************\n"
934 "!! Reallocation error !!\n"
935 "************************";
936 }
937 }
938 mapLength = newLength;
939 }
940 }
941 }
942
943 /* Write results into new dataset */
944 #ifdef DEBUG_PLUGIN_REORDER
945 printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
946 #endif
947 for(kk = 0; kk < mapLength; kk++) {
948 EDIT_substitute_brick(new_dset, kk, MRI_short, sData[kk]);
949 }
950 break;
951
952 case MRI_float:
953 #ifdef DEBUG_PLUGIN_REORDER
954 printf("\n[Reorder] Collating/float data...\n");
955 #endif
956
957 /* Reorder the time-course of each voxel */
958 for(kk = 0; kk < mapLength; kk++) {
959 if(NULL == memcpy((void *)fData[kk]
960 , (void *)fPtr[map[kk]]
961 , sizeof(float)*nvox)) {
962 THD_delete_3dim_dataset(new_dset, False);
963 FREEUP(fPtr);
964 FREEUP(fData);
965 FREE_WORKSPACE;
966 return "***********************************\n"
967 "!! Error performing memory copy !!\n"
968 "**********************************";
969 }
970 perc = (100 * kk) / mapLength; /* display percentage done */
971 PLUTO_set_meter(plint, perc); /* on the progress meter */
972 }
973
974 PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
975
976 if(dupaction) { /* new_dset may change in length again */
977 /* Verify that averaging is possible:
978 All duplicates must have the same length.
979 */
980 error = 0; /* start with no error */
981 for(ii = 0; ii < classCount; ii++) {
982 for(kk = 1; kk < (classCount-1); kk++) {
983 if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
984 if(classKRH[ii].length != classKRH[kk].length) {
985 error = 1; /* lengths differ! */
986 break;
987 }
988 dupsToAverage = True; /* there is at least one duplicate to average */
989 }
990 }
991 if(error) {
992 PLUTO_popup_message(plint
993 ,"*********************************************************\n"
994 "* Duplicate stimulus regions in map file must have the *\n"
995 "* same length for averaging. Returning collated regions.*\n"
996 "*********************************************************");
997 break;
998 }
999 }
1000
1001 if(!error) { /* average the crap */
1002 #ifdef DEBUG_PLUGIN_REORDER
1003 printf("[Reorder] All duplicate epochs have the same length...\n");
1004 #endif
1005
1006 if(dupsToAverage) {
1007 #ifdef DEBUG_PLUGIN_REORDER
1008 printf("[Reorder] Averaging duplicate epochs...\n");
1009 printf("[Reorder] Initial time-length is %d...\n", mapLength);
1010 #endif
1011
1012 newLength = mapLength; /* length may change after averaging, track new length */
1013
1014 /* Perform averaging */
1015 /* For each class, do: */
1016 for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
1017 #ifdef DEBUG_PLUGIN_REORDER
1018 printf("[Reorder] Processing instances of class %c...\n"
1019 , classKRH[ii].classKRH);
1020 #endif
1021
1022 /* Get the number of duplicates of the current class */
1023 for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
1024 if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
1025 break;
1026 }
1027 }
1028
1029 #ifdef DEBUG_PLUGIN_REORDER
1030 printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
1031 #endif
1032
1033 if(instances > 1) {
1034 #ifdef DEBUG_PLUGIN_REORDER
1035 printf(" Averaging %d contiguous epochs...\n", instances);
1036 #endif
1037
1038 divisor = (double)instances;
1039
1040 /* For each time point in the duplicated classes: */
1041 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
1042 /* For each voxel: */
1043 for(kk = 0; kk < nvox; kk++) {
1044 /* Average the duplicate time-points for the current voxel */
1045 for(mm = 0, sum = 0.0; mm < instances; mm++) {
1046 sum = sum + (double)fData[jj+(mm*classKRH[ii].length)][kk];
1047 }
1048 fData[timePoint][kk] = (float)((sum / divisor) + 0.5);
1049 }
1050 }
1051
1052 /* Update final length of time-course */
1053 newLength = newLength - (classKRH[ii].length * (instances - 1));
1054
1055 /* Jump over the indices of all instances of the current class */
1056 startIndex = startIndex + (instances *classKRH[ii].length);
1057
1058 /* Go to the next class to process */
1059 ii += instances;
1060 }
1061 else { /* just copy the data points */
1062 #ifdef DEBUG_PLUGIN_REORDER
1063 printf(" Copying single epoch...\n");
1064 #endif
1065
1066 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
1067 for(kk = 0; kk < nvox; kk++) {
1068 fData[timePoint][kk] = fData[jj][kk];
1069 }
1070 }
1071
1072 /* Jump to the index of the next class in the time-course */
1073 startIndex = startIndex + classKRH[ii].length;
1074
1075 /* eat up another class */
1076 ++ii;
1077 }
1078
1079 /* Update progress meter with the percentage done */
1080 perc = (100 * timePoint) / mapLength;
1081 PLUTO_set_meter(plint, perc);
1082 }
1083
1084 /* Set progress meter to 100% */
1085 PLUTO_set_meter(plint, 100);
1086
1087 #ifdef DEBUG_PLUGIN_REORDER
1088 printf("[Reorder] Final time-length is %d...\n", newLength);
1089 #endif
1090 ii = EDIT_dset_items(
1091 new_dset,
1092 ADN_nvals, newLength, /* # time points */
1093 ADN_ntt, newLength, /* # time points */
1094 ADN_none);
1095 if(0 != ii) {
1096 THD_delete_3dim_dataset(new_dset, False);
1097 FREEUP(fPtr);
1098 FREEUP(fData);
1099 FREE_WORKSPACE;
1100 return "******************************************\n"
1101 "!! Error while creating output dataset !!\n"
1102 "*****************************************";
1103 }
1104
1105 /* Free unused subBRIKS */
1106 if(newLength < mapLength) {
1107 #ifdef DEBUG_PLUGIN_REORDER
1108 printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
1109 #endif
1110 for(ii = newLength; ii < mapLength; ii++) {
1111 FREEUP(fData[ii]);
1112 }
1113 fData = (float **)realloc((void *)fData, sizeof(float *)*newLength);
1114 if(NULL == fData) {
1115 THD_delete_3dim_dataset(new_dset, False);
1116 FREEUP(fPtr);
1117 FREE_WORKSPACE;
1118 return "*************************\n"
1119 "!! Reallocation error !!\n"
1120 "************************";
1121 }
1122 }
1123 mapLength = newLength;
1124 }
1125 }
1126 }
1127
1128 /* Write results into new dataset */
1129 #ifdef DEBUG_PLUGIN_REORDER
1130 printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
1131 #endif
1132 for(kk = 0; kk < mapLength; kk++) {
1133 EDIT_substitute_brick(new_dset, kk, MRI_float, fData[kk]);
1134 }
1135 break;
1136
1137 } /* end of switch over input type */
1138
1139 DSET_unload(old_dset);
1140
1141 /*************** Cleanup and go home *****************/
1142
1143 PLUTO_add_dset(plint, new_dset, DSET_ACTION_MAKE_CURRENT);
1144
1145 #ifdef DEBUG_PLUGIN_REORDER
1146 printf("[Reorder] Cleaning up workspace...\n");
1147 #endif
1148 FREE_WORKSPACE;
1149
1150 return NULL; /* NULL string returned means all was OK */
1151 }
1152
1153 #endif
1154
1155