Doxygen Source Code Documentation
plug_reorder.c File Reference
#include "afni.h"#include <unistd.h>#include "plug_reorder_parseMap.c"Go to the source code of this file.
Data Structures | |
| struct | ClassInfo |
Defines | |
| #define | MAIN_PLUGIN_REORDER |
| #define | NUM_DUPACTION_STRINGS (sizeof(dupaction_strings)/sizeof(char *)) |
| #define | FREEUP(thing) if((thing) != NULL) {free((thing)); (thing) = NULL;} |
| #define | FREE_WORKSPACE |
Functions | |
| char * | REORDER_main (PLUGIN_interface *) |
| int | compare_class (const void *i, const void *j) |
| DEFINE_PLUGIN_PROTOTYPE PLUGIN_interface * | PLUGIN_init (int ncall) |
Variables | |
| char | helpstring [] |
| char * | dupaction_strings [] |
Define Documentation
|
|
Value: Definition at line 242 of file plug_reorder.c. Referenced by REORDER_main(). |
|
|
Definition at line 240 of file plug_reorder.c. |
|
|
Definition at line 28 of file plug_reorder.c. |
|
|
Definition at line 104 of file plug_reorder.c. Referenced by PLUGIN_init(), and REORDER_main(). |
Function Documentation
|
||||||||||||
|
Definition at line 125 of file plug_reorder.c. References ClassInfo::classKRH, and i. Referenced by REORDER_main().
|
|
|
Definition at line 160 of file plug_reorder.c. References ANAT_ALL_MASK, dupaction_strings, helpstring, NUM_DUPACTION_STRINGS, PLUTO_set_sequence(), and REORDER_main().
00161 {
00162 PLUGIN_interface *plint; /* will be the output of this routine */
00163
00164 if(ncall >= 1) { /* only one interface */
00165 return NULL;
00166 }
00167
00168 /*---------------- set titles and call point ----------------*/
00169
00170 plint = PLUTO_new_interface("Reorder",
00171 "Reorders voxel time-courses within a 3D+Time BRIK",
00172 helpstring,
00173 PLUGIN_CALL_VIA_MENU,
00174 REORDER_main);
00175
00176 PLUTO_set_sequence( plint , "z:Kummer" ) ;
00177
00178 /*--------- 1st line: Input dataset ---------*/
00179
00180 PLUTO_add_option(plint,
00181 "Input",/* label at left of input line */
00182 "Input",/* tag to return to plugin */
00183 TRUE /* is this mandatory? */
00184 );
00185
00186 PLUTO_add_dataset(plint,
00187 "---->>", /* label next to button */
00188 ANAT_ALL_MASK, /* take any anat datasets */
00189 0, /* don't allow fim funcs*/
00190 DIMEN_4D_MASK | /* need 3D+time datasets*/
00191 BRICK_ALLTYPE_MASK /* handle any type */
00192 );
00193
00194 /*---------- 2nd line: Output dataset prefix ----------*/
00195
00196 PLUTO_add_option(plint,
00197 "Output", /* label at left of input line */
00198 "Output", /* tag to return to plugin */
00199 TRUE /* is this mandatory? */
00200 );
00201
00202 PLUTO_add_string(plint,
00203 "Prefix", /* label next to textfield */
00204 0,NULL, /* no fixed strings to choose among */
00205 19 /* 19 spaces for typing in value */
00206 );
00207
00208 /*---------- 3rd line: Map file name ----------*/
00209
00210 PLUTO_add_option(plint,
00211 "Map file", /* label at left of input line */
00212 "MapFile", /* tag to return to plugin */
00213 TRUE /* is this mandatory? */
00214 );
00215
00216 PLUTO_add_string(plint,
00217 "FileName", /* label next to textfield */
00218 0,NULL, /* no fixed strings to choose among */
00219 19 /* 19 spaces for typing in value */
00220 );
00221
00222 PLUTO_add_string(plint,
00223 "Duplicates", /* label next to textfield */
00224 NUM_DUPACTION_STRINGS, /* number of strings to choose among */
00225 dupaction_strings, /* fixed strings to choose among */
00226 0 /* index of default (collate) */
00227 );
00228
00229 /*--------- done with interface setup ---------*/
00230
00231 return plint;
00232 }
|
|
|
Definition at line 251 of file plug_reorder.c. References ADN_malloc_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, Bool, calloc, ClassInfo::classKRH, compare_class(), DATABLOCK_MEM_MALLOC, THD_3dim_dataset::daxes, DSET_ARRAY, DSET_BRICK_TYPE, DSET_load, DSET_NUM_TIMES, DSET_unload, dupaction_strings, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), free, FREE_WORKSPACE, FREEUP, ClassInfo::length, NUM_DUPACTION_STRINGS, THD_dataxes::nxx, THD_dataxes::nyy, THD_dataxes::nzz, PLUTO_add_dset(), PLUTO_commandstring(), PLUTO_find_dset(), PLUTO_popup_meter(), PLUTO_prefix_ok(), PLUTO_set_meter(), PLUTO_string_index(), realloc, REORDER_parseMap(), THD_delete_3dim_dataset(), tross_Append_History(), and tross_Copy_History(). Referenced by PLUGIN_init().
00252 {
00253 MCW_idcode *idc; /* input dataset idcode */
00254 THD_3dim_dataset *old_dset; /* input dataset */
00255 THD_3dim_dataset *new_dset; /* output dataset */
00256
00257 ClassInfo *classKRH; /* array of ClassInfo structures, returned by reference from 'parseMap' */
00258
00259 char *new_prefix; /* pointer to the output dataset prefix */
00260 char *mapFile; /* pointer to the name of map file to open */
00261 char *dupstr; /* pointer to duplicate action label string */
00262 char currentClass; /* character sentinel */
00263
00264 register int ii; /* register counter variables */
00265 register int jj;
00266 register int kk;
00267 register int mm;
00268
00269 int perc; /* percentage value for progress meter */
00270 int ninp; /* number of time points */
00271 int nvox; /* number of voxels */
00272 int old_datum; /* data type of input dataset */
00273 int error; /* general status variable */
00274 int startIndex; /* index sentinel */
00275 int currentLength; /* length of the current epoch */
00276 int newLength; /* length of the time-course after averaging duplicates */
00277 int instances; /* number of instances of the current epoch class */
00278 int timePoint; /* current time point being processed */
00279 int *map; /* array of indices for remapping time points, returned from 'parseMap' */
00280 int mapLength; /* length of 'map' array, returned by reference from 'parseMap' */
00281 int dupaction; /* 0: collate duplicate epochs, !0: average duplicate epochs */
00282 int classCount; /* length of 'ClassInfo' array, returned by reference from 'parseMap' */
00283
00284 double sum; /* sum variable for averaging */
00285 double divisor; /* divisor for calculating the mean, this will equal 'instances' */
00286
00287 byte **bPtr = NULL; /* for byte data, will point to the 3D+time data to be reordered */
00288 short **sPtr = NULL; /* for short data, will point to the 3D+time data to be reordered */
00289 float **fPtr = NULL; /* for float data, will point to the 3D+time data to be reordered */
00290
00291 byte **bData = NULL; /* for byte data, will point to the reordered 3D+time data */
00292 short **sData = NULL; /* for short data, will point to the reordered 3D+time data */
00293 float **fData = NULL; /* for float data, will point to the reordered 3D+time data */
00294
00295 Bool dupsToAverage = False; /* flag to indicate whether there are segments to average */
00296
00297 /**********************************************************************/
00298 /****** Check inputs from AFNI to see if they are reasonable-ish ******/
00299
00300 /*--------- go to first input line ---------*/
00301
00302 PLUTO_next_option(plint);
00303
00304 idc = PLUTO_get_idcode(plint); /* get dataset item */
00305
00306 old_dset = PLUTO_find_dset(idc); /* get Ptr to dataset */
00307 if(NULL == old_dset) {
00308 return "*************************\n"
00309 "Cannot find Input Dataset\n"
00310 "*************************";
00311 }
00312
00313 /*--------- go to second input line ---------*/
00314
00315 PLUTO_next_option(plint);
00316
00317 new_prefix = PLUTO_get_string(plint); /* get string item (the output prefix) */
00318 if(! PLUTO_prefix_ok(new_prefix)) { /* check if it is OK */
00319 return "************************\n"
00320 "Output Prefix is illegal\n"
00321 "************************";
00322 }
00323
00324 /*--------- go to next input lines ---------*/
00325
00326 PLUTO_next_option(plint); /* skip to next line */
00327
00328 mapFile = PLUTO_get_string(plint); /* get the map filename */
00329 if(NULL == mapFile || access(mapFile, R_OK) < 0) { /* check if it is OK and readable */
00330 return "****************************\n"
00331 "Epoch map file is unreadable\n"
00332 "****************************";
00333 }
00334
00335 dupstr = PLUTO_get_string(plint); /* get the string item */
00336 dupaction = PLUTO_string_index(dupstr, /* find it in the list it came from */
00337 NUM_DUPACTION_STRINGS,
00338 dupaction_strings);
00339
00340 /********************************************************/
00341 /*********** At this point, the inputs are OK ***********/
00342
00343 PLUTO_popup_meter(plint); /* popup a progress meter */
00344
00345 /*--------- set up pointers to each sub-brick in the input dataset ---------*/
00346
00347 #ifdef DEBUG_PLUGIN_REORDER
00348 printf("[Reorder] Loading dataset into memory...\n");
00349 #endif
00350 DSET_load(old_dset); /* must be in memory before we get pointers to it */
00351
00352 nvox = /* number of voxels per sub-brick */
00353 old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
00354
00355 #ifdef DEBUG_PLUGIN_REORDER
00356 printf("[Reorder] %d voxels per subBRIK [%dx%dx%d].\n"
00357 , nvox, old_dset->daxes->nxx, old_dset->daxes->nyy, old_dset->daxes->nzz);
00358 #endif
00359
00360 ninp = DSET_NUM_TIMES(old_dset); /* get number of time points in old dataset */
00361 #ifdef DEBUG_PLUGIN_REORDER
00362 printf("[Reorder] %d time points [subBRIK].\n", ninp);
00363 #endif
00364
00365 /********************************************************/
00366 /*************** Parse the epoch map file ***************/
00367
00368 mapLength = ninp; /* initialize with the expected length for mapFile entries */
00369 if(NULL == (map = REORDER_parseMap(mapFile
00370 , &mapLength
00371 , &classKRH
00372 , &classCount))) {
00373 return "*******************************\n"
00374 "Critical error parsing map file\n"
00375 "*******************************";
00376 }
00377
00378 #ifdef DEBUG_PLUGIN_REORDER0
00379 printf("\n[Reorder] Indices for epoch remapping:\n");
00380 for(ii = 0; ii < mapLength; ii++) {
00381 printf("%d\n", map[ii]);
00382 }
00383
00384 printf("\n[Reorder] Meta-sequence of epoch classes is:\n");
00385 for(ii = 0; ii < classCount; ii++) {
00386 printf("[%d] Class %c [Width %d TRs]\n", ii, classKRH[ii].classKRH, classKRH[ii].length);
00387 }
00388 #endif
00389
00390 /* Sort 'class' to reflect the final order
00391 to aid in any epoch averaging required */
00392 qsort((void *)classKRH, classCount, sizeof(ClassInfo), compare_class);
00393
00394 #ifdef DEBUG_PLUGIN_REORDER
00395 printf("\n[Reorder] Sorted meta-sequence of epoch classes is:\n");
00396 for(ii = 0; ii < classCount; ii++) {
00397 printf("[%d] Class %c [Width %d TRs]\n", ii, classKRH[ii].classKRH, classKRH[ii].length);
00398 }
00399 #endif
00400
00401 old_datum /* get old dataset datum type */
00402 = DSET_BRICK_TYPE(old_dset, 0);
00403
00404 switch(old_datum) { /* pointer type depends on input datum type */
00405
00406 default:
00407 return "*******************************\n"
00408 "Illegal datum in Input Dataset\n"
00409 "******************************";
00410
00411 /**create array of pointers into old dataset sub-bricks **/
00412 /**Note that we skip the first 'ignore' sub-bricks here **/
00413
00414 /*---------- input is bytes -----------*/
00415 /* voxel #i at time #k is bPtr[k][i] */
00416 /* for i = 0..nvox-1 and k = 0..ninp-1.*/
00417
00418 case MRI_byte:
00419 /* Workspace for 3D+time data */
00420 bPtr = (byte **)calloc(sizeof(byte *), ninp);
00421 if(NULL == bPtr) {
00422 return "************************\n"
00423 "!! Allocation Failure !!\n"
00424 "************************";
00425 }
00426
00427 for(kk = 0; kk < ninp; kk++) {
00428 bPtr[kk] = (byte *)DSET_ARRAY(old_dset,kk);
00429 }
00430
00431 /* Array of pointers for reordering subBRIKS */
00432 bData = (byte **)calloc(sizeof(byte *), mapLength);
00433 if(NULL == bData) {
00434 FREEUP(bPtr);
00435 return "************************\n"
00436 "!! Allocation Failure !!\n"
00437 "************************";
00438 }
00439 /* Allocate memory for subBRIKS */
00440 for(kk = 0; kk < mapLength; kk++) {
00441 bData[kk] = (byte *)calloc(sizeof(byte), nvox);
00442 if(NULL == bData[kk]) {
00443 FREEUP(bPtr);
00444 FREEUP(bData);
00445 return "************************\n"
00446 "!! Allocation Failure !!\n"
00447 "************************";
00448 }
00449 }
00450 break;
00451
00452 /*---------- input is shorts ----------*/
00453 /* voxel #i at time #k is sPtr[k][i] */
00454 /* for i = 0..nvox-1 and k = 0..ninp-1.*/
00455
00456 case MRI_short:
00457 /* Workspace for 3D+time data */
00458 sPtr = (short **)calloc(sizeof(short *), ninp);
00459 if(NULL == sPtr) {
00460 return "************************\n"
00461 "!! Allocation Failure !!\n"
00462 "************************";
00463 }
00464
00465 for(kk = 0; kk < ninp; kk++) {
00466 sPtr[kk] = (short *)DSET_ARRAY(old_dset,kk);
00467 }
00468
00469 /* Array of pointers for reordering subBRIKS */
00470 sData = (short **)calloc(sizeof(short *), mapLength);
00471 if(NULL == sData) {
00472 FREEUP(sPtr);
00473 return "************************\n"
00474 "!! Allocation Failure !!\n"
00475 "************************";
00476 }
00477 /* Allocate memory for subBRIKS */
00478 for(kk = 0; kk < mapLength; kk++) {
00479 sData[kk] = (short *)calloc(sizeof(short), nvox);
00480 if(NULL == sData[kk]) {
00481 FREEUP(sPtr);
00482 FREEUP(sData);
00483 return "************************\n"
00484 "!! Allocation Failure !!\n"
00485 "************************";
00486 }
00487 }
00488 break;
00489
00490 /*---------- input is floats ----------*/
00491 /* voxel #i at time #k is fPtr[k][i] */
00492 /* for i = 0..nvox-1 and k = 0..ninp-1.*/
00493
00494 case MRI_float:
00495 /* Workspace for 3D+time data */
00496 fPtr = (float **)calloc(sizeof(float *), ninp);
00497 if(NULL == fPtr) {
00498 return "************************\n"
00499 "!! Allocation Failure !!\n"
00500 "************************";
00501 }
00502
00503 for(kk = 0; kk < ninp; kk++) {
00504 fPtr[kk] = (float *)DSET_ARRAY(old_dset,kk);
00505 }
00506
00507 /* Array of pointers for reordering subBRIKS */
00508 fData = (float **)calloc(sizeof(float *), mapLength);
00509 if(NULL == fData) {
00510 FREEUP(fPtr);
00511 return "************************\n"
00512 "!! Allocation Failure !!\n"
00513 "************************";
00514 }
00515 /* Allocate memory for subBRIKS */
00516 for(kk = 0; kk < mapLength; kk++) {
00517 fData[kk] = (float *)calloc(sizeof(float), nvox);
00518 if(NULL == fData[kk]) {
00519 FREEUP(fPtr);
00520 FREEUP(fData);
00521 return "************************\n"
00522 "!! Allocation Failure !!\n"
00523 "************************";
00524 }
00525 }
00526 break;
00527
00528 } /* end of switch on input type */
00529
00530 /*********************** Make a new dataset ***********************/
00531
00532 new_dset /* start with copy of old one */
00533 = EDIT_empty_copy(old_dset);
00534
00535 { char * his = PLUTO_commandstring(plint) ;
00536 tross_Copy_History(old_dset,new_dset) ;
00537 tross_Append_History(new_dset,his) ; free(his);
00538 }
00539
00540 /*-- edit some of its internal parameters --*/
00541
00542 ii = EDIT_dset_items(
00543 new_dset,
00544 ADN_prefix , new_prefix, /* filename prefix */
00545 ADN_malloc_type, DATABLOCK_MEM_MALLOC, /* store in memory */
00546 ADN_nvals , mapLength, /* # time points */
00547 ADN_ntt , mapLength, /* # time points */
00548 ADN_none);
00549 if(0 != ii) {
00550 THD_delete_3dim_dataset(new_dset, False);
00551 switch(old_datum) {
00552 case MRI_byte:
00553 FREEUP(bPtr);
00554 FREEUP(bData);
00555 break;
00556 case MRI_short:
00557 FREEUP(sPtr);
00558 FREEUP(sData);
00559 break;
00560 case MRI_float:
00561 FREEUP(fPtr);
00562 FREEUP(fData);
00563 break;
00564 }
00565 FREE_WORKSPACE;
00566 return "***********************************\n"
00567 "Error while creating output dataset\n"
00568 "***********************************";
00569 }
00570
00571 /*****************************************************/
00572 /****** Setup has ended. Now do some real work. ******/
00573
00574 switch(old_datum) {
00575
00576 /*****************************************/
00577 /* voxel #i at time #k is array[kk][ii] */
00578 /* for ii = 0..nvox-1 and kk = 0..ninp-1.*/
00579
00580 case MRI_byte:
00581 #ifdef DEBUG_PLUGIN_REORDER
00582 printf("\n[Reorder] Collating/byte data...\n");
00583 #endif
00584
00585 /* Reorder the time-course of each voxel */
00586 for(kk = 0; kk < mapLength; kk++) {
00587 if(NULL == memcpy((void *)bData[kk]
00588 , (void *)bPtr[map[kk]]
00589 , sizeof(byte)*nvox)) {
00590 THD_delete_3dim_dataset(new_dset, False);
00591 FREEUP(bPtr);
00592 FREEUP(bData);
00593 FREE_WORKSPACE;
00594 return "***********************************\n"
00595 "!! Error performing memory copy !!\n"
00596 "**********************************";
00597 }
00598 perc = (100 * kk) / mapLength; /* display percentage done */
00599 PLUTO_set_meter(plint, perc); /* on the progress meter */
00600 }
00601
00602 PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
00603
00604 if(dupaction) { /* new_dset may change in length again */
00605 /* Verify that averaging is possible:
00606 All duplicates must have the same length.
00607 */
00608 error = 0; /* start with no error */
00609 for(ii = 0; ii < classCount; ii++) {
00610 for(kk = 1; kk < (classCount-1); kk++) {
00611 if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
00612 if(classKRH[ii].length != classKRH[kk].length) {
00613 error = 1; /* lengths differ! */
00614 break;
00615 }
00616 dupsToAverage = True; /* there is at least one duplicate to average */
00617 }
00618 }
00619 if(error) {
00620 PLUTO_popup_message(plint
00621 ,"*********************************************************\n"
00622 "* Duplicate stimulus regions in map file must have the *\n"
00623 "* same length for averaging. Returning collated regions.*\n"
00624 "*********************************************************");
00625 break;
00626 }
00627 }
00628
00629 if(!error) { /* average the crap */
00630 #ifdef DEBUG_PLUGIN_REORDER
00631 printf("[Reorder] All duplicate epochs have the same length...\n");
00632 #endif
00633
00634 if(dupsToAverage) {
00635 #ifdef DEBUG_PLUGIN_REORDER
00636 printf("[Reorder] Averaging duplicate epochs...\n");
00637 printf("[Reorder] Initial time-length is %d...\n", mapLength);
00638 #endif
00639
00640 newLength = mapLength; /* length may change after averaging, track new length */
00641
00642 /* Perform averaging */
00643 /* For each class, do: */
00644 for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
00645 #ifdef DEBUG_PLUGIN_REORDER
00646 printf("[Reorder] Processing instances of class %c...\n"
00647 , classKRH[ii].classKRH);
00648 #endif
00649
00650 /* Get the number of duplicates of the current class */
00651 for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
00652 if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
00653 break;
00654 }
00655 }
00656
00657 #ifdef DEBUG_PLUGIN_REORDER
00658 printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
00659 #endif
00660
00661 if(instances > 1) {
00662 #ifdef DEBUG_PLUGIN_REORDER
00663 printf(" Averaging %d contiguous epochs...\n", instances);
00664 #endif
00665
00666 divisor = (double)instances;
00667
00668 /* For each time point in the duplicated classes: */
00669 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
00670 /* For each voxel: */
00671 for(kk = 0; kk < nvox; kk++) {
00672 /* Average the duplicate time-points for the current voxel */
00673 for(mm = 0, sum = 0.0; mm < instances; mm++) {
00674 sum = sum + (double)bData[jj+(mm*classKRH[ii].length)][kk];
00675 }
00676 bData[timePoint][kk] = (byte)((sum / divisor) + 0.5);
00677 }
00678 }
00679
00680 /* Update final length of time-course */
00681 newLength = newLength - (classKRH[ii].length * (instances - 1));
00682
00683 /* Jump over the indices of all instances of the current class */
00684 startIndex = startIndex + (instances *classKRH[ii].length);
00685
00686 /* Go to the next class to process */
00687 ii += instances;
00688 }
00689 else { /* just copy the data points */
00690 #ifdef DEBUG_PLUGIN_REORDER
00691 printf(" Copying single epoch...\n");
00692 #endif
00693
00694 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
00695 for(kk = 0; kk < nvox; kk++) {
00696 bData[timePoint][kk] = bData[jj][kk];
00697 }
00698 }
00699
00700 /* Jump to the index of the next class in the time-course */
00701 startIndex = startIndex + classKRH[ii].length;
00702
00703 /* eat up another class */
00704 ++ii;
00705 }
00706
00707 /* Update progress meter with the percentage done */
00708 perc = (100 * timePoint) / mapLength;
00709 PLUTO_set_meter(plint, perc);
00710 }
00711
00712 /* Set progress meter to 100% */
00713 PLUTO_set_meter(plint, 100);
00714
00715 #ifdef DEBUG_PLUGIN_REORDER
00716 printf("[Reorder] Final time-length is %d...\n", newLength);
00717 #endif
00718 ii = EDIT_dset_items(
00719 new_dset,
00720 ADN_nvals, newLength, /* # time points */
00721 ADN_ntt, newLength, /* # time points */
00722 ADN_none);
00723 if(0 != ii) {
00724 THD_delete_3dim_dataset(new_dset, False);
00725 FREEUP(bPtr);
00726 FREEUP(bData);
00727 FREE_WORKSPACE;
00728 return "******************************************\n"
00729 "!! Error while creating output dataset !!\n"
00730 "*****************************************";
00731 }
00732
00733 /* Free unused subBRIKS */
00734 if(newLength < mapLength) {
00735 #ifdef DEBUG_PLUGIN_REORDER
00736 printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
00737 #endif
00738 for(ii = newLength; ii < mapLength; ii++) {
00739 FREEUP(bData[ii]);
00740 }
00741 bData = (byte **)realloc((void *)bData, sizeof(byte *)*newLength);
00742 if(NULL == bData) {
00743 THD_delete_3dim_dataset(new_dset, False);
00744 FREEUP(bPtr);
00745 FREE_WORKSPACE;
00746 return "************************\n"
00747 "!! Reallocation error !!\n"
00748 "************************";
00749 }
00750 }
00751 mapLength = newLength;
00752 }
00753 }
00754 }
00755
00756 /* Write results into new dataset */
00757 #ifdef DEBUG_PLUGIN_REORDER
00758 printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
00759 #endif
00760 for(kk = 0; kk < mapLength; kk++) {
00761 EDIT_substitute_brick(new_dset, kk, MRI_byte, bData[kk]);
00762 }
00763 break;
00764
00765 case MRI_short:
00766 #ifdef DEBUG_PLUGIN_REORDER
00767 printf("\n[Reorder] Collating/short data...\n");
00768 #endif
00769
00770 /* Reorder the time-course of each voxel */
00771 for(kk = 0; kk < mapLength; kk++) {
00772 if(NULL == memcpy((void *)sData[kk]
00773 , (void *)sPtr[map[kk]]
00774 , sizeof(short)*nvox)) {
00775 THD_delete_3dim_dataset(new_dset, False);
00776 FREEUP(sPtr);
00777 FREEUP(sData);
00778 FREE_WORKSPACE;
00779 return "***********************************\n"
00780 "!! Error performing memory copy !!\n"
00781 "**********************************";
00782 }
00783 perc = (100 * kk) / mapLength; /* display percentage done */
00784 PLUTO_set_meter(plint, perc); /* on the progress meter */
00785 }
00786
00787 PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
00788
00789 if(dupaction) { /* new_dset may change in length again */
00790 /* Verify that averaging is possible:
00791 All duplicates must have the same length.
00792 */
00793 error = 0; /* start with no error */
00794 for(ii = 0; ii < classCount; ii++) {
00795 for(kk = 1; kk < (classCount-1); kk++) {
00796 if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
00797 if(classKRH[ii].length != classKRH[kk].length) {
00798 error = 1; /* lengths differ! */
00799 break;
00800 }
00801 dupsToAverage = True; /* there is at least one duplicate to average */
00802 }
00803 }
00804 if(error) {
00805 PLUTO_popup_message(plint
00806 ,"*********************************************************\n"
00807 "* Duplicate stimulus regions in map file must have the *\n"
00808 "* same length for averaging. Returning collated regions.*\n"
00809 "*********************************************************");
00810 break;
00811 }
00812 }
00813
00814 if(!error) { /* average the crap */
00815 #ifdef DEBUG_PLUGIN_REORDER
00816 printf("[Reorder] All duplicate epochs have the same length...\n");
00817 #endif
00818
00819 if(dupsToAverage) {
00820 #ifdef DEBUG_PLUGIN_REORDER
00821 printf("[Reorder] Averaging duplicate epochs...\n");
00822 printf("[Reorder] Initial time-length is %d...\n", mapLength);
00823 #endif
00824
00825 newLength = mapLength; /* length may change after averaging, track new length */
00826
00827 /* Perform averaging */
00828 /* For each class, do: */
00829 for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
00830 #ifdef DEBUG_PLUGIN_REORDER
00831 printf("[Reorder] Processing instances of class %c...\n"
00832 , classKRH[ii].classKRH);
00833 #endif
00834
00835 /* Get the number of duplicates of the current class */
00836 for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
00837 if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
00838 break;
00839 }
00840 }
00841
00842 #ifdef DEBUG_PLUGIN_REORDER
00843 printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
00844 #endif
00845
00846 if(instances > 1) {
00847 #ifdef DEBUG_PLUGIN_REORDER
00848 printf(" Averaging %d contiguous epochs...\n", instances);
00849 #endif
00850
00851 divisor = (double)instances;
00852
00853 /* For each time point in the duplicated classes: */
00854 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
00855 /* For each voxel: */
00856 for(kk = 0; kk < nvox; kk++) {
00857 /* Average the duplicate time-points for the current voxel */
00858 for(mm = 0, sum = 0.0; mm < instances; mm++) {
00859 sum = sum + (double)sData[jj+(mm*classKRH[ii].length)][kk];
00860 }
00861 sData[timePoint][kk] = (short)((sum / divisor) + 0.5);
00862 }
00863 }
00864
00865 /* Update final length of time-course */
00866 newLength = newLength - (classKRH[ii].length * (instances - 1));
00867
00868 /* Jump over the indices of all instances of the current class */
00869 startIndex = startIndex + (instances *classKRH[ii].length);
00870
00871 /* Go to the next class to process */
00872 ii += instances;
00873 }
00874 else { /* just copy the data points */
00875 #ifdef DEBUG_PLUGIN_REORDER
00876 printf(" Copying single epoch...\n");
00877 #endif
00878
00879 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
00880 for(kk = 0; kk < nvox; kk++) {
00881 sData[timePoint][kk] = sData[jj][kk];
00882 }
00883 }
00884
00885 /* Jump to the index of the next class in the time-course */
00886 startIndex = startIndex + classKRH[ii].length;
00887
00888 /* eat up another class */
00889 ++ii;
00890 }
00891
00892 /* Update progress meter with the percentage done */
00893 perc = (100 * timePoint) / mapLength;
00894 PLUTO_set_meter(plint, perc);
00895 }
00896
00897 /* Set progress meter to 100% */
00898 PLUTO_set_meter(plint, 100);
00899
00900 #ifdef DEBUG_PLUGIN_REORDER
00901 printf("[Reorder] Final time-length is %d...\n", newLength);
00902 #endif
00903 ii = EDIT_dset_items(
00904 new_dset,
00905 ADN_nvals, newLength, /* # time points */
00906 ADN_ntt, newLength, /* # time points */
00907 ADN_none);
00908 if(0 != ii) {
00909 THD_delete_3dim_dataset(new_dset, False);
00910 FREEUP(sPtr);
00911 FREEUP(sData);
00912 FREE_WORKSPACE;
00913 return "******************************************\n"
00914 "!! Error while creating output dataset !!\n"
00915 "*****************************************";
00916 }
00917
00918 /* Free unused subBRIKS */
00919 if(newLength < mapLength) {
00920 #ifdef DEBUG_PLUGIN_REORDER
00921 printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
00922 #endif
00923 for(ii = newLength; ii < mapLength; ii++) {
00924 FREEUP(sData[ii]);
00925 }
00926 sData = (short **)realloc((void *)sData, sizeof(short *)*newLength);
00927 if(NULL == sData) {
00928 THD_delete_3dim_dataset(new_dset, False);
00929 FREEUP(sPtr);
00930 FREE_WORKSPACE;
00931 return "*************************\n"
00932 "!! Reallocation error !!\n"
00933 "************************";
00934 }
00935 }
00936 mapLength = newLength;
00937 }
00938 }
00939 }
00940
00941 /* Write results into new dataset */
00942 #ifdef DEBUG_PLUGIN_REORDER
00943 printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
00944 #endif
00945 for(kk = 0; kk < mapLength; kk++) {
00946 EDIT_substitute_brick(new_dset, kk, MRI_short, sData[kk]);
00947 }
00948 break;
00949
00950 case MRI_float:
00951 #ifdef DEBUG_PLUGIN_REORDER
00952 printf("\n[Reorder] Collating/float data...\n");
00953 #endif
00954
00955 /* Reorder the time-course of each voxel */
00956 for(kk = 0; kk < mapLength; kk++) {
00957 if(NULL == memcpy((void *)fData[kk]
00958 , (void *)fPtr[map[kk]]
00959 , sizeof(float)*nvox)) {
00960 THD_delete_3dim_dataset(new_dset, False);
00961 FREEUP(fPtr);
00962 FREEUP(fData);
00963 FREE_WORKSPACE;
00964 return "***********************************\n"
00965 "!! Error performing memory copy !!\n"
00966 "**********************************";
00967 }
00968 perc = (100 * kk) / mapLength; /* display percentage done */
00969 PLUTO_set_meter(plint, perc); /* on the progress meter */
00970 }
00971
00972 PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
00973
00974 if(dupaction) { /* new_dset may change in length again */
00975 /* Verify that averaging is possible:
00976 All duplicates must have the same length.
00977 */
00978 error = 0; /* start with no error */
00979 for(ii = 0; ii < classCount; ii++) {
00980 for(kk = 1; kk < (classCount-1); kk++) {
00981 if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
00982 if(classKRH[ii].length != classKRH[kk].length) {
00983 error = 1; /* lengths differ! */
00984 break;
00985 }
00986 dupsToAverage = True; /* there is at least one duplicate to average */
00987 }
00988 }
00989 if(error) {
00990 PLUTO_popup_message(plint
00991 ,"*********************************************************\n"
00992 "* Duplicate stimulus regions in map file must have the *\n"
00993 "* same length for averaging. Returning collated regions.*\n"
00994 "*********************************************************");
00995 break;
00996 }
00997 }
00998
00999 if(!error) { /* average the crap */
01000 #ifdef DEBUG_PLUGIN_REORDER
01001 printf("[Reorder] All duplicate epochs have the same length...\n");
01002 #endif
01003
01004 if(dupsToAverage) {
01005 #ifdef DEBUG_PLUGIN_REORDER
01006 printf("[Reorder] Averaging duplicate epochs...\n");
01007 printf("[Reorder] Initial time-length is %d...\n", mapLength);
01008 #endif
01009
01010 newLength = mapLength; /* length may change after averaging, track new length */
01011
01012 /* Perform averaging */
01013 /* For each class, do: */
01014 for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
01015 #ifdef DEBUG_PLUGIN_REORDER
01016 printf("[Reorder] Processing instances of class %c...\n"
01017 , classKRH[ii].classKRH);
01018 #endif
01019
01020 /* Get the number of duplicates of the current class */
01021 for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
01022 if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
01023 break;
01024 }
01025 }
01026
01027 #ifdef DEBUG_PLUGIN_REORDER
01028 printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
01029 #endif
01030
01031 if(instances > 1) {
01032 #ifdef DEBUG_PLUGIN_REORDER
01033 printf(" Averaging %d contiguous epochs...\n", instances);
01034 #endif
01035
01036 divisor = (double)instances;
01037
01038 /* For each time point in the duplicated classes: */
01039 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
01040 /* For each voxel: */
01041 for(kk = 0; kk < nvox; kk++) {
01042 /* Average the duplicate time-points for the current voxel */
01043 for(mm = 0, sum = 0.0; mm < instances; mm++) {
01044 sum = sum + (double)fData[jj+(mm*classKRH[ii].length)][kk];
01045 }
01046 fData[timePoint][kk] = (float)((sum / divisor) + 0.5);
01047 }
01048 }
01049
01050 /* Update final length of time-course */
01051 newLength = newLength - (classKRH[ii].length * (instances - 1));
01052
01053 /* Jump over the indices of all instances of the current class */
01054 startIndex = startIndex + (instances *classKRH[ii].length);
01055
01056 /* Go to the next class to process */
01057 ii += instances;
01058 }
01059 else { /* just copy the data points */
01060 #ifdef DEBUG_PLUGIN_REORDER
01061 printf(" Copying single epoch...\n");
01062 #endif
01063
01064 for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
01065 for(kk = 0; kk < nvox; kk++) {
01066 fData[timePoint][kk] = fData[jj][kk];
01067 }
01068 }
01069
01070 /* Jump to the index of the next class in the time-course */
01071 startIndex = startIndex + classKRH[ii].length;
01072
01073 /* eat up another class */
01074 ++ii;
01075 }
01076
01077 /* Update progress meter with the percentage done */
01078 perc = (100 * timePoint) / mapLength;
01079 PLUTO_set_meter(plint, perc);
01080 }
01081
01082 /* Set progress meter to 100% */
01083 PLUTO_set_meter(plint, 100);
01084
01085 #ifdef DEBUG_PLUGIN_REORDER
01086 printf("[Reorder] Final time-length is %d...\n", newLength);
01087 #endif
01088 ii = EDIT_dset_items(
01089 new_dset,
01090 ADN_nvals, newLength, /* # time points */
01091 ADN_ntt, newLength, /* # time points */
01092 ADN_none);
01093 if(0 != ii) {
01094 THD_delete_3dim_dataset(new_dset, False);
01095 FREEUP(fPtr);
01096 FREEUP(fData);
01097 FREE_WORKSPACE;
01098 return "******************************************\n"
01099 "!! Error while creating output dataset !!\n"
01100 "*****************************************";
01101 }
01102
01103 /* Free unused subBRIKS */
01104 if(newLength < mapLength) {
01105 #ifdef DEBUG_PLUGIN_REORDER
01106 printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
01107 #endif
01108 for(ii = newLength; ii < mapLength; ii++) {
01109 FREEUP(fData[ii]);
01110 }
01111 fData = (float **)realloc((void *)fData, sizeof(float *)*newLength);
01112 if(NULL == fData) {
01113 THD_delete_3dim_dataset(new_dset, False);
01114 FREEUP(fPtr);
01115 FREE_WORKSPACE;
01116 return "*************************\n"
01117 "!! Reallocation error !!\n"
01118 "************************";
01119 }
01120 }
01121 mapLength = newLength;
01122 }
01123 }
01124 }
01125
01126 /* Write results into new dataset */
01127 #ifdef DEBUG_PLUGIN_REORDER
01128 printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
01129 #endif
01130 for(kk = 0; kk < mapLength; kk++) {
01131 EDIT_substitute_brick(new_dset, kk, MRI_float, fData[kk]);
01132 }
01133 break;
01134
01135 } /* end of switch over input type */
01136
01137 DSET_unload(old_dset);
01138
01139 /*************** Cleanup and go home *****************/
01140
01141 PLUTO_add_dset(plint, new_dset, DSET_ACTION_MAKE_CURRENT);
01142
01143 #ifdef DEBUG_PLUGIN_REORDER
01144 printf("[Reorder] Cleaning up workspace...\n");
01145 #endif
01146 FREE_WORKSPACE;
01147
01148 return NULL; /* NULL string returned means all was OK */
01149 }
|
Variable Documentation
|
|
Initial value:
{ "Collate", "Average" }Definition at line 101 of file plug_reorder.c. Referenced by PLUGIN_init(), and REORDER_main(). |
|
|
Definition at line 42 of file plug_reorder.c. Referenced by PLUGIN_init(). |
