00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #ifndef __MAIN_PLUGIN_REORDER_C__
00008 #define __MAIN_PLUGIN_REORDER_C__
00009 
00010 /***************************************************
00011 *  File:   plug_reorder.c
00012 *
00013 *  Author: Jay Brian Kummer, Medical College of WI
00014 *          [kummer@its.mcw.edu]
00015 *
00016 *  Date:   May 8, 1997
00017 *
00018 *  Purpose:
00019 *
00020 *    Plugin to reorder the epochs (i.e., stimulus
00021 *    presentation regions) in a 3D+Time data set.
00022 ****************************************************/
00023 
00024 /***********************************************************************/
00025 /* 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 */
00026 /***********************************************************************/
00027 
00028 #define MAIN_PLUGIN_REORDER
00029 /* Uncomment to make this plugin verbose:
00030 #define DEBUG_PLUGIN_REORDER
00031 */
00032 
00033 /*********************************/
00034 /* D A T A   S T R U C T U R E S */
00035 /*********************************/
00036 
00037 typedef struct { /* returned from 'plug_reorder_parseMap' */
00038         char classKRH;  /* the class (i.e., [a-zA-Z]) of a time-course segment */
00039         int length;  /* the length of the class in time points */
00040         } ClassInfo;
00041 
00042 static char helpstring[] = /* Help string for plugin interface */
00043                 "Purpose: Reorders voxel time-courses within a 3D+time BRIK.\n"
00044                 "Input  : Name of the 3D+time dataset to reorder.\n"
00045                 "\n"
00046                 "Output : Prefix of filename for the new dataset.\n"
00047                 "\n"
00048                 "Map file : Text file with instructions a reordering the input dataset.\n"
00049                 "       Sample map file:\n"
00050                 "               # Experiment: Bogus\n"
00051                 "               # Date: April 1, 1997\n"
00052                 "               -\n"
00053                 "               -\n"
00054                 "               D1\n"
00055                 "               D2\n"
00056                 "               D3\n"
00057                 "               B1\n"
00058                 "               B2\n"
00059                 "               B3\n"
00060                 "               B4\n"
00061                 "               -\n"
00062                 "               a1\n"
00063                 "               a2\n"
00064                 "               a3\n"
00065                 "               A1\n"
00066                 "               A2\n"
00067                 "               A3\n"
00068                 "               A1\n"
00069                 "               A2\n"
00070                 "               A3\n"
00071                 "               -\n"
00072                 "\n"
00073                 "       # : Comment line (ignored)\n"
00074                 "       - : Point to be omitted from the output time-series. In the sample map,\n"
00075                 "               time points 1, 2, 10, and 20 will be omitted.\n"
00076                 "       Each non-comment line in a map file corresponds to a point in the time-series of\n"
00077                 "       an input 3D+time dataset. Codes for each time point define the plug-in action for\n"
00078                 "       those points. For example, the code in the fifth non-comment line of the sample\n"
00079                 "       above, 'D3', means that the fifth point in the time-series is the third point in\n"
00080                 "       the presentation of a stimulus named 'D'. \n"
00081                 "               Stimulus codes are single alphabetic characters and are case-sensitive, with\n"
00082                 "       uppercase preceding lowercase in the reordering (i.e., A-Z then a-z). Stimulus\n"
00083                 "       regions are numbered starting from 1 (i.e., [ X0 X1 X2 ] and [ X1 X3 X5 ] will\n"
00084                 "       cause the plug-in to abort), and, as shown in the sample with the code 'A', this\n"
00085                 "       delimits duplicate instances of the same stimulus code. For the sample map, images\n"
00086                 "       would be reordered so that the output time-series is:\n"
00087                 "                                       [ A1 A2 A3 A1 A2 A3 B1 B2 B3 B4 C1 C2 C3 a1 a2 a3 ]\n"
00088                 "\n"
00089                 "Duplicates : If there are duplicates of a stimulus presentation within a map (e.g.,\n"
00090                 "[ A1 A2 A3 ] in the sample), they can be:\n"
00091                 "       1) Collated: Duplicate instances of a stimulus will be placed in the order in\n"
00092                 "               which they occur in the map. In the case of the sample, this results in the\n"
00093                 "               ordering [ A1 A2 A3 A1 A2 A3 ... ].\n"
00094                 "       2) Averaged: Duplicate instances of a stimulus will be averaged. In the sample,\n"
00095                 "               this would result in the ordering [ A1' A2' A3' ... ], where Ai' is the mean\n"
00096                 "               of the i-th point of all of the stimulus A's present in the map file.\n"
00097                 "       **Duplicate stimulus regions must have the same length! For example,\n"
00098                 "       [ D1 D2 X1 D1 D2 D3 ] will cause the plug-in to abort.\n"
00099                 ;
00100 
00101 static char *dupaction_strings[] = /* strings for duplicate action in plugin interface */ 
00102                 { "Collate", "Average" };
00103 
00104 #define NUM_DUPACTION_STRINGS (sizeof(dupaction_strings)/sizeof(char *))
00105 
00106 /***************************/
00107 /* H E A D E R   F I L E S */
00108 /***************************/
00109 
00110 #include "afni.h"
00111 
00112 #include <unistd.h>
00113 #include "plug_reorder_parseMap.c"
00114 
00115 /***********************************************************/
00116 /* 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 */
00117 /***********************************************************/
00118 
00119 char *REORDER_main(PLUGIN_interface *);/* the entry point */
00120 
00121 /*************************************/
00122 /* S U P P O R T   F U N C T I O N S */
00123 /*************************************/
00124 
00125 int compare_class(const void *i, const void *j)
00126 /* qsort compare routine for
00127    sorting 'ClassInfo' array
00128 */
00129 {
00130 ClassInfo *iTmp = (ClassInfo *)i;
00131 ClassInfo *jTmp = (ClassInfo *)j;
00132 
00133 if(iTmp->classKRH == jTmp->classKRH) {
00134         return(0);
00135         }
00136 else if(iTmp->classKRH < jTmp->classKRH) {
00137         return(-1);
00138         }
00139 else {
00140         return(1);
00141         }
00142 }
00143 
00144 /***********************************************************************
00145 *  Set up the interface to the user:
00146 *  1) Create a new interface using "PLUTO_new_interface";
00147 *  
00148 *  2) For each line of inputs, create the line with "PLUTO_add_option"
00149 *       (this line of inputs can be optional or mandatory);
00150 *  
00151 *  3) For each item on the line, create the item with
00152 *       "PLUTO_add_dataset" for a dataset chooser,
00153 *       "PLUTO_add_string"for a string chooser,
00154 *       "PLUTO_add_number"for a number chooser.
00155 ************************************************************************/
00156 
00157 
00158 DEFINE_PLUGIN_PROTOTYPE
00159 
00160 PLUGIN_interface *PLUGIN_init(int ncall)
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 }
00233 
00234 /*************************************/
00235 /* M A C R O   D E F I N I T I O N S */
00236 /*************************************/
00237 
00238 /* Free vectors  */
00239 #undef  FREEUP
00240 #define FREEUP(thing) if((thing) != NULL) {free((thing)); (thing) = NULL;}
00241 
00242 #define FREE_WORKSPACE do{ FREEUP(bPtr); FREEUP(sPtr); FREEUP(fPtr); \
00243                                                    FREEUP(map); FREEUP(classKRH); } while(0)
00244 
00245 /*****************************************************************
00246 *  Main routine for this plugin (will be called from AFNI).
00247 *  If the return string is not NULL, some error transpired, and
00248 *  AFNI will popup the return string in a message box.
00249 ******************************************************************/
00250 
00251 char *REORDER_main(PLUGIN_interface *plint)
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 }
01150 
01151 #endif
01152