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 #include "afni.h"
00008 #include "afni_plugin.h"
00009 
00010 #ifndef ALLOW_PLUGINS
00011 #  error "Plugins not properly set up -- see machdep.h"
00012 #endif
00013 
00014 /* -----------------------START---------------------------------*/
00015 /* definition and decleration part to suport the main algorithm */
00016 
00017 #include <stdlib.h>
00018 #include <stdio.h>
00019 #include <math.h>
00020 
00021 
00022 /*------------------ macros to return workspace at exit -------------------*/
00023 
00024 #define ZFREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
00025 
00026 #undef  ZFREE_WORKSPACE
00027 #define ZFREE_WORKSPACE                              \
00028   do{ ZFREEUP(bptr) ; ZFREEUP(sptr) ; ZFREEUP(fptr) ;  \
00029       ZFREEUP(cptr) ; ZFREEUP(fxar) ; ZFREEUP(fac)  ;  \
00030       ZFREEUP(dtr)  ;   \
00031     } while(0)
00032 /*-------------------------------------------------------------------------*/
00033 
00034 
00035 /*      Definitions of prototypes and declaration of support functions 
00036         this is taken from the list of include files that I use in the original code*/ 
00037 
00038 
00039 /*-------------------------------------------------------------------*/
00040 /* COMPLEX STRUCTURE */
00041 
00042 typedef struct {
00043     float x, y, z;
00044 } fXYZ;
00045 
00046 
00047 /***********************************************************************
00048   Plugin to extract 3D+time time courses whos index or xyz corrdinates 
00049   match a certain criterion
00050 ************************************************************************/
00051 typedef struct 
00052         {
00053                   int nxx;                      /* number of voxels in the x direction */
00054                   int nyy;                      /* number of voxels in the y direction */
00055                   int nzz;                      /* number of voxels in the z direction */
00056                   char *dsetname; /* prefix of data set analyzed */
00057                   int ignore;                   /* number ofpoints to ignore from time courses */
00058                   int ln;                       /* length of FMRI vector */
00059                   int dtrnd;            /* remove linear trend or just the mean */
00060                   int errcode;          /* error code number returned from function */
00061                   int out;                      /* flag for writing to a file */
00062                   int outts;            /* flag for writing time series to a file */
00063                   int format;
00064                   int iloc;
00065                   int xloc;
00066                   int yloc;
00067                   int zloc;
00068                   int ncols;
00069                   int nrows;
00070                   float * indvect;      /* vector that will hold the list of indices */
00071                   fXYZ * xyzvect;                       /* vecor that will hold the list of xyz coordinates */
00072                   char * strout;
00073                   char * strin;
00074                   FILE * outwritets;
00075                   FILE * outlogfile;
00076                   char outname[PLUGIN_MAX_STRING_RANGE]; /* output data file name */
00077         }extract_data;
00078 
00079 /*--------------------- string to 'help' the user --------------------*/
00080 
00081 static char helpstring[] =
00082   "                               4D Dump Plugin\n\n"
00083   "This plugin allows the extraction of selected FMRI time series into an ascii file. \n"
00084   "The time series to be extracted are specified either by their AFNI index or \n"
00085   "by their AFNI XYZ corrdinates.\n\n"
00086   "Plugin Inputs\n"
00087   "   1- Data :\n"
00088   "      3D+time    -> Chooser for 3D+time data set.\n"
00089   "      Ignore     -> Specify the number of points to ignore from the beginning \n"
00090   "                    of each voxel time series.\n"
00091   "      Dtrnd      -> (Y/N) Apply detrending to time series. \n"
00092   "                    If the dtrnd option is off, time series do not have zero mean.\n\n"
00093   "   2- Mask file :\n"
00094   "      Mask File  -> The index or X Y Z coordinates that specify wich voxel time series\n"
00095   "                    should be output are specified in an ascii file.\n"
00096   "                    Each line in the ascii file holds the index or X Y Z coordinates \n"
00097   "                    of the voxel time series to be extracted. Each line of the ascii \n"
00098   "                    file can hold many values, the user chooses which ones correspond \n"
00099   "                    to i or the X Y Z triplet. The X Y Z coordinates used for masking\n"
00100   "                    correspond to those shown in the AFNI window and NOT the graph window. \n"
00101   "                    You can see these coordinates in the upper left side of the AFNI window.\n"
00102   "                    To do so, you must switch the voxel coordinate units from mm to slice coordinates.\n"
00103   "                      Define Datamode -> Misc -> Voxel Coords ?\n"
00104   "                    Masking files could be the ascii output files of other plugins such as\n" 
00105   "                   '3Ddump' and 'Hilber Delay98'.\n"
00106   "      N Columns  -> Total number of values on each line in the ascii mask file.\n\n"
00107   "   3- Index Mask : (optional)\n"
00108   "      i col.     -> Column number where the AFNI indices for the time series to be extracted \n"
00109   "                    are stored in the ascii mask file.\n"
00110   "   4- XYZ Mask : (optional)\n"
00111   "      x,y,z col. -> Specify the column numbers where X Y Z triplets are stored\n\n"
00112   "   5- Output :\n"
00113   "      Filename   -> Name of the ascii file where the time series are written.\n"
00114   "                    If no output filename is specified, the ouput is the prefix\n"
00115   "                    of the 3D+time input brick with the extenstion '.4Ddump' appended to it.\n"
00116   "                    A LOG file, 'Filename.log' is also written to disk.\n"
00117   "                    The log file contains all the parameters settings used\n"
00118   "                    for generating 'Filename'.\n"
00119   "      Format     -> ts[1] ts[2] .... -> one time series per line.\n"
00120   "                    i x y z ts[1] ts[2] .... -> index, X Y Z coordinates and time series per line.\n\n"
00121   "   PS: If all the voxels specified in the mask file exist in the data file then each line\n"
00122   "       in the output file correspond to the line in the mask file. Otherwise you should use\n"
00123   "       the option i x y z ts[1] ts[2] .... to know which time series came from which voxel\n\n"
00124   "If you have/find questions/comments/bugs about the plugin, \n"
00125   "send me an E-mail: ziad@image.bien.mu.edu\n\n"
00126   "                          Ziad Saad Jan 6 97, lastest update Feb 23 98.\n\n"
00127 ;
00128 
00129 /*--------------------- strings for output format --------------------*/
00130 
00131 static char * yn_strings[] = { "n" , "y" }; 
00132 static char * format_strings[] = { "i x y z ts[1] ..." , "ts[1] ts[2] ..." };
00133 
00134 #define NUM_YN_STRINGS (sizeof(yn_strings)/sizeof(char *))
00135 #define NUM_FORMAT_STRINGS (sizeof(yn_strings)/sizeof(char *))
00136 
00137 
00138 #define YUP  1
00139 #define NOPE 0
00140 
00141 #define ERROR_NOSUCHVOXEL       1                               /* Nothing to do in function */
00142 #define ERROR_FILEREAD          2
00143 #define ERROR_FILEWRITE         2
00144 #define ERROR_OPTIONS           17
00145 #define ERROR_OUTCONFLICT       19
00146 
00147 /*----------------- prototypes for internal routines -----------------*/
00148 
00149 static char * EXTRACT_main( PLUGIN_interface * ) ;  /* the entry point */
00150 
00151 static void EXTRACT_tsfunc( double T0 , double TR ,
00152                             int npts , float ts[] , double ts_mean , double ts_slope ,
00153                             void * udp) ;
00154 
00155 static void show_ud (extract_data* ud);
00156 
00157 static void write_ud (extract_data* ud);
00158 
00159 static void indexTOxyz (extract_data* ud, int ncall, int *xpos , int *ypos , int *zpos);
00160 
00161 static void error_report (extract_data* ud, int ncall );
00162 
00163 static void writets (extract_data* ud,float * ts, int ncall);
00164 
00165 static float * extract_index (char *fname, int ind_col_loc, int ncols, int *nrows, int *Err);
00166 
00167 static fXYZ * extract_xyz (char *fname, int x_col_loc, int y_col_loc, int z_col_loc, int ncols, int *nrows, int *Err);
00168 
00169 static void disp_vect (float *v,int l);
00170 
00171 static int filexists (char *f_name);
00172 
00173 static int f_file_size (char *f_name);
00174 
00175 static int * PLUTO_4D_to_nothing (THD_3dim_dataset * old_dset , int ignore , int detrend ,
00176                          generic_func * user_func , void * user_data );
00177 
00178 /*---------------------------- global data ---------------------------*/
00179 
00180 static PLUGIN_interface * global_plint = NULL ;
00181 
00182 /***********************************************************************
00183    Set up the interface to the user:
00184     1) Create a new interface using "PLUTO_new_interface";
00185 
00186     2) For each line of inputs, create the line with "PLUTO_add_option"
00187          (this line of inputs can be optional or mandatory);
00188 
00189     3) For each item on the line, create the item with
00190         "PLUTO_add_dataset" for a dataset chooser,
00191         "PLUTO_add_string"  for a string chooser,
00192         "PLUTO_add_number"  for a number chooser.
00193 ************************************************************************/
00194 
00195 DEFINE_PLUGIN_PROTOTYPE
00196 
00197 PLUGIN_interface * PLUGIN_init( int ncall )
00198 {
00199    PLUGIN_interface * plint ;     /* will be the output of this routine */
00200 
00201    if( ncall > 0 ) return NULL ;  /* only one interface */
00202 
00203    /*---------------- set titles and call point ----------------*/
00204 
00205    plint = PLUTO_new_interface( "4D Dump" ,
00206                                 "Extract voxel time courses given their index or XYZ coordinates" ,
00207                                 helpstring ,
00208                                 PLUGIN_CALL_VIA_MENU , EXTRACT_main  ) ;
00209 
00210    global_plint = plint ;  /* make global copy */
00211 
00212    /*--------- 1st line: Input dataset and mask files ---------*/
00213 
00214    PLUTO_add_option( plint ,
00215                      "Data" ,  /* label at left of input line */
00216                      "Data" ,  /* tag to return to plugin */
00217                      TRUE       /* is this mandatory? */
00218                    ) ;
00219 
00220    PLUTO_add_dataset(  plint ,
00221                        "3D+time" ,        /* label next to button   */
00222                        ANAT_ALL_MASK ,    /* take only EPI datasets */
00223                        FUNC_ALL_MASK ,          /*  No fim funcs   */
00224                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00225                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00226                     ) ;
00227   
00228    PLUTO_add_number( plint ,
00229                     "Ignore" ,  /* label next to chooser */
00230                     0 ,         /* smallest possible value */
00231                     50 ,        /* largest possible value (inactivated for now)*/
00232                     0 ,         /* decimal shift (none in this case) */
00233                     0 ,         /* default value */
00234                     FALSE       /* allow user to edit value? */
00235                   ) ;
00236         PLUTO_add_string( plint ,
00237                      "Dtrnd" ,  /* label next to textfield */
00238                      2,yn_strings,    /* strings to choose among */
00239                      1         /* Default option */
00240                    ) ;
00241         
00242         /*---------- 2nd line: Mask file info  ----------*/
00243    PLUTO_add_option( plint ,
00244                      "Mask file" ,  /* label at left of input line */
00245                      "Mask" ,  /* tag to return to plugin */
00246                      TRUE       /* is this mandatory? */
00247                    ) ;
00248    
00249    PLUTO_add_string( plint , "Mask File" , 0 , NULL , 19 ) ;
00250    
00251    PLUTO_add_number( plint ,
00252                     "N Columns" ,  /* label next to chooser */
00253                     1 ,         /* smallest possible value */
00254                     1000 ,        /* largest possible value (inactivated for now)*/
00255                     0 ,         /* decimal shift (none in this case) */
00256                     3 ,         /* default value */
00257                     TRUE       /* allow user to edit value? */
00258                   ) ;
00259    
00260 
00261    /*---------- 3rd line: index mask location ----------*/
00262    
00263    PLUTO_add_option( plint ,
00264                      "Index Mask ?" ,  /* label at left of input line */
00265                      "Index" ,  /* tag to return to plugin */
00266                      FALSE       /* is this mandatory? */
00267                    ) ;
00268    
00269    PLUTO_add_number( plint ,
00270                     "i col." ,  /* label next to chooser */
00271                     1 ,         /* smallest possible value */
00272                     1000 ,        /* largest possible value (inactivated for now)*/
00273                     0 ,         /* decimal shift (none in this case) */
00274                     1 ,         /* default value */
00275                     TRUE       /* allow user to edit value? */
00276                   ) ;
00277         
00278                    
00279    /*---------- 4th line: xyz mask location ----------*/
00280 
00281    PLUTO_add_option( plint ,
00282                      "XYZ Mask ?" ,  /* label at left of input line */
00283                      "XYZ" ,  /* tag to return to plugin */
00284                      FALSE       /* is this mandatory? */
00285                    ) ;
00286 
00287    PLUTO_add_number( plint ,
00288                     "x col." ,  /* label next to chooser */
00289                     1 ,         /* smallest possible value */
00290                     1000 ,        /* largest possible value */
00291                     0 ,         /* decimal shift (none in this case) */
00292                     2 ,         /* default value */
00293                     TRUE       /* allow user to edit value? */
00294                   ) ;
00295         
00296         PLUTO_add_number( plint ,
00297                     "y col." ,  /* label next to chooser */
00298                     1 ,         /* smallest possible value */
00299                     1000 ,        /* largest possible value */
00300                     0 ,         /* decimal shift (none in this case) */
00301                     3 ,         /* default value */
00302                     TRUE       /* allow user to edit value? */
00303                   ) ;
00304                   
00305 
00306         PLUTO_add_number( plint ,
00307                     "z col." ,  /* label next to chooser */
00308                     1 ,         /* smallest possible value */
00309                     1000 ,      /* largest possible value */
00310                     0 ,         /* decimal shift (none in this case) */
00311                     4 ,         /* default value */
00312                     TRUE       /* allow user to edit value? */
00313                   ) ;
00314                   
00315    /*---------- 5th line: output stuff ----------*/
00316 
00317    PLUTO_add_option( plint ,
00318                      "Output" ,  /* label at left of input line */
00319                      "Output" ,  /* tag to return to plugin */
00320                      TRUE        /* is this mandatory? */
00321                    ) ;
00322 
00323    
00324    PLUTO_add_string( plint , "Filename" , 0 , NULL , 19 ) ;
00325               
00326    PLUTO_add_string( plint ,
00327                      "Format" ,  /* label next to textfield */
00328                      2,format_strings,    /* strings to choose among */
00329                      0          /* Default option */
00330                    ) ;
00331    
00332    return plint ;
00333 }
00334 
00335 /***************************************************************************
00336   Main routine for this plugin (will be called from AFNI).
00337   If the return string is not NULL, some error transpired, and
00338   AFNI will popup the return string in a message box.
00339 ****************************************************************************/
00340 
00341 static char * EXTRACT_main( PLUGIN_interface * plint )
00342 {
00343    extract_data uda,*ud;
00344    MRI_IMAGE * tsim;
00345    MCW_idcode * idc ;                          /* input dataset idcode */
00346    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00347    char *tmpstr , * str ;                 
00348    int   ntime, nvec , Err , itmp, nprfx;
00349         float * vec , fs , T ;
00350         char * tag;                     /* plugin option tag */ 
00351         
00352         /* Allocate as much character space as Bob specifies in afni.h + a bit more */
00353         
00354         tmpstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
00355         
00356         if (tmpstr == NULL) 
00357                                                                           return "********************\n"
00358                                                                                                 "Could not Allocate\n"
00359                                                                                                 "a teeni weeni bit of\n"
00360                                                                                                 "Memory ! \n"
00361                                                                                                 "********************\n";
00362                                                                                                 
00363         ud = &uda;              /* ud now points to an allocated space */
00364         ud->errcode = 0;        /*reset error flag */
00365         
00366    /*--------------------------------------------------------------------*/
00367    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00368 
00369    /*--------- go to first input line ---------*/
00370                 
00371    tag = PLUTO_get_optiontag(plint) ;
00372    
00373    if (tag == NULL)
00374         {
00375                 return "************************\n"
00376              "Bad 1st line option \n"
00377              "************************"  ;
00378         }       
00379 
00380    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00381    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00382    if( old_dset == NULL )
00383       return "*************************\n"
00384              "Cannot find Input Dataset\n"
00385              "*************************"  ;
00386    
00387    ud->dsetname = DSET_PREFIX (old_dset);
00388         
00389         ud->ignore = PLUTO_get_number(plint) ;    /* get number item */
00390         
00391         str = PLUTO_get_string(plint) ;                                 
00392         ud->dtrnd = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
00393         
00394         
00395         
00396         /*--------- loop over remaining options ---------*/
00397         
00398                 
00399         ud->iloc = -1;
00400         ud->xloc = -1;
00401         ud->yloc = -1;
00402         ud->zloc = -1;
00403         do
00404                 {
00405                         tag = PLUTO_get_optiontag(plint) ;
00406                         if (tag == NULL) break;
00407                         if (strcmp (tag, "Mask") == 0)
00408                                 {
00409                                         ud->strin = PLUTO_get_string(plint) ; 
00410                                         ud->ncols = PLUTO_get_number(plint) ;
00411                                         continue;
00412                                 }
00413                         
00414                         if (strcmp (tag, "Index") == 0)
00415                                 {
00416                                         ud->iloc = PLUTO_get_number(plint) ;    /* get number item */
00417                                         continue;
00418                                 }
00419                 
00420                         if (strcmp (tag, "XYZ") == 0)
00421                                 {
00422                                                 ud->xloc = PLUTO_get_number(plint) ;    /* get number item */
00423                                                 ud->yloc = PLUTO_get_number(plint) ;    /* get number item */
00424                                                 ud->zloc = PLUTO_get_number(plint) ;    /* get number item */
00425                                                 continue;
00426                                 }
00427 
00428                         if (strcmp (tag, "Output") == 0)
00429                                         { 
00430                                         ud->strout = PLUTO_get_string(plint) ; 
00431 
00432                                         str = PLUTO_get_string(plint) ;                                 
00433                                                 ud->format = (int)PLUTO_string_index( str , NUM_FORMAT_STRINGS , format_strings );
00434                                                 continue;
00435                                         }
00436                         
00437                 } while (1);
00438         /* ------------------ Check for some errorsor inconsistencies ------------- */
00439                 
00440         if (ud->iloc == -1 && ud->xloc == -1)
00441                 {
00442                         return "**************************\n"
00443                                          "At least iloc or x/y/zloc\n"
00444                                          "must be specified\n"
00445                                          "**************************\n"
00446                                          ;
00447                 }
00448         
00449         if (ud->iloc != -1 && ud->xloc != -1)
00450                 {
00451                         return "***************************\n"
00452                                          "iloc AND x/y/zloc can not\n"
00453                                          "be simultaneously specified\n"
00454                                          "***************************\n"
00455                                          ;
00456                 }
00457         
00458         
00459         /* ------------------Done with user parameters ---------------------------- */
00460         
00461         /* Now loadup that index list or the xyz list */
00462         if (ud->iloc != -1)
00463                 {       
00464                         itmp = 0;  /* might want to give option of setting it to number of rows if*/ 
00465                     /* the users know it, otherwise, it is automatically determined*/    
00466                         ud->indvect = extract_index (ud->strin, ud->iloc, ud->ncols, &itmp, &Err);
00467                 }
00468         else            /* assuming the only other case is that x y z are specified */
00469                 {
00470                         itmp = 0; 
00471                         ud->xyzvect = extract_xyz (ud->strin , ud->xloc , ud->yloc , ud->zloc , ud->ncols, &itmp, &Err);
00472                 }
00473         
00474                 
00475         ud->nrows = itmp;
00476         
00477         switch (Err)
00478         {
00479                 case (0):
00480                         break;
00481                 case (1):
00482                         return "****************************************\n"
00483                                "index location should be > 0 and < ncols\n"
00484                                "****************************************\n";
00485                 case (2):
00486                         return "********************************************\n"
00487                 "file size and number of columns do not match\n"
00488                                "********************************************\n";
00489                 case (3):
00490                         return "**********************\n"
00491                 "Can't find matrix file\n"
00492                                "**********************\n";
00493                 case (4):
00494                         return "*****************\n"
00495                 "ncols must be > 0\n"
00496                                "*****************\n";
00497                 case (5):
00498                         return "****************************************\n"
00499                 "x/y/z column numbers can NOT be the same\n"
00500                                "****************************************\n";
00501                 default:
00502                         return "****************************************\n"
00503                 "Should not have gotten here .....\n"
00504                                "****************************************\n";
00505                 
00506         }
00507         
00508         /* check to see if the field is empty */
00509         if (ud->strout == NULL) nprfx = 0;
00510                 else nprfx = 1;
00511         
00512         /* check if the size is larger than 0. I did not want to check for this unless it's allocated */
00513         if (nprfx == 1 && (int)strlen (ud->strout) == 0) nprfx = 0;
00514                 
00515                 
00516         if (nprfx == 0)   /* no output file is specified */ 
00517                 {
00518                         sprintf ( tmpstr , "%s.4Ddump" , DSET_PREFIX (old_dset));
00519                         ud->strout = tmpstr;
00520                 }
00521         
00522         
00523         if (filexists(ud->strout) == 1)
00524                 {
00525                         return "*******************************\n"
00526                                          "Outfile exists, won't overwrite\n"
00527                                          "*******************************\n";   
00528                 }
00529         ud->outwritets = fopen (ud->strout ,"w");       
00530         
00531         sprintf ( tmpstr , "%s.log" , ud->strout);
00532         if (filexists(tmpstr) == 1)
00533                 {
00534                         return "*******************************\n"
00535                                          "Outfile.log exists, won't overwrite\n"
00536                                          "*******************************\n";   
00537                 }
00538         
00539         ud->outlogfile = fopen (tmpstr ,"w"); 
00540                 
00541         if ((ud->outwritets == NULL) || (ud->outlogfile == NULL) )
00542                                                 {
00543                                                         ud->errcode = ERROR_FILEWRITE; 
00544                                                         
00545                                                         return "***********************\n"
00546                                                                          "Could Not Write Outfile\n"
00547                                                                          "***********************\n";
00548                                                 }                               
00549         
00550         
00551         ud->nxx = (int)old_dset->daxes->nxx;                            /* get data set dimensions */
00552         ud->nyy = (int)old_dset->daxes->nyy;
00553         ud->nzz = (int)old_dset->daxes->nzz;
00554         
00555    /* ready to dump the log file */
00556    write_ud (ud);
00557    
00558    /*------------- ready to compute new dataset -----------*/
00559 
00560         PLUTO_4D_to_nothing (old_dset , ud->ignore , ud->dtrnd ,
00561                              EXTRACT_tsfunc , (void *)ud );
00562 
00563         fclose (ud->outlogfile);
00564         fclose (ud->outwritets);
00565         
00566         if (tmpstr) free (tmpstr);              
00567    return NULL ;  /* null string returned means all was OK */
00568 }
00569 
00570 /**********************************************************************
00571    Function that does the real work
00572 ***********************************************************************/
00573 
00574 static void EXTRACT_tsfunc( double T0 , double TR ,
00575                    int npts , float ts[] , double ts_mean , double ts_slope ,
00576                    void * udp)
00577 {
00578    static int nvox = -1, ncall = -1;
00579         extract_data uda,*ud;
00580         float xcor=0.0 ,  tmp=0.0 , tmp2 = 0.0 ,  dtx = 0.0 ,\
00581                          slp = 0.0 , vts = 0.0 , vrvec = 0.0 , rxcorCoef = 0.0;
00582         int i , is_ts_null , status , opt , actv , zpos , ypos , xpos , hit;
00583         
00584         ud = &uda;
00585         ud = (extract_data *) udp;
00586         
00587         
00588    /** is this a "notification"? **/
00589 
00590         
00591    if( ncall == -1 || ncall == nvox)
00592         {
00593       if( npts > 0 ){  /* the "start notification" */
00594          PLUTO_popup_meter( global_plint ) ;  /* progress meter  */
00595          nvox  = npts ;                       /* keep track of   */
00596          ncall = 0 ;                          /* number of calls */
00597                         
00598       } else {  /* the "end notification" */
00599          PLUTO_set_meter( global_plint , 100 ) ; /* set meter to 100% */
00600 
00601       }
00602       return ;
00603    }
00604 
00605         
00606         hit = 0;
00607         ud->ln = npts;
00608         
00609         if (ud->iloc != -1)                     /* for each ncall, find out if this index is wanted*/
00610                 {
00611                         for (i=0;i<ud->nrows;++i)
00612                                 {
00613                                         if (ud->indvect[i] == (float)ncall) 
00614                                                 {
00615                                                         hit = 1;
00616                                                         writets (ud,ts,ncall); 
00617                                                         i = ud->nrows;
00618                                                 }
00619                                 }
00620                 }
00621         else
00622                 {
00623                         indexTOxyz (ud, ncall, &xpos , &ypos , &zpos);
00624                         for (i=0;i<ud->nrows;++i)
00625                                 {
00626                                         if (ud->xyzvect[i].x == (float)xpos)
00627                                                 {
00628                                                         if (ud->xyzvect[i].y == (float)ypos)
00629                                                                 {
00630                                                                         if (ud->xyzvect[i].z == (float)zpos)
00631                                                                                 {
00632                                                                                         hit = 1;
00633                                                                                         writets (ud,ts,ncall); 
00634                                                                                         i = ud->nrows;
00635                                                                                 }
00636                                                                 }
00637                                                 }               
00638                                 }
00639                         
00640                 }
00641         
00642         
00643    if (ud->errcode == 0)                                /* if there are no errors, proceed */   
00644                 {/*                                             */
00645         }/* ud->errcode == 0 outer loop */
00646    
00647    /** set the progress meter to the % of completion **/
00648    ncall++ ;
00649    
00650    PLUTO_set_meter( global_plint , (100*ncall)/nvox ) ;
00651    
00652    ud->errcode = 0;     /* Reset error to nothing */
00653    
00654    return ;
00655 }
00656 
00657 /* ************************************************************ */
00658 /* function to display user data input (debugging stuff)        */
00659 /* ************************************************************ */
00660 
00661 static void show_ud (extract_data* ud)
00662         {
00663                 printf ("\n\nUser Data Values at location :\n");
00664                 printf ("ud->dsetname= %s\n",ud->dsetname);
00665                 printf ("ud->strin= %s\n",ud->strin);
00666                 printf ("ud->strout= %s\n",ud->strout);
00667                 printf ("ud->nxx= %d\n",ud->nxx);
00668                 printf ("ud->nyy= %d\n",ud->nyy);
00669                 printf ("ud->nzz= %d\n",ud->nzz);
00670                 printf ("ud->iloc= %d\n",ud->iloc);
00671                 printf ("ud->xloc= %d\n",ud->xloc);
00672                 printf ("ud->yloc= %d\n",ud->yloc);
00673                 printf ("ud->zloc= %d\n",ud->zloc);
00674                 printf ("ud->ncols= %d\n",ud->ncols);
00675                 printf ("ud->nrows= %d\n",ud->nrows);
00676                 printf ("ud->ignore= %d\n",ud->ignore);
00677                 printf ("ud->dtrnd= %d\n",ud->dtrnd);
00678                 printf ("ud->format= %d\n",ud->format);
00679                 printf ("Hit enter to continue\a\n\n");
00680                 getchar ();
00681                 return;
00682         }
00683 
00684 /* ************************************************************ */
00685 /* function to write user data input to log file        */
00686 /* ************************************************************ */
00687 
00688 static void write_ud (extract_data* ud)
00689         {
00690                 fprintf (ud->outlogfile,"\n\nUser Data Values \n");
00691                 fprintf (ud->outlogfile,"ud->dsetname= %s\n",ud->dsetname);
00692                 fprintf (ud->outlogfile,"ud->strin= %s\n",ud->strin);
00693                 fprintf (ud->outlogfile,"ud->strout= %s\n",ud->strout);
00694                 fprintf (ud->outlogfile,"ud->nxx= %d\n",ud->nxx);
00695                 fprintf (ud->outlogfile,"ud->nyy= %d\n",ud->nyy);
00696                 fprintf (ud->outlogfile,"ud->nzz= %d\n",ud->nzz);
00697                 fprintf (ud->outlogfile,"ud->iloc= %d\n",ud->iloc);
00698                 fprintf (ud->outlogfile,"ud->xloc= %d\n",ud->xloc);
00699                 fprintf (ud->outlogfile,"ud->yloc= %d\n",ud->yloc);
00700                 fprintf (ud->outlogfile,"ud->zloc= %d\n",ud->zloc);
00701                 fprintf (ud->outlogfile,"ud->ncols= %d\n",ud->ncols);
00702                 fprintf (ud->outlogfile,"ud->nrows= %d\n",ud->nrows);
00703                 fprintf (ud->outlogfile,"ud->ignore= %d\n",ud->ignore);
00704                 fprintf (ud->outlogfile,"ud->dtrnd= %d\n",ud->dtrnd);
00705                 fprintf (ud->outlogfile,"ud->format= %d\n",ud->format);
00706                 fprintf (ud->outlogfile,"\nThe format for the output file is the following:\n");
00707                 switch (ud->format)
00708                         {
00709                                 case (0):
00710                                         fprintf (ud->outlogfile,"ncall\txpos\typos\tzpos\tts[0]\tts[1]\t...\n");
00711                                         break;
00712                                 case (1):
00713                                         fprintf (ud->outlogfile,"ts[0]\tts[1]\t...\n");
00714                                         break;
00715                                         
00716                         }
00717                 
00718                 return;
00719         }
00720 
00721 /* ************************************************************ */ 
00722 /* function to compute x, y, z coordinates from the index       */
00723 /* ************************************************************ */ 
00724 
00725 static void indexTOxyz (extract_data* ud, int ncall, int *xpos , int *ypos , int *zpos)         
00726         {
00727                 *zpos = (int)ncall / (int)(ud->nxx*ud->nyy);
00728                 *ypos = (int)(ncall - *zpos * ud->nxx * ud->nyy) / (int)ud->nxx;
00729                 *xpos = ncall - ( *ypos * ud->nxx ) - ( *zpos * ud->nxx * ud->nyy ) ;
00730                 return;
00731         }
00732         
00733 /* ************************************************************ */
00734 /* function to report errors encountered to the logfile         */
00735 /* Only errors that happen during runtime  are handeled here,   */
00736 /* the others are handeled instantaneously, and need not be     */
00737 /* logged                                                                                                                                                */
00738 /* ************************************************************ */
00739 
00740 void error_report (extract_data* ud, int ncall ) 
00741         {
00742                 int xpos,ypos,zpos;
00743                 
00744                 indexTOxyz (ud, ncall, &xpos , &ypos , &zpos); 
00745                 
00746                 switch (ud->errcode)
00747                         {
00748                                 
00749                                 default:
00750                                         fprintf (ud->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",ud->errcode);
00751                                         break;
00752                         }       
00753                 fprintf (ud->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos  );
00754                 return;
00755         }
00756         
00757 /* *************************************************************** */
00758 /* function to write the time course into a line in the given file */
00759 /* *************************************************************** */
00760 
00761 void writets (extract_data * ud,float * ts, int ncall)
00762 
00763         {       
00764                 int i,xpos,ypos,zpos;
00765                 
00766                 switch (ud->format)
00767                         {
00768                                 case (0):
00769                                         indexTOxyz (ud,ncall,&xpos , &ypos , &zpos); 
00770                                         fprintf (ud->outwritets, "%d\t%d\t%d\t%d\t",ncall,xpos, ypos,zpos);
00771                                         break;
00772                                 case (1):
00773                                         break;
00774                                 default:
00775                                         break;
00776                         }
00777                 
00778                 for (i=0;i<ud->ln;++i)
00779                         {
00780                                 fprintf (ud->outwritets, "%f\t",ts[i]);
00781                         }
00782                 fprintf (ud->outwritets,"\n");
00783         }
00784 
00785 /* *************************************************************** */
00786 /* function to extract x y z colums form a matrix format  file */
00787 /* *************************************************************** */
00788 static fXYZ * extract_xyz (char *fname, int x_col_loc, int y_col_loc, int z_col_loc, int ncols, int *nrows, int *Err)
00789 {/*extract_xyz*/
00790         
00791         float tmp, *tmpX;
00792         int sz,i,indx,tst;
00793         div_t divstuff,tempx,tempy,tempz;
00794         FILE * INFILE;
00795         fXYZ * xyzvect=NULL ;
00796         
00797         /* ncols must be > 0 */
00798         if (ncols <= 0)
00799                 
00800                 {
00801                         *Err = 4;
00802                         return (xyzvect); /* ncols <= 0 !*/
00803                 }
00804                 
00805                 
00806         /* ind_col_loc must be > 0  (1st column is 1, and it must be less than ncols) */
00807         if (x_col_loc <= 0 || x_col_loc > ncols || y_col_loc <= 0 || y_col_loc > ncols || z_col_loc <= 0 || z_col_loc > ncols )
00808                 {
00809                         *Err = 1;
00810                         return (xyzvect); /* ?_col_loc should be >0 and <ncols*/
00811                 }
00812         
00813         /* if the number of rows is given, compute required size */
00814         if (*nrows > 0)
00815                 {
00816                         sz = *nrows * ncols;
00817                 }
00818         else
00819                 { /* dtermine the size and compute the number of rows, check if it's an integer*/
00820                         sz = f_file_size (fname);
00821                         if (sz == -1)
00822            {
00823                    *Err = 3;
00824                    return (xyzvect);
00825             }
00826                         divstuff = div (sz,ncols);
00827                         if (divstuff.rem != 0)
00828                                 {
00829                                         *Err = 2;
00830                                         return (xyzvect); /* size of file and number of columns don't match */
00831                                 }
00832                         else 
00833                                 {
00834                                         *nrows = divstuff.quot;
00835                                 }
00836                 }
00837         
00838         tst = (x_col_loc - y_col_loc) * (x_col_loc - z_col_loc) * (y_col_loc - z_col_loc);
00839         if (tst == 0)
00840                 {
00841                         *Err = 5;
00842                         return (xyzvect); /* columns specified are the same */
00843                 } 
00844         
00845         /* Allocate and check for necessary space */
00846         xyzvect = (fXYZ *) calloc (sz+2,sizeof(fXYZ));
00847         
00848          if (xyzvect == NULL)
00849                                 {
00850                                         printf ("\nFatal Error : Failed to Allocate memory\a\n");
00851                                         printf ("Abandon Lab Immediately !\n\n");
00852                                         return NULL ;
00853                                 };
00854                                 
00855         INFILE = fopen (fname,"r");
00856         if (INFILE == NULL)
00857                 {
00858                         *Err = 3; /* can't find file */
00859                         return (xyzvect);
00860                 }
00861 
00862         tempx = div (x_col_loc,ncols);
00863         tempy = div (y_col_loc,ncols);
00864         tempz = div (z_col_loc,ncols);
00865         
00866         for (i=0;i<sz;++i)
00867                 {
00868                         fscanf (INFILE,"%f",&tmp);
00869                         divstuff = div ((i+1),ncols);
00870                         if (divstuff.rem != 0)
00871                                 indx = divstuff.quot;
00872                         else 
00873                                 indx = divstuff.quot - 1;
00874                         
00875                         if (divstuff.rem == tempx.rem)
00876                                         xyzvect[indx].x = tmp;
00877                                 
00878                                 else if (divstuff.rem == tempy.rem)
00879                                                 xyzvect[indx].y = tmp;
00880                                         
00881                                         else if (divstuff.rem == tempz.rem)
00882                                                         xyzvect[indx].z = tmp;
00883                 }
00884 
00885         
00886         *Err = 0;
00887         return (xyzvect);
00888         
00889 }/*extract_xyz*/
00890 
00891 
00892 /* *************************************************************** */
00893 /* function to extract an indices columnform a matrix format   file */
00894 /* *************************************************************** */
00895 
00896 
00897 static float * extract_index (char *fname, int ind_col_loc, int ncols, int *nrows, int *Err)
00898 {/*extract_index*/
00899         
00900         float tmp, *indxvect=NULL;
00901         int sz,i;
00902         div_t divstuff,temp;
00903         FILE * INFILE;
00904         
00905         /* ncols must be > 0 */
00906         if (ncols <= 0)
00907                 
00908                 {
00909                         *Err = 4;
00910                         return (indxvect); /* ncols <= 0 !*/
00911                 }
00912                 
00913         
00914         /* ind_col_loc must be > 0  (1st column is 1, and it must be less than ncols) */
00915         if (ind_col_loc <= 0 || ind_col_loc > ncols)
00916                 {
00917                         *Err = 1;
00918                         return (indxvect); /* ind_col_loc should be >0 and <ncols*/
00919                 }
00920         
00921         /* if the number of rows is given, compute required size */
00922         if (*nrows > 0)
00923                 {
00924                         sz = *nrows * ncols;
00925                 }
00926         else
00927                 { /* dtermine the size and compute the number of rows, check if it's an integer*/
00928                         sz = f_file_size (fname);
00929                         if (sz == -1)
00930             {
00931                *Err = 3;
00932               return (indxvect);      
00933             }
00934                         divstuff = div (sz,ncols);
00935                         if (divstuff.rem != 0)
00936                                 {
00937                                         *Err = 2;
00938                                         return (indxvect); /* size of file and number of columns don't match */
00939                                 }
00940                         else 
00941                                 {
00942                                         *nrows = divstuff.quot;
00943                                 }
00944                 }
00945         
00946         /* Allocate and check for necessary space */
00947         indxvect = (float *) calloc (sz+2,sizeof(float));
00948         
00949          if (indxvect == NULL)
00950                                 {
00951                                         printf ("\nFatal Error : Failed to Allocate memory\a\n");
00952                                         printf ("Abandon Lab Immediately !\n\n");
00953                                         return NULL ;
00954                                 }; 
00955         
00956         INFILE = fopen (fname,"r");
00957         if (INFILE == NULL)
00958                 {
00959                         *Err = 3; /* can't find file */
00960                         return (indxvect);
00961                 }
00962 
00963         temp = div (ind_col_loc,ncols);
00964         
00965         for (i=0;i<sz;++i)
00966                 {
00967                         fscanf (INFILE,"%f",&tmp);
00968                         divstuff = div ((i+1),ncols);
00969                         if (divstuff.rem == temp.rem)
00970                                 {
00971                                         if (divstuff.rem != 0) indxvect[divstuff.quot] = tmp;
00972                                                 else indxvect[divstuff.quot-1] = tmp;
00973                                 }
00974                 }
00975 
00976         
00977         *Err = 0;
00978         return (indxvect);
00979         
00980 }/*extract_index*/
00981 
00982 /* *************************************************************** */
00983 /* function to return the size of a float numbers   file */
00984 /* *************************************************************** */
00985 
00986 static int f_file_size (char *f_name)
00987    
00988     { 
00989       
00990 
00991      int cnt=0,ex;
00992      float buf;
00993      
00994      FILE*internal_file;
00995      
00996      internal_file = fopen (f_name,"r");
00997      if (internal_file == NULL) {
00998                                                                 return (-1);
00999                                                         }
01000      ex = fscanf (internal_file,"%f",&buf);                                             
01001      while (ex != EOF)
01002       {
01003         ++cnt;
01004         ex = fscanf (internal_file,"%f",&buf);
01005       }
01006       
01007       
01008       fclose (internal_file);
01009       return (cnt);                                                          
01010    }
01011  
01012 /* *************************************************************** */
01013 /* function to test if a file exists                                                                        */
01014 /* *************************************************************** */
01015 
01016 static int filexists (char *f_name)
01017 {/*filexists*/
01018         FILE *outfile;
01019         
01020         outfile = fopen (f_name,"r");
01021         if (outfile == NULL)
01022                 return (0);
01023         else 
01024                 fclose (outfile);
01025                 return (1);
01026                 
01027 }/*filexists*/
01028 
01029 
01030 /* *************************************************************** */
01031 /* function to display a vector (debugging)                                                              */
01032 /* *************************************************************** */
01033 
01034 static void disp_vect (float *v,int l)
01035         {
01036                 int i;
01037 
01038                 printf ("\n");
01039                 if ((l-1) == 0)
01040                         {
01041                                 printf ("V = %f\n",*v);
01042                         }
01043                 else 
01044                 {
01045                         for (i=0;i<l;++i)
01046                         {
01047                                 printf ("V[%d] = %f\t",i,v[i]);
01048                         }
01049                         printf ("\n");
01050                 }
01051                 return;
01052 
01053         }
01054 
01055 /* *************************************************************** */
01056 /* function to retrieve 4D data from AFNI                          */
01057 /*    (adapted from BDW's function PLUTO_4D_to_typed_fith                        */
01058 /* *************************************************************** */
01059 
01060 static int * PLUTO_4D_to_nothing (THD_3dim_dataset * old_dset , int ignore , int detrend ,
01061                          generic_func * user_func, void * user_data )
01062 {
01063 
01064    byte    ** bptr = NULL ;  /* one of these will be the array of */
01065    short   ** sptr = NULL ;  /* pointers to input dataset sub-bricks */
01066    float   ** fptr = NULL ;  /* (depending on input datum type) */
01067    complex ** cptr = NULL ;
01068 
01069    float * fxar = NULL ;  /* array loaded from input dataset */
01070    float * fac  = NULL ;  /* array of brick scaling factors */
01071    float * dtr  = NULL ;  /* will be array of detrending coeff */
01072 
01073    float val , d0fac , d1fac , x0,x1;
01074    double tzero , tdelta , ts_mean , ts_slope ;
01075    int   ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox ;
01076    static int retval;
01077         register int kk ;
01078 
01079    /*----------------------------------------------------------*/
01080    /*----- Check inputs to see if they are reasonable-ish -----*/
01081 
01082    if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;
01083 
01084    if( user_func == NULL ) return NULL ;
01085 
01086    if( ignore < 0 ) ignore = 0 ;
01087 
01088    /*--------- set up pointers to each sub-brick in the input dataset ---------*/
01089 
01090    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;   /* get old dataset datum */
01091    nuse      = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */
01092    if( nuse < 2 ) return NULL ;
01093 
01094    DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */
01095 
01096    kk = THD_count_databricks( old_dset->dblk ) ;  /* check if it was */
01097    if( kk < DSET_NVALS(old_dset) ){               /* loaded correctly */
01098       DSET_unload( old_dset ) ;
01099       return NULL ;
01100    }
01101 
01102    switch( old_datum ){  /* pointer type depends on input datum type */
01103 
01104       default:                      /** don't know what to do **/
01105          DSET_unload( old_dset ) ;
01106          return NULL ;
01107 
01108       /** create array of pointers into old dataset sub-bricks **/
01109 
01110       /*--------- input is bytes ----------*/
01111       /* voxel #i at time #k is bptr[k][i] */
01112       /* for i=0..nvox-1 and k=0..nuse-1.  */
01113 
01114       case MRI_byte:
01115          bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
01116          if( bptr == NULL ) return NULL ;
01117          for( kk=0 ; kk < nuse ; kk++ )
01118             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
01119       break ;
01120 
01121       /*--------- input is shorts ---------*/
01122       /* voxel #i at time #k is sptr[k][i] */
01123       /* for i=0..nvox-1 and k=0..nuse-1.  */
01124 
01125       case MRI_short:
01126          sptr = (short **) malloc( sizeof(short *) * nuse ) ;
01127          if( sptr == NULL ) return NULL ;
01128          for( kk=0 ; kk < nuse ; kk++ )
01129             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
01130       break ;
01131 
01132       /*--------- input is floats ---------*/
01133       /* voxel #i at time #k is fptr[k][i] */
01134       /* for i=0..nvox-1 and k=0..nuse-1.  */
01135 
01136       case MRI_float:
01137          fptr = (float **) malloc( sizeof(float *) * nuse ) ;
01138          if( fptr == NULL ) return NULL ;
01139          for( kk=0 ; kk < nuse ; kk++ )
01140             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
01141       break ;
01142 
01143       /*--------- input is complex ---------*/
01144       /* voxel #i at time #k is cptr[k][i]  */
01145       /* for i=0..nvox-1 and k=0..nuse-1.   */
01146 
01147       case MRI_complex:
01148          cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;
01149          if( cptr == NULL ) return NULL ;
01150          for( kk=0 ; kk < nuse ; kk++ )
01151             cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;
01152       break ;
01153 
01154    } /* end of switch on input type */
01155 
01156         nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
01157 
01158    
01159    /*---- allocate space for 1 voxel timeseries ----*/
01160 
01161    fxar = (float *) malloc( sizeof(float) * nuse ) ;   /* voxel timeseries */
01162    if( fxar == NULL ){ ZFREE_WORKSPACE ; return NULL ; }
01163 
01164    /*--- get scaling factors for sub-bricks ---*/
01165 
01166    fac = (float *) malloc( sizeof(float) * nuse ) ;   /* factors */
01167    if( fac == NULL ){ ZFREE_WORKSPACE ; return NULL ; }
01168 
01169    use_fac = 0 ;
01170    for( kk=0 ; kk < nuse ; kk++ ){
01171       fac[kk] = DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
01172       if( fac[kk] != 0.0 ) use_fac++ ;
01173       else                 fac[kk] = 1.0 ;
01174    }
01175    if( !use_fac ) ZFREEUP(fac) ;
01176 
01177    /*--- setup for detrending ---*/
01178 
01179    dtr = (float *) malloc( sizeof(float) * nuse ) ;
01180    if( dtr == NULL ){ ZFREE_WORKSPACE ; return NULL ; }
01181 
01182    d0fac = 1.0 / nuse ;
01183    d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
01184    for( kk=0 ; kk < nuse ; kk++ )
01185       dtr[kk] = kk - 0.5 * (nuse-1) ;  /* linear trend, orthogonal to 1 */
01186 
01187 
01188    /*----- set up to find time at each voxel -----*/
01189 
01190    tdelta = old_dset->taxis->ttdel ;
01191    if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
01192    if( tdelta == 0.0 ) tdelta = 1.0 ;
01193 
01194    izold  = -666 ;
01195    nxy    = old_dset->daxes->nxx * old_dset->daxes->nyy ;
01196 
01197    /*----------------------------------------------------*/
01198    /*----- Setup has ended.  Now do some real work. -----*/
01199 
01200    /* start notification */
01201 #if 0
01202    user_func(  0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data ) ;
01203 #else
01204    { void (*uf)(double,double,int,float *,double,double,void *) =
01205      (void (*)(double,double,int,float *,double,double,void *))(user_func) ;
01206      uf( 0.0l,0.0l , nvox , NULL , 0.0l,0.0l , user_data ) ;
01207    }
01208 #endif
01209 
01210    /***** loop over voxels *****/   
01211    for( ii=0 ; ii < nvox ; ii++  ){  /* 1 time series at a time */
01212                 
01213       /*** load data from input dataset, depending on type ***/
01214 
01215       switch( old_datum ){
01216 
01217          /*** input = bytes ***/
01218 
01219          case MRI_byte:
01220             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
01221          break ;
01222 
01223          /*** input = shorts ***/
01224 
01225          case MRI_short:
01226             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
01227          break ;
01228 
01229          /*** input = floats ***/
01230 
01231          case MRI_float:
01232             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
01233          break ;
01234 
01235          /*** input = complex (note we use absolute value) ***/
01236 
01237          case MRI_complex:
01238             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
01239          break ;
01240 
01241       } /* end of switch over input type */
01242 
01243       /*** scale? ***/
01244      if( use_fac )
01245          for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
01246 
01247       /** compute mean and slope **/
01248 
01249       x0 = x1 = 0.0 ;
01250       for( kk=0 ; kk < nuse ; kk++ ){
01251          x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
01252       }
01253 
01254       x0 *= d0fac ; x1 *= d1fac ;  /* factors to remove mean and trend */
01255 
01256       ts_mean  = x0 ;
01257       ts_slope = x1 / tdelta ;
01258  
01259       /** detrend? **/
01260 
01261       if( detrend )
01262          for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
01263 
01264       /** compute start time of this timeseries **/
01265                 /* The info computed here is not being used in this version*/
01266       iz = ii / nxy ;    /* which slice am I in? */
01267 
01268       if( iz != izold ){          /* in a new slice? */
01269          tzero = THD_timeof( ignore ,
01270                              old_dset->daxes->zzorg
01271                            + iz*old_dset->daxes->zzdel , old_dset->taxis ) ;
01272          izold = iz ;
01273 
01274          if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
01275       }
01276 
01277       /*** Send data to user function ***/
01278 #if 0
01279       user_func( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data) ;
01280 #else
01281      { void (*uf)(double,double,int,float *,double,double,void *) =
01282        (void (*)(double,double,int,float *,double,double,void *))(user_func) ;
01283        uf( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data) ;
01284      }
01285 #endif
01286 
01287       
01288 
01289    } /* end of outer loop over 1 voxels at a time */
01290 
01291    DSET_unload( old_dset ) ;  
01292 
01293    /* end notification */
01294 #if 0
01295    user_func( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data ) ;
01296 #else
01297    { void (*uf)(double,double,int,float *,double,double,void *) =
01298      (void (*)(double,double,int,float *,double,double,void *))(user_func) ;
01299      uf( 0.0l,0.0l, 0 , NULL,0.0l,0.0l, user_data ) ;
01300    }
01301 #endif
01302 
01303    
01304    /*-------------- Cleanup and go home ----------------*/
01305    
01306    ZFREE_WORKSPACE ;
01307         retval = 0;
01308         return &retval; /* this value is not used for now .... */
01309 
01310 }
01311