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