/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#ifndef __MAIN_PLUGIN_REORDER_C__
#define __MAIN_PLUGIN_REORDER_C__
/***************************************************
* File: plug_reorder.c
*
* Author: Jay Brian Kummer, Medical College of WI
* [kummer@its.mcw.edu]
*
* Date: May 8, 1997
*
* Purpose:
*
* Plugin to reorder the epochs (i.e., stimulus
* presentation regions) in a 3D+Time data set.
****************************************************/
/***********************************************************************/
/* 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 */
/***********************************************************************/
#define MAIN_PLUGIN_REORDER
/* Uncomment to make this plugin verbose:
#define DEBUG_PLUGIN_REORDER
*/
/*********************************/
/* D A T A S T R U C T U R E S */
/*********************************/
typedef struct { /* returned from 'plug_reorder_parseMap' */
char classKRH; /* the class (i.e., [a-zA-Z]) of a time-course segment */
int length; /* the length of the class in time points */
} ClassInfo;
static char helpstring[] = /* Help string for plugin interface */
"Purpose: Reorders voxel time-courses within a 3D+time BRIK.\n"
"Input : Name of the 3D+time dataset to reorder.\n"
"\n"
"Output : Prefix of filename for the new dataset.\n"
"\n"
"Map file : Text file with instructions a reordering the input dataset.\n"
" Sample map file:\n"
" # Experiment: Bogus\n"
" # Date: April 1, 1997\n"
" -\n"
" -\n"
" D1\n"
" D2\n"
" D3\n"
" B1\n"
" B2\n"
" B3\n"
" B4\n"
" -\n"
" a1\n"
" a2\n"
" a3\n"
" A1\n"
" A2\n"
" A3\n"
" A1\n"
" A2\n"
" A3\n"
" -\n"
"\n"
" # : Comment line (ignored)\n"
" - : Point to be omitted from the output time-series. In the sample map,\n"
" time points 1, 2, 10, and 20 will be omitted.\n"
" Each non-comment line in a map file corresponds to a point in the time-series of\n"
" an input 3D+time dataset. Codes for each time point define the plug-in action for\n"
" those points. For example, the code in the fifth non-comment line of the sample\n"
" above, 'D3', means that the fifth point in the time-series is the third point in\n"
" the presentation of a stimulus named 'D'. \n"
" Stimulus codes are single alphabetic characters and are case-sensitive, with\n"
" uppercase preceding lowercase in the reordering (i.e., A-Z then a-z). Stimulus\n"
" regions are numbered starting from 1 (i.e., [ X0 X1 X2 ] and [ X1 X3 X5 ] will\n"
" cause the plug-in to abort), and, as shown in the sample with the code 'A', this\n"
" delimits duplicate instances of the same stimulus code. For the sample map, images\n"
" would be reordered so that the output time-series is:\n"
" [ A1 A2 A3 A1 A2 A3 B1 B2 B3 B4 C1 C2 C3 a1 a2 a3 ]\n"
"\n"
"Duplicates : If there are duplicates of a stimulus presentation within a map (e.g.,\n"
"[ A1 A2 A3 ] in the sample), they can be:\n"
" 1) Collated: Duplicate instances of a stimulus will be placed in the order in\n"
" which they occur in the map. In the case of the sample, this results in the\n"
" ordering [ A1 A2 A3 A1 A2 A3 ... ].\n"
" 2) Averaged: Duplicate instances of a stimulus will be averaged. In the sample,\n"
" this would result in the ordering [ A1' A2' A3' ... ], where Ai' is the mean\n"
" of the i-th point of all of the stimulus A's present in the map file.\n"
" **Duplicate stimulus regions must have the same length! For example,\n"
" [ D1 D2 X1 D1 D2 D3 ] will cause the plug-in to abort.\n"
;
static char *dupaction_strings[] = /* strings for duplicate action in plugin interface */
{ "Collate", "Average" };
#define NUM_DUPACTION_STRINGS (sizeof(dupaction_strings)/sizeof(char *))
/***************************/
/* H E A D E R F I L E S */
/***************************/
#include "afni.h"
#include <unistd.h>
#include "plug_reorder_parseMap.c"
/***********************************************************/
/* 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 */
/***********************************************************/
char *REORDER_main(PLUGIN_interface *);/* the entry point */
/*************************************/
/* S U P P O R T F U N C T I O N S */
/*************************************/
int compare_class(const void *i, const void *j)
/* qsort compare routine for
sorting 'ClassInfo' array
*/
{
ClassInfo *iTmp = (ClassInfo *)i;
ClassInfo *jTmp = (ClassInfo *)j;
if(iTmp->classKRH == jTmp->classKRH) {
return(0);
}
else if(iTmp->classKRH < jTmp->classKRH) {
return(-1);
}
else {
return(1);
}
}
/***********************************************************************
* Set up the interface to the user:
* 1) Create a new interface using "PLUTO_new_interface";
*
* 2) For each line of inputs, create the line with "PLUTO_add_option"
* (this line of inputs can be optional or mandatory);
*
* 3) For each item on the line, create the item with
* "PLUTO_add_dataset" for a dataset chooser,
* "PLUTO_add_string"for a string chooser,
* "PLUTO_add_number"for a number chooser.
************************************************************************/
DEFINE_PLUGIN_PROTOTYPE
PLUGIN_interface *PLUGIN_init(int ncall)
{
PLUGIN_interface *plint; /* will be the output of this routine */
if(ncall >= 1) { /* only one interface */
return NULL;
}
/*---------------- set titles and call point ----------------*/
plint = PLUTO_new_interface("Reorder",
"Reorders voxel time-courses within a 3D+Time BRIK",
helpstring,
PLUGIN_CALL_VIA_MENU,
REORDER_main);
PLUTO_set_sequence( plint , "z:Kummer" ) ;
/*--------- 1st line: Input dataset ---------*/
PLUTO_add_option(plint,
"Input",/* label at left of input line */
"Input",/* tag to return to plugin */
TRUE /* is this mandatory? */
);
PLUTO_add_dataset(plint,
"---->>", /* label next to button */
ANAT_ALL_MASK, /* take any anat datasets */
0, /* don't allow fim funcs*/
DIMEN_4D_MASK | /* need 3D+time datasets*/
BRICK_ALLTYPE_MASK /* handle any type */
);
/*---------- 2nd line: Output dataset prefix ----------*/
PLUTO_add_option(plint,
"Output", /* label at left of input line */
"Output", /* tag to return to plugin */
TRUE /* is this mandatory? */
);
PLUTO_add_string(plint,
"Prefix", /* label next to textfield */
0,NULL, /* no fixed strings to choose among */
19 /* 19 spaces for typing in value */
);
/*---------- 3rd line: Map file name ----------*/
PLUTO_add_option(plint,
"Map file", /* label at left of input line */
"MapFile", /* tag to return to plugin */
TRUE /* is this mandatory? */
);
PLUTO_add_string(plint,
"FileName", /* label next to textfield */
0,NULL, /* no fixed strings to choose among */
19 /* 19 spaces for typing in value */
);
PLUTO_add_string(plint,
"Duplicates", /* label next to textfield */
NUM_DUPACTION_STRINGS, /* number of strings to choose among */
dupaction_strings, /* fixed strings to choose among */
0 /* index of default (collate) */
);
/*--------- done with interface setup ---------*/
return plint;
}
/*************************************/
/* M A C R O D E F I N I T I O N S */
/*************************************/
/* Free vectors */
#undef FREEUP
#define FREEUP(thing) if((thing) != NULL) {free((thing)); (thing) = NULL;}
#define FREE_WORKSPACE do{ FREEUP(bPtr); FREEUP(sPtr); FREEUP(fPtr); \
FREEUP(map); FREEUP(classKRH); } while(0)
/*****************************************************************
* Main routine for this plugin (will be called from AFNI).
* If the return string is not NULL, some error transpired, and
* AFNI will popup the return string in a message box.
******************************************************************/
char *REORDER_main(PLUGIN_interface *plint)
{
MCW_idcode *idc; /* input dataset idcode */
THD_3dim_dataset *old_dset; /* input dataset */
THD_3dim_dataset *new_dset; /* output dataset */
ClassInfo *classKRH; /* array of ClassInfo structures, returned by reference from 'parseMap' */
char *new_prefix; /* pointer to the output dataset prefix */
char *mapFile; /* pointer to the name of map file to open */
char *dupstr; /* pointer to duplicate action label string */
char currentClass; /* character sentinel */
register int ii; /* register counter variables */
register int jj;
register int kk;
register int mm;
int perc; /* percentage value for progress meter */
int ninp; /* number of time points */
int nvox; /* number of voxels */
int old_datum; /* data type of input dataset */
int error; /* general status variable */
int startIndex; /* index sentinel */
int currentLength; /* length of the current epoch */
int newLength; /* length of the time-course after averaging duplicates */
int instances; /* number of instances of the current epoch class */
int timePoint; /* current time point being processed */
int *map; /* array of indices for remapping time points, returned from 'parseMap' */
int mapLength; /* length of 'map' array, returned by reference from 'parseMap' */
int dupaction; /* 0: collate duplicate epochs, !0: average duplicate epochs */
int classCount; /* length of 'ClassInfo' array, returned by reference from 'parseMap' */
double sum; /* sum variable for averaging */
double divisor; /* divisor for calculating the mean, this will equal 'instances' */
byte **bPtr = NULL; /* for byte data, will point to the 3D+time data to be reordered */
short **sPtr = NULL; /* for short data, will point to the 3D+time data to be reordered */
float **fPtr = NULL; /* for float data, will point to the 3D+time data to be reordered */
byte **bData = NULL; /* for byte data, will point to the reordered 3D+time data */
short **sData = NULL; /* for short data, will point to the reordered 3D+time data */
float **fData = NULL; /* for float data, will point to the reordered 3D+time data */
Bool dupsToAverage = False; /* flag to indicate whether there are segments to average */
/**********************************************************************/
/****** Check inputs from AFNI to see if they are reasonable-ish ******/
/*--------- go to first input line ---------*/
PLUTO_next_option(plint);
idc = PLUTO_get_idcode(plint); /* get dataset item */
old_dset = PLUTO_find_dset(idc); /* get Ptr to dataset */
if(NULL == old_dset) {
return "*************************\n"
"Cannot find Input Dataset\n"
"*************************";
}
/*--------- go to second input line ---------*/
PLUTO_next_option(plint);
new_prefix = PLUTO_get_string(plint); /* get string item (the output prefix) */
if(! PLUTO_prefix_ok(new_prefix)) { /* check if it is OK */
return "************************\n"
"Output Prefix is illegal\n"
"************************";
}
/*--------- go to next input lines ---------*/
PLUTO_next_option(plint); /* skip to next line */
mapFile = PLUTO_get_string(plint); /* get the map filename */
if(NULL == mapFile || access(mapFile, R_OK) < 0) { /* check if it is OK and readable */
return "****************************\n"
"Epoch map file is unreadable\n"
"****************************";
}
dupstr = PLUTO_get_string(plint); /* get the string item */
dupaction = PLUTO_string_index(dupstr, /* find it in the list it came from */
NUM_DUPACTION_STRINGS,
dupaction_strings);
/********************************************************/
/*********** At this point, the inputs are OK ***********/
PLUTO_popup_meter(plint); /* popup a progress meter */
/*--------- set up pointers to each sub-brick in the input dataset ---------*/
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Loading dataset into memory...\n");
#endif
DSET_load(old_dset); /* must be in memory before we get pointers to it */
nvox = /* number of voxels per sub-brick */
old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] %d voxels per subBRIK [%dx%dx%d].\n"
, nvox, old_dset->daxes->nxx, old_dset->daxes->nyy, old_dset->daxes->nzz);
#endif
ninp = DSET_NUM_TIMES(old_dset); /* get number of time points in old dataset */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] %d time points [subBRIK].\n", ninp);
#endif
/********************************************************/
/*************** Parse the epoch map file ***************/
mapLength = ninp; /* initialize with the expected length for mapFile entries */
if(NULL == (map = REORDER_parseMap(mapFile
, &mapLength
, &classKRH
, &classCount))) {
return "*******************************\n"
"Critical error parsing map file\n"
"*******************************";
}
#ifdef DEBUG_PLUGIN_REORDER0
printf("\n[Reorder] Indices for epoch remapping:\n");
for(ii = 0; ii < mapLength; ii++) {
printf("%d\n", map[ii]);
}
printf("\n[Reorder] Meta-sequence of epoch classes is:\n");
for(ii = 0; ii < classCount; ii++) {
printf("[%d] Class %c [Width %d TRs]\n", ii, classKRH[ii].classKRH, classKRH[ii].length);
}
#endif
/* Sort 'class' to reflect the final order
to aid in any epoch averaging required */
qsort((void *)classKRH, classCount, sizeof(ClassInfo), compare_class);
#ifdef DEBUG_PLUGIN_REORDER
printf("\n[Reorder] Sorted meta-sequence of epoch classes is:\n");
for(ii = 0; ii < classCount; ii++) {
printf("[%d] Class %c [Width %d TRs]\n", ii, classKRH[ii].classKRH, classKRH[ii].length);
}
#endif
old_datum /* get old dataset datum type */
= DSET_BRICK_TYPE(old_dset, 0);
switch(old_datum) { /* pointer type depends on input datum type */
default:
return "*******************************\n"
"Illegal datum in Input Dataset\n"
"******************************";
/**create array of pointers into old dataset sub-bricks **/
/**Note that we skip the first 'ignore' sub-bricks here **/
/*---------- input is bytes -----------*/
/* voxel #i at time #k is bPtr[k][i] */
/* for i = 0..nvox-1 and k = 0..ninp-1.*/
case MRI_byte:
/* Workspace for 3D+time data */
bPtr = (byte **)calloc(sizeof(byte *), ninp);
if(NULL == bPtr) {
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
for(kk = 0; kk < ninp; kk++) {
bPtr[kk] = (byte *)DSET_ARRAY(old_dset,kk);
}
/* Array of pointers for reordering subBRIKS */
bData = (byte **)calloc(sizeof(byte *), mapLength);
if(NULL == bData) {
FREEUP(bPtr);
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
/* Allocate memory for subBRIKS */
for(kk = 0; kk < mapLength; kk++) {
bData[kk] = (byte *)calloc(sizeof(byte), nvox);
if(NULL == bData[kk]) {
FREEUP(bPtr);
FREEUP(bData);
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
}
break;
/*---------- input is shorts ----------*/
/* voxel #i at time #k is sPtr[k][i] */
/* for i = 0..nvox-1 and k = 0..ninp-1.*/
case MRI_short:
/* Workspace for 3D+time data */
sPtr = (short **)calloc(sizeof(short *), ninp);
if(NULL == sPtr) {
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
for(kk = 0; kk < ninp; kk++) {
sPtr[kk] = (short *)DSET_ARRAY(old_dset,kk);
}
/* Array of pointers for reordering subBRIKS */
sData = (short **)calloc(sizeof(short *), mapLength);
if(NULL == sData) {
FREEUP(sPtr);
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
/* Allocate memory for subBRIKS */
for(kk = 0; kk < mapLength; kk++) {
sData[kk] = (short *)calloc(sizeof(short), nvox);
if(NULL == sData[kk]) {
FREEUP(sPtr);
FREEUP(sData);
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
}
break;
/*---------- input is floats ----------*/
/* voxel #i at time #k is fPtr[k][i] */
/* for i = 0..nvox-1 and k = 0..ninp-1.*/
case MRI_float:
/* Workspace for 3D+time data */
fPtr = (float **)calloc(sizeof(float *), ninp);
if(NULL == fPtr) {
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
for(kk = 0; kk < ninp; kk++) {
fPtr[kk] = (float *)DSET_ARRAY(old_dset,kk);
}
/* Array of pointers for reordering subBRIKS */
fData = (float **)calloc(sizeof(float *), mapLength);
if(NULL == fData) {
FREEUP(fPtr);
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
/* Allocate memory for subBRIKS */
for(kk = 0; kk < mapLength; kk++) {
fData[kk] = (float *)calloc(sizeof(float), nvox);
if(NULL == fData[kk]) {
FREEUP(fPtr);
FREEUP(fData);
return "************************\n"
"!! Allocation Failure !!\n"
"************************";
}
}
break;
} /* end of switch on input type */
/*********************** Make a new dataset ***********************/
new_dset /* start with copy of old one */
= EDIT_empty_copy(old_dset);
{ char * his = PLUTO_commandstring(plint) ;
tross_Copy_History(old_dset,new_dset) ;
tross_Append_History(new_dset,his) ; free(his);
}
/*-- edit some of its internal parameters --*/
ii = EDIT_dset_items(
new_dset,
ADN_prefix , new_prefix, /* filename prefix */
ADN_malloc_type, DATABLOCK_MEM_MALLOC, /* store in memory */
ADN_nvals , mapLength, /* # time points */
ADN_ntt , mapLength, /* # time points */
ADN_none);
if(0 != ii) {
THD_delete_3dim_dataset(new_dset, False);
switch(old_datum) {
case MRI_byte:
FREEUP(bPtr);
FREEUP(bData);
break;
case MRI_short:
FREEUP(sPtr);
FREEUP(sData);
break;
case MRI_float:
FREEUP(fPtr);
FREEUP(fData);
break;
}
FREE_WORKSPACE;
return "***********************************\n"
"Error while creating output dataset\n"
"***********************************";
}
/*****************************************************/
/****** Setup has ended. Now do some real work. ******/
switch(old_datum) {
/*****************************************/
/* voxel #i at time #k is array[kk][ii] */
/* for ii = 0..nvox-1 and kk = 0..ninp-1.*/
case MRI_byte:
#ifdef DEBUG_PLUGIN_REORDER
printf("\n[Reorder] Collating/byte data...\n");
#endif
/* Reorder the time-course of each voxel */
for(kk = 0; kk < mapLength; kk++) {
if(NULL == memcpy((void *)bData[kk]
, (void *)bPtr[map[kk]]
, sizeof(byte)*nvox)) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(bPtr);
FREEUP(bData);
FREE_WORKSPACE;
return "***********************************\n"
"!! Error performing memory copy !!\n"
"**********************************";
}
perc = (100 * kk) / mapLength; /* display percentage done */
PLUTO_set_meter(plint, perc); /* on the progress meter */
}
PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
if(dupaction) { /* new_dset may change in length again */
/* Verify that averaging is possible:
All duplicates must have the same length.
*/
error = 0; /* start with no error */
for(ii = 0; ii < classCount; ii++) {
for(kk = 1; kk < (classCount-1); kk++) {
if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
if(classKRH[ii].length != classKRH[kk].length) {
error = 1; /* lengths differ! */
break;
}
dupsToAverage = True; /* there is at least one duplicate to average */
}
}
if(error) {
PLUTO_popup_message(plint
,"*********************************************************\n"
"* Duplicate stimulus regions in map file must have the *\n"
"* same length for averaging. Returning collated regions.*\n"
"*********************************************************");
break;
}
}
if(!error) { /* average the crap */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] All duplicate epochs have the same length...\n");
#endif
if(dupsToAverage) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Averaging duplicate epochs...\n");
printf("[Reorder] Initial time-length is %d...\n", mapLength);
#endif
newLength = mapLength; /* length may change after averaging, track new length */
/* Perform averaging */
/* For each class, do: */
for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Processing instances of class %c...\n"
, classKRH[ii].classKRH);
#endif
/* Get the number of duplicates of the current class */
for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
break;
}
}
#ifdef DEBUG_PLUGIN_REORDER
printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
#endif
if(instances > 1) {
#ifdef DEBUG_PLUGIN_REORDER
printf(" Averaging %d contiguous epochs...\n", instances);
#endif
divisor = (double)instances;
/* For each time point in the duplicated classes: */
for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
/* For each voxel: */
for(kk = 0; kk < nvox; kk++) {
/* Average the duplicate time-points for the current voxel */
for(mm = 0, sum = 0.0; mm < instances; mm++) {
sum = sum + (double)bData[jj+(mm*classKRH[ii].length)][kk];
}
bData[timePoint][kk] = (byte)((sum / divisor) + 0.5);
}
}
/* Update final length of time-course */
newLength = newLength - (classKRH[ii].length * (instances - 1));
/* Jump over the indices of all instances of the current class */
startIndex = startIndex + (instances *classKRH[ii].length);
/* Go to the next class to process */
ii += instances;
}
else { /* just copy the data points */
#ifdef DEBUG_PLUGIN_REORDER
printf(" Copying single epoch...\n");
#endif
for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
for(kk = 0; kk < nvox; kk++) {
bData[timePoint][kk] = bData[jj][kk];
}
}
/* Jump to the index of the next class in the time-course */
startIndex = startIndex + classKRH[ii].length;
/* eat up another class */
++ii;
}
/* Update progress meter with the percentage done */
perc = (100 * timePoint) / mapLength;
PLUTO_set_meter(plint, perc);
}
/* Set progress meter to 100% */
PLUTO_set_meter(plint, 100);
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Final time-length is %d...\n", newLength);
#endif
ii = EDIT_dset_items(
new_dset,
ADN_nvals, newLength, /* # time points */
ADN_ntt, newLength, /* # time points */
ADN_none);
if(0 != ii) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(bPtr);
FREEUP(bData);
FREE_WORKSPACE;
return "******************************************\n"
"!! Error while creating output dataset !!\n"
"*****************************************";
}
/* Free unused subBRIKS */
if(newLength < mapLength) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
#endif
for(ii = newLength; ii < mapLength; ii++) {
FREEUP(bData[ii]);
}
bData = (byte **)realloc((void *)bData, sizeof(byte *)*newLength);
if(NULL == bData) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(bPtr);
FREE_WORKSPACE;
return "************************\n"
"!! Reallocation error !!\n"
"************************";
}
}
mapLength = newLength;
}
}
}
/* Write results into new dataset */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
#endif
for(kk = 0; kk < mapLength; kk++) {
EDIT_substitute_brick(new_dset, kk, MRI_byte, bData[kk]);
}
break;
case MRI_short:
#ifdef DEBUG_PLUGIN_REORDER
printf("\n[Reorder] Collating/short data...\n");
#endif
/* Reorder the time-course of each voxel */
for(kk = 0; kk < mapLength; kk++) {
if(NULL == memcpy((void *)sData[kk]
, (void *)sPtr[map[kk]]
, sizeof(short)*nvox)) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(sPtr);
FREEUP(sData);
FREE_WORKSPACE;
return "***********************************\n"
"!! Error performing memory copy !!\n"
"**********************************";
}
perc = (100 * kk) / mapLength; /* display percentage done */
PLUTO_set_meter(plint, perc); /* on the progress meter */
}
PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
if(dupaction) { /* new_dset may change in length again */
/* Verify that averaging is possible:
All duplicates must have the same length.
*/
error = 0; /* start with no error */
for(ii = 0; ii < classCount; ii++) {
for(kk = 1; kk < (classCount-1); kk++) {
if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
if(classKRH[ii].length != classKRH[kk].length) {
error = 1; /* lengths differ! */
break;
}
dupsToAverage = True; /* there is at least one duplicate to average */
}
}
if(error) {
PLUTO_popup_message(plint
,"*********************************************************\n"
"* Duplicate stimulus regions in map file must have the *\n"
"* same length for averaging. Returning collated regions.*\n"
"*********************************************************");
break;
}
}
if(!error) { /* average the crap */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] All duplicate epochs have the same length...\n");
#endif
if(dupsToAverage) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Averaging duplicate epochs...\n");
printf("[Reorder] Initial time-length is %d...\n", mapLength);
#endif
newLength = mapLength; /* length may change after averaging, track new length */
/* Perform averaging */
/* For each class, do: */
for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Processing instances of class %c...\n"
, classKRH[ii].classKRH);
#endif
/* Get the number of duplicates of the current class */
for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
break;
}
}
#ifdef DEBUG_PLUGIN_REORDER
printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
#endif
if(instances > 1) {
#ifdef DEBUG_PLUGIN_REORDER
printf(" Averaging %d contiguous epochs...\n", instances);
#endif
divisor = (double)instances;
/* For each time point in the duplicated classes: */
for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
/* For each voxel: */
for(kk = 0; kk < nvox; kk++) {
/* Average the duplicate time-points for the current voxel */
for(mm = 0, sum = 0.0; mm < instances; mm++) {
sum = sum + (double)sData[jj+(mm*classKRH[ii].length)][kk];
}
sData[timePoint][kk] = (short)((sum / divisor) + 0.5);
}
}
/* Update final length of time-course */
newLength = newLength - (classKRH[ii].length * (instances - 1));
/* Jump over the indices of all instances of the current class */
startIndex = startIndex + (instances *classKRH[ii].length);
/* Go to the next class to process */
ii += instances;
}
else { /* just copy the data points */
#ifdef DEBUG_PLUGIN_REORDER
printf(" Copying single epoch...\n");
#endif
for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
for(kk = 0; kk < nvox; kk++) {
sData[timePoint][kk] = sData[jj][kk];
}
}
/* Jump to the index of the next class in the time-course */
startIndex = startIndex + classKRH[ii].length;
/* eat up another class */
++ii;
}
/* Update progress meter with the percentage done */
perc = (100 * timePoint) / mapLength;
PLUTO_set_meter(plint, perc);
}
/* Set progress meter to 100% */
PLUTO_set_meter(plint, 100);
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Final time-length is %d...\n", newLength);
#endif
ii = EDIT_dset_items(
new_dset,
ADN_nvals, newLength, /* # time points */
ADN_ntt, newLength, /* # time points */
ADN_none);
if(0 != ii) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(sPtr);
FREEUP(sData);
FREE_WORKSPACE;
return "******************************************\n"
"!! Error while creating output dataset !!\n"
"*****************************************";
}
/* Free unused subBRIKS */
if(newLength < mapLength) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
#endif
for(ii = newLength; ii < mapLength; ii++) {
FREEUP(sData[ii]);
}
sData = (short **)realloc((void *)sData, sizeof(short *)*newLength);
if(NULL == sData) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(sPtr);
FREE_WORKSPACE;
return "*************************\n"
"!! Reallocation error !!\n"
"************************";
}
}
mapLength = newLength;
}
}
}
/* Write results into new dataset */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
#endif
for(kk = 0; kk < mapLength; kk++) {
EDIT_substitute_brick(new_dset, kk, MRI_short, sData[kk]);
}
break;
case MRI_float:
#ifdef DEBUG_PLUGIN_REORDER
printf("\n[Reorder] Collating/float data...\n");
#endif
/* Reorder the time-course of each voxel */
for(kk = 0; kk < mapLength; kk++) {
if(NULL == memcpy((void *)fData[kk]
, (void *)fPtr[map[kk]]
, sizeof(float)*nvox)) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(fPtr);
FREEUP(fData);
FREE_WORKSPACE;
return "***********************************\n"
"!! Error performing memory copy !!\n"
"**********************************";
}
perc = (100 * kk) / mapLength; /* display percentage done */
PLUTO_set_meter(plint, perc); /* on the progress meter */
}
PLUTO_set_meter(plint, 100); /* set progress meter to 100% */
if(dupaction) { /* new_dset may change in length again */
/* Verify that averaging is possible:
All duplicates must have the same length.
*/
error = 0; /* start with no error */
for(ii = 0; ii < classCount; ii++) {
for(kk = 1; kk < (classCount-1); kk++) {
if(classKRH[ii].classKRH == classKRH[kk].classKRH) {
if(classKRH[ii].length != classKRH[kk].length) {
error = 1; /* lengths differ! */
break;
}
dupsToAverage = True; /* there is at least one duplicate to average */
}
}
if(error) {
PLUTO_popup_message(plint
,"*********************************************************\n"
"* Duplicate stimulus regions in map file must have the *\n"
"* same length for averaging. Returning collated regions.*\n"
"*********************************************************");
break;
}
}
if(!error) { /* average the crap */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] All duplicate epochs have the same length...\n");
#endif
if(dupsToAverage) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Averaging duplicate epochs...\n");
printf("[Reorder] Initial time-length is %d...\n", mapLength);
#endif
newLength = mapLength; /* length may change after averaging, track new length */
/* Perform averaging */
/* For each class, do: */
for(ii = 0, startIndex = 0, timePoint = 0; ii < classCount;) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Processing instances of class %c...\n"
, classKRH[ii].classKRH);
#endif
/* Get the number of duplicates of the current class */
for(jj = ii + 1, instances = 1; jj < classCount; jj++, instances++) {
if(classKRH[ii].classKRH != classKRH[jj].classKRH) {
break;
}
}
#ifdef DEBUG_PLUGIN_REORDER
printf(" %d instances of class %c...\n", instances, classKRH[ii].classKRH);
#endif
if(instances > 1) {
#ifdef DEBUG_PLUGIN_REORDER
printf(" Averaging %d contiguous epochs...\n", instances);
#endif
divisor = (double)instances;
/* For each time point in the duplicated classes: */
for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
/* For each voxel: */
for(kk = 0; kk < nvox; kk++) {
/* Average the duplicate time-points for the current voxel */
for(mm = 0, sum = 0.0; mm < instances; mm++) {
sum = sum + (double)fData[jj+(mm*classKRH[ii].length)][kk];
}
fData[timePoint][kk] = (float)((sum / divisor) + 0.5);
}
}
/* Update final length of time-course */
newLength = newLength - (classKRH[ii].length * (instances - 1));
/* Jump over the indices of all instances of the current class */
startIndex = startIndex + (instances *classKRH[ii].length);
/* Go to the next class to process */
ii += instances;
}
else { /* just copy the data points */
#ifdef DEBUG_PLUGIN_REORDER
printf(" Copying single epoch...\n");
#endif
for(jj = startIndex; jj < (classKRH[ii].length+startIndex); jj++, timePoint++) {
for(kk = 0; kk < nvox; kk++) {
fData[timePoint][kk] = fData[jj][kk];
}
}
/* Jump to the index of the next class in the time-course */
startIndex = startIndex + classKRH[ii].length;
/* eat up another class */
++ii;
}
/* Update progress meter with the percentage done */
perc = (100 * timePoint) / mapLength;
PLUTO_set_meter(plint, perc);
}
/* Set progress meter to 100% */
PLUTO_set_meter(plint, 100);
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Final time-length is %d...\n", newLength);
#endif
ii = EDIT_dset_items(
new_dset,
ADN_nvals, newLength, /* # time points */
ADN_ntt, newLength, /* # time points */
ADN_none);
if(0 != ii) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(fPtr);
FREEUP(fData);
FREE_WORKSPACE;
return "******************************************\n"
"!! Error while creating output dataset !!\n"
"*****************************************";
}
/* Free unused subBRIKS */
if(newLength < mapLength) {
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Freeing %d unused subBRIKs...\n", mapLength-newLength);
#endif
for(ii = newLength; ii < mapLength; ii++) {
FREEUP(fData[ii]);
}
fData = (float **)realloc((void *)fData, sizeof(float *)*newLength);
if(NULL == fData) {
THD_delete_3dim_dataset(new_dset, False);
FREEUP(fPtr);
FREE_WORKSPACE;
return "*************************\n"
"!! Reallocation error !!\n"
"************************";
}
}
mapLength = newLength;
}
}
}
/* Write results into new dataset */
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Writing %d subBRIKS into new data set...\n", mapLength);
#endif
for(kk = 0; kk < mapLength; kk++) {
EDIT_substitute_brick(new_dset, kk, MRI_float, fData[kk]);
}
break;
} /* end of switch over input type */
DSET_unload(old_dset);
/*************** Cleanup and go home *****************/
PLUTO_add_dset(plint, new_dset, DSET_ACTION_MAKE_CURRENT);
#ifdef DEBUG_PLUGIN_REORDER
printf("[Reorder] Cleaning up workspace...\n");
#endif
FREE_WORKSPACE;
return NULL; /* NULL string returned means all was OK */
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1