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 /*      Definitions of prototypes and declaration of support functions 
00022         this is taken from the list of include files that I use in the original code*/ 
00023 
00024 /*-------------------------------------------------------------------*/
00025 /* taken from #include "/usr/people/ziad/Programs/C/DSP_in_C/dft.h" */
00026 /* dft.h - function prototypes and structures for dft and fft functions */
00027 
00028 
00029 
00030 
00031 /* COMPLEX STRUCTURE */
00032 
00033 typedef struct {
00034     float real, imag;
00035 } COMPLEX;
00036 
00037 /* function prototypes for dft and inverse dft functions */
00038 
00039 static void fft(COMPLEX *,int);
00040 static void ifft(COMPLEX *,int);
00041 static void dft(COMPLEX *,COMPLEX *,int);
00042 static void idft(COMPLEX *,COMPLEX *,int);
00043 static void rfft(float *,COMPLEX *,int);
00044 static void ham(COMPLEX *,int);
00045 static void han(COMPLEX *,int);
00046 static void triang(COMPLEX *,int);
00047 static void black(COMPLEX *,int);
00048 static void harris(COMPLEX *,int);
00049 
00050 #include "plug_delay_V2.h"
00051 
00052 
00053 /***********************************************************************
00054   Plugin to compute a 3D+time dataset voxelwise delay with respect to 
00055   a reference waveform
00056 ************************************************************************/
00057 typedef struct 
00058         {
00059                   int nxx;                      /* number of voxels in the x direction */
00060                   int nyy;                      /* number of voxels in the y direction */
00061                   int nzz;                      /* number of voxels in the z direction */
00062                   char *dsetname; /* prefix of data set analyzed */
00063                   char *refname;        /* reference vector filename */
00064                   float *rvec;          /* reference vetor */
00065                   float fs;                     /* Sampling frequency */
00066                                                                 /* it is only used for the log file in this version*/
00067                                                                 /* the ts_func, gives TR automatically */
00068                   float T;                      /* Stimulus period */
00069                   float co;                     /* Correlation Coefficient Threshold*/
00070                   float dmask;  /* delays to be masked will take this value ( Not used anymore )*/
00071                   int unt;                      /* Delay units */
00072                   int wrp;                      /* flag for Polar Wrap */
00073                   int Navg;                     /* number of data sets averaged to obtain the brick (for statistical stuff) */
00074                   int Nort;             /* Number of nuisance parameters (orts) (for statistical stuff) */
00075                   int Nfit;                     /* Number of fit parameters (for statistical stuff) */
00076                   int Nseg;                     /* Number of segments */
00077                   int nsamp;            /* number of points in time series */
00078                   int ignore;           /* number ofpoints to ignore from time courses */
00079                   int Pover;            /* Percent overlap */
00080                   int ln;                       /* length of FMRI vector */
00081                   int dtrnd;            /* remove linear trend or just the mean */
00082                   int biasrem;          /* flag for removing delay bias */
00083                   int Dsamp;            /* flag for correction of non uniform sampling start time */
00084                   int errcode;          /* error code number returned from hdelay */
00085                   int out;                      /* flag for writing delays to a file */
00086                   int outts;            /* flag for writing time series to a file */
00087                   char * new_prefix; /* new prefix for data set */
00088                   char * strout;
00089                   FILE * outwrite;
00090                   FILE * outwritets;
00091                   FILE * outlogfile;
00092                   char outname[PLUGIN_MAX_STRING_RANGE]; /* output data file name */
00093         }hilbert_data_V2;
00094 
00095 /*--------------------- string to 'help' the user --------------------*/
00096 
00097 static char helpstring[] =
00098   "            Hilbert Delay98 Plugin \n\n"
00099   "The Plugin estimates the time delay between each voxel time series in a 3D+time dataset \n"
00100   "and a reference time series[1][2]. The estimated delays are relative to the reference time series.\n"
00101   "For example, a delay of 4 seconds means that the voxel time series is delayed by 4 seconds with respect\n"
00102   "to the reference time series.\n\n"
00103   "Plugin Inputs:\n\n"
00104   "   1- Data : \n"
00105   "      3d+time    -> 3D+time dataset to analyze.\n"
00106   "      Nort       -> Number of orts or nuisance parameters. \n"
00107   "                    It is set to 2 by default because the mean \n"
00108   "                    and linear trends are always removed from the time series.\n"
00109   "                    This value is used for estimating the p-value at the cross\n"
00110   "                    correlation threshold selected with AFNI's threshold slider.\n"
00111   "                    The p-value is estimated only when the cross correlation coefficient\n"
00112   "                    subbrick is used for thresholding.\n\n"
00113   "   2- Ref. :\n"
00114   "      Ref. Vect. -> 1D Reference time series vector. \n"
00115   "                    The reference time series vector is stored in an ascii file.\n"
00116   "                    The plugin assumes that there is one value per line and that all\n"
00117   "                    values in the file are part of the reference vector.\n"
00118   "                    PS: Unlike with 3dfim, and FIM in AFNI, values over 33333 are treated \n"
00119   "                        as part of the time series.\n" 
00120   "                    The reference vectors must be placed in a directory specified in \n"
00121   "                    AFNI_TSPATH environment variable. Read AFNI documentation for more info.\n"
00122   "      Ignore     -> Number of samples to ignore from the beginning of each voxel's time series.\n"
00123   "                    Ignore is NOT applied to the reference time series.\n"
00124   "      Dsamp      -> (Y/N) Correct a voxel's estimated delay by the time at which the slice\n"
00125   "                    containing that voxel was acquired.\n\n"
00126   "   3- Sig. :\n"
00127   "      fs in Hz   -> Sampling frequency in Hz. of data time series. \n"  
00128   "      Tstim sec  -> Stimulus period in seconds. \n"
00129   "                    If the stimulus is not periodic, you can set Tstim to 0.\n"
00130   "      C-Off      -> Cross Correlation Coefficient threshold value.\n"
00131   "                    This is only used to censor the ascii output (see below).\n"
00132   "      No-bias    -> (y/n) Correct for the bias in the cross correlation coefficient estimate [1][2].\n\n" 
00133   "   4- Alg. :\n"
00134   "      N seg.     -> Number of segments used to estimate the periodogram.\n"
00135   "                    (you can't modify this parameter in this version)\n"
00136   "      % ovrlp    -> Percent segment overlap when estimating periodogram.\n"
00137   "                    (you can't modify this parameter in this version)\n" 
00138   "      Units      -> (Seconds/Degrees/Radians) Units for delay estimates.\n"
00139   "                    You can't use Degrees or Radians as units unless you specify\n"
00140   "                    a value for Tstim > 0.\n"
00141   "      Phz Wrp    -> (Y/N) Delay (or phase) wrap.\n"
00142   "                    This switch maps delays from: \n" 
00143   "                    (Seconds) 0->T/2 to 0->T/2 and T/2->T to -T/2->0\n"  
00144   "                    (Degrees) 0->180 to 0->180 and 180->360 to -180->0\n"   
00145   "                    (Radians) 0->pi to 0->pi and pi->2pi to -pi->0\n" 
00146   "                    You can't use this option unless you specify a value for Tstim > 0.\n\n"
00147   "   5- Output :\n"
00148   "      AFNI Prfx  -> Prefix for output brick of the bucket type (fbuc).\n"
00149   "                    The first subbrick is for Delay.\n"
00150   "                    The second subbrick is for Covariance, which is an estimate\n" 
00151   "                    of the power in voxel time series at the frequencies present \n"
00152   "                    in the reference time series.\n"
00153   "                    The third subbrick is for the Cross Correlation Coefficients between\n"
00154   "                    FMRI time series and reference time series.\n"
00155   "                    The fourth subbrick contains estimates of the Variance of voxel time series.\n"
00156   "                    If you leave the field empty, a default prefix is used.\n"
00157   "                    The default prefix is the prefix of the input 3D+time brick \n"
00158   "                    with a '.DEL' extension appended to it.\n"
00159   "      Write      -> (Y/N) Write the results to an ascii file for voxels with \n"
00160   "                    Cross Correlation Coefficients larger than C-Off.\n"
00161   "                    Each line in the output file contains information for one voxel.\n"
00162   "                    There a 9 columns in the file which hold the following values:\n"
00163   "                    1- Voxel Index (VI) : Each voxel in an AFNI brick has a unique index.\n"
00164   "                                          Indices map directly to XYZ coordinates.\n"
00165   "                                          See AFNI plugin documentations for more info.\n"
00166   "                    2..4- Voxel coordinates (X Y Z): Those are the voxel slice coordinates.\n"
00167   "                                                     You can see these coordinates in the upper left side\n"
00168   "                                                     of the AFNI window. To do so, you must first switch the\n"
00169   "                                                     voxel coordinate units from mm to slice coordinates.\n"
00170   "                                                     Define Datamode -> Misc -> Voxel Coords ?\n"
00171   "                                                     PS: The coords that show up in the graph window\n"
00172   "                                                     could be different from those in the upper left side \n"
00173   "                                                     of AFNI's main window.\n"
00174   "                    5- Duff : A value of no interest to you. It is preserved for it's historical value.\n"
00175   "                    6- Delay (Del) : The estimated voxel delay.\n"
00176   "                    7- Covariance (Cov) : Covariance estimate.\n"
00177   "                    8- Cross Correlation Coefficient (xCorCoef) : Cross Correlation Coefficient estimate.\n"
00178   "                    9- Variance (VTS) : Variance of voxel's time series.\n"
00179   "                    This output file can be used as an input to two other plugins:\n"
00180   "                    '4Ddump' and '3D+t Extract'\n\n"
00181   "                    If Write is used, a log file is written too.\n"
00182   "                    The log file contains all parameter settings used for generating the output brick.\n"
00183   "                    The name of the log file is the same as 'Filename' (see below) with a '.log' extension\n"
00184   "                    appended at the end. The log file also holds any warnings generated by the plugin.\n"
00185   "                    Some warnings, such as 'null time series ...' , or 'Could not find zero crossing ...'\n"
00186   "                    are harmless, I might remove them in future versions.\n"
00187   "      Filename   -> Name of the ascii file to write results to.\n"
00188   "                    If the field is left empty, a default name similar \n"
00189   "                    to the default output prefix is used.\n"
00190   "      Write ts   -> (Y/N) Write the time series to an ascii file for voxels with \n"
00191   "                    Cross Correlation Coefficients larger than C-Off.\n" 
00192   "                    The file name is that used in 'Filename' with a '.ts' extension appended at the end.\n"
00193   "                    A line (L) in the file 'Filename.ts' contains the time series of the voxel whose\n"
00194   "                    results are written on line (L) in the file 'Filename'.\n"
00195   "                    The time series written to 'Filename.ts' do not contain the ignored samples,\n"
00196   "                    they are detrended and have zero mean.\n\n"
00197   "Random Comments/Advice:\n"
00198   "Of course, the longer you time series, the better. It is generally recomended that\n"
00199   "the largest delay be less than N/10, N being the length of the time series.\n"
00200   "The algorithm does go all the way to N/2, but that's not too good.\n\n"
00201   "Disclaimer: \n"
00202   "(Blaaa bla bla bla bla .... --> legal terms you probably wouldn't read) \n"
00203   "             I am not responsible for anything bad.\n\n"
00204   "If you have/find questions/comments/bugs about the plugin, \n"
00205   "send me an E-mail: ziad@image.bien.mu.edu\n\n"
00206   "                          Ziad Saad June 20 97, lastest update Aug 21 00.\n\n"
00207   "[1] : Bendat, J. S. (1985). The Hilbert transform and applications to correlation measurements,\n"
00208   "       Bruel and Kjaer Instruments Inc.\n"
00209   "[2] : Bendat, J. S. and G. A. Piersol (1986). Random Data analysis and measurement procedures, \n"
00210   "       John Wiley & Sons.\n"
00211 ;
00212 
00213 /*--------------------- strings for output format --------------------*/
00214 /* do not change the order in this string*/
00215 static char * method_strings[] = { "Seconds" , "Degrees" , "Radians"} ;
00216 static char * yn_strings[] = { "n" , "y" }; 
00217 
00218 #define NUM_METHOD_STRINGS (sizeof(method_strings)/sizeof(char *))
00219 #define NUM_YN_STRINGS (sizeof(yn_strings)/sizeof(char *))
00220 
00221 /* do not change these three*/
00222 #define METH_SECONDS 0
00223 #define METH_DEGREES 1
00224 #define METH_RADIANS 2
00225 
00226 #undef  DELAY
00227 #define DELAY    0
00228 #define XCOR     1
00229 #define XCORCOEF 2
00230 #ifndef NOWAYXCORCOEF
00231         #define NOWAYXCORCOEF 0                                 /* A flag value indicating that something lethal went on */
00232 #endif
00233 
00234 #define NBUCKETS 4                              /* Number of values per voxel in Buket data set */
00235 #define DELINDX 0                                       /* index of delay value in results vector */
00236 #define COVINDX 1                                       /* index of covariance value in results vector */
00237 #define COFINDX 2                                       /* index of cross correlation coefficient value in results vector */
00238 #define VARINDX 3                                       /* index of FMRI time course variance value in results vector */
00239 
00240 #define YUP  1
00241 #define NOPE 0
00242 
00243 #define ERROR_NOTHINGTODO       1                               /* Nothing to do in hilbertdelay_V2 function */
00244 #define ERROR_LARGENSEG         2                               /* Too many segments specified in hilbertdelay_V2 function */
00245 #define ERROR_LONGDELAY         3                               /* Could not detect zero crossing before half of time course was gone */
00246 #define ERROR_WRONGUNIT         8                               /* Wrong units selected to pass to the delay functions */
00247 #define ERROR_WARPVALUES        9
00248 #define ERROR_FSVALUES          10
00249 #define ERROR_TVALUES           11
00250 #define ERROR_TaUNITVALUES      12
00251 #define ERROR_TaWRAPVALUES      13
00252 #define ERROR_FILEOPEN          15
00253 #define ERROR_SERIESLENGTH      16
00254 #define ERROR_OPTIONS           17
00255 #define ERROR_NULLTIMESERIES    18
00256 #define ERROR_OUTCONFLICT       19
00257 #define ERROR_BADLENGTH         20
00258 
00259 /*----------------- prototypes for internal routines -----------------*/
00260 
00261 static char * DELAY_main( PLUGIN_interface * ) ;  /* the entry point */
00262 
00263 static void DELAY_tsfuncV2( double T0 , double TR ,
00264                    int npts , float ts[] , double ts_mean , double ts_slope ,
00265                    void * udp , int nbrick , float * buckar) ;
00266 
00267 static void show_ud (hilbert_data_V2* ud,int loc);
00268 
00269 static void write_ud (hilbert_data_V2* ud);
00270 
00271 static void indexTOxyz (hilbert_data_V2* ud, int ncall, int *xpos , int *ypos , int *zpos);
00272 
00273 static void error_report (hilbert_data_V2* ud, int ncall );
00274 
00275 static void writets (hilbert_data_V2* ud,float * ts);
00276 
00277 /*---------------------------- global data ---------------------------*/
00278 
00279 static PLUGIN_interface * global_plint = NULL ;
00280 
00281 /***********************************************************************
00282    Set up the interface to the user:
00283     1) Create a new interface using "PLUTO_new_interface";
00284 
00285     2) For each line of inputs, create the line with "PLUTO_add_option"
00286          (this line of inputs can be optional or mandatory);
00287 
00288     3) For each item on the line, create the item with
00289         "PLUTO_add_dataset" for a dataset chooser,
00290         "PLUTO_add_string"  for a string chooser,
00291         "PLUTO_add_number"  for a number chooser.
00292 ************************************************************************/
00293 
00294 
00295 DEFINE_PLUGIN_PROTOTYPE
00296 
00297 PLUGIN_interface * PLUGIN_init( int ncall )
00298 {
00299    PLUGIN_interface * plint ;     /* will be the output of this routine */
00300 
00301    if( ncall > 0 ) return NULL ;  /* only one interface */
00302 
00303    /*---------------- set titles and call point ----------------*/
00304 
00305    plint = PLUTO_new_interface( "Hilbert Delay98" ,
00306                                 "Time delay between FMRI and reference time series" ,
00307                                 helpstring ,
00308                                 PLUGIN_CALL_VIA_MENU , DELAY_main  ) ;
00309 
00310    global_plint = plint ;  /* make global copy */
00311 
00312    /*--------- 1st line: Input dataset ---------*/
00313 
00314    PLUTO_add_option( plint ,
00315                      "Data" ,  /* label at left of input line */
00316                      "Data" ,  /* tag to return to plugin */
00317                      TRUE       /* is this mandatory? */
00318                    ) ;
00319 
00320    PLUTO_add_dataset(  plint ,
00321                        "3D+time" ,        /* label next to button   */
00322                        ANAT_ALL_MASK ,    /* take only EPI datasets */
00323                        FUNC_ALL_MASK ,    /*  No fim funcs   */
00324                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00325                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00326                     ) ;
00327                                                  
00328         PLUTO_add_number( plint ,
00329                     "Nort" ,  /* label next to chooser */
00330                     1 ,         /* smallest possible value */
00331                     100 ,        /* largest possible value (inactivated for now)*/
00332                     0 ,         /* decimal shift (none in this case) */
00333                     2 ,         /* default value */
00334                     FALSE       /* allow user to edit value? */
00335                   ) ;
00336         
00337    /*---------- 2nd line: Input time series ----------*/
00338    
00339    PLUTO_add_option( plint ,
00340                      "Ref." ,  /* label at left of input line */
00341                      "Ref." ,  /* tag to return to plugin */
00342                      TRUE       /* is this mandatory? */
00343                    ) ;
00344 
00345    PLUTO_add_timeseries(plint,"Ref. Vect."); 
00346    
00347    PLUTO_add_number( plint ,
00348                     "Ignore" ,  /* label next to chooser */
00349                     0 ,         /* smallest possible value */
00350                     50 ,        /* largest possible value (inactivated for now)*/
00351                     0 ,         /* decimal shift (none in this case) */
00352                     0 ,         /* default value */
00353                     FALSE       /* allow user to edit value? */
00354                   ) ;
00355         
00356         PLUTO_add_string( plint ,
00357                      "Dsamp" ,  /*label next to textfield */
00358                      2,yn_strings,  /*   strings to choose among */
00359                      1          /* Default option */
00360                    ) ; 
00361                    
00362    /*---------- 3rd line: sampling frequency ----------*/
00363 
00364    PLUTO_add_option( plint ,
00365                      "Sig." ,  /* label at left of input line */
00366                      "Sig." ,  /* tag to return to plugin */
00367                      TRUE       /* is this mandatory? */
00368                    ) ;
00369 
00370    PLUTO_add_number( plint ,
00371                     "fs in Hz" ,  /* label next to chooser */
00372                     0 ,         /* smallest possible value */
00373                     2000 ,        /* largest possible value */
00374                     1 ,         /* decimal shift (none in this case) */
00375                     5 ,         /* default value */
00376                     TRUE       /* allow user to edit value? */
00377                   ) ;
00378         
00379         PLUTO_add_number( plint ,
00380                     "Tstim sec" ,  /* label next to chooser */
00381                     0.0 ,         /* smallest possible value */
00382                     500 ,        /* largest possible value */
00383                     0 ,         /* decimal shift (none in this case) */
00384                     40 ,         /* default value */
00385                     TRUE       /* allow user to edit value? */
00386                   ) ;
00387 
00388         PLUTO_add_number( plint ,
00389                     "C-Off" ,  /* label next to chooser */
00390                     -10 ,         /* smallest possible value */
00391                     10 ,        /* largest possible value */
00392                     1 ,         /* decimal shift  */
00393                     5 ,         /* default value */
00394                     TRUE       /* allow user to edit value? */
00395                   ) ;
00396    
00397    
00398    PLUTO_add_string( plint ,
00399                      "No-bias" ,  /*label next to textfield */
00400                      2,yn_strings,  /*   strings to choose among */
00401                      1          /* Default option */
00402                    ) ; 
00403                   
00404 
00405         
00406    /*---------- 4th line: Delay Units ----------*/
00407 
00408    PLUTO_add_option( plint ,
00409                      "Alg." ,  /* label at left of input line */
00410                      "Alg." ,  /* tag to return to plugin */
00411                      TRUE        /* is this mandatory? */
00412                    ) ;
00413 
00414    PLUTO_add_number( plint ,
00415                     "N seg." ,  /* label next to chooser */
00416                     1 ,         /* smallest possible value */
00417                     1 ,        /* largest possible value (turned Off for the moment, supporting code is present)*/
00418                     0 ,         /* decimal shift (none in this case) */
00419                     1 ,         /* default value */
00420                     FALSE       /* allow user to edit value? */
00421                   ) ;
00422         
00423         PLUTO_add_number( plint ,
00424                     "% ovrlp" ,  /* label next to chooser */
00425                     0 ,         /* smallest possible value */
00426                     0 ,        /* largest possible value (not implemented)*/
00427                     0 ,         /* decimal shift (none in this case) */
00428                     0 ,         /* default value */
00429                     FALSE       /* allow user to edit value? */
00430                   ) ;
00431 
00432         
00433    PLUTO_add_string( plint ,
00434                      "Units" ,  /* label next to textfield */
00435                      3,method_strings,    /* strings to choose among */
00436                      0          /* Default option */
00437                    ) ;
00438    
00439    PLUTO_add_string( plint ,
00440                      "Phz Wrp" ,  /* label next to textfield */
00441                      2,yn_strings,    /* strings to choose among */
00442                      0          /* Default option */
00443                    ) ;
00444                   
00445 
00446    /*---------- 5th line: Output dataset ----------*/
00447 
00448    PLUTO_add_option( plint ,
00449                      "Output" ,  /* label at left of input line */
00450                      "Output" ,  /* tag to return to plugin */
00451                      TRUE        /* is this mandatory? */
00452                    ) ;
00453 
00454    PLUTO_add_string( plint ,
00455                      "AFNI Prfx" ,  /* label next to textfield */
00456                      0,NULL ,    /* no fixed strings to choose among */
00457                      19          /* 19 spaces for typing in value */
00458                    ) ;
00459         
00460         PLUTO_add_string( plint ,
00461                      "Write" ,  /* label next to textfield */
00462                      2,yn_strings ,    
00463                      1          
00464                    ) ;
00465                    
00466    PLUTO_add_string( plint , "Filename" , 0 , NULL , 19 ) ;
00467    
00468    PLUTO_add_string( plint ,
00469                      "Write ts" ,  /* label next to textfield */
00470                      2,yn_strings ,    
00471                      1          
00472                    ) ;
00473 
00474    /*--------- done with interface setup ---------*/
00475    return plint ;
00476 }
00477 
00478 /***************************************************************************
00479   Main routine for this plugin (will be called from AFNI).
00480   If the return string is not NULL, some error transpired, and
00481   AFNI will popup the return string in a message box.
00482 ****************************************************************************/
00483 
00484 static char * DELAY_main( PLUGIN_interface * plint )
00485 {
00486    hilbert_data_V2 uda,*ud;
00487    MRI_IMAGE * tsim;
00488    MCW_idcode * idc ;                          /* input dataset idcode */
00489    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00490    char *tmpstr , * str , *nprfxstr;                 /* strings from user */
00491    int   ntime, nvec ,nprfx, i;
00492         float * vec , fs , T ;
00493                 
00494         /* Allocate as much character space as Bob specifies in afni.h + a bit more */
00495         
00496         tmpstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
00497         nprfxstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
00498         
00499         if (tmpstr == NULL || nprfxstr == NULL) 
00500                                                                           return "********************\n"
00501                                                                                                 "Could not Allocate\n"
00502                                                                                                 "a teeni weeni bit of\n"
00503                                                                                                 "Memory ! \n"
00504                                                                                                 "********************\n";
00505                                                                                                                 
00506         ud = &uda;              /* ud now points to an allocated space */
00507         ud->errcode = 0;        /*reset error flag */
00508         
00509    /*--------------------------------------------------------------------*/
00510    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00511 
00512    /*--------- go to first input line ---------*/
00513                 
00514    PLUTO_next_option(plint) ;
00515 
00516    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
00517    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
00518    if( old_dset == NULL )
00519       return "*************************\n"
00520              "Cannot find Input Dataset\n"
00521              "*************************"  ;
00522    
00523    ud->dsetname = DSET_FILECODE (old_dset);
00524         ud->nsamp = DSET_NUM_TIMES (old_dset);
00525         ud->Navg = 1 ;    /* Navg does not play a role for the p value, averaging increases sensitivity */
00526         ud->Nort = PLUTO_get_number(plint) ; /* Should be two by default, for mean and linear trend */
00527         ud->Nfit = 2 ;  /* Always 2 for phase and amplitude for this plugin */
00528         /*--------- go to 2nd input line, input time series ---------*/
00529                 
00530         PLUTO_next_option(plint) ;
00531         
00532         tsim = PLUTO_get_timeseries(plint);
00533         if (tsim == NULL) return "No Timeseries Input";
00534         
00535         ud->ln = (int)tsim -> nx;                                                                       /* number of points in each vector */
00536         nvec    = tsim -> ny;                                                                   /* number of vectors */
00537         ud->rvec   = (float *) MRI_FLOAT_PTR(tsim);     /* vec[i+j*nx] = ith point of jth vector */
00538                                                                                                                                 /* for i=0 .. ntime-1 and j=0 .. nvec-1 */
00539         
00540         if (is_vect_null (ud->rvec,ud->ln) == 1)        /* check if ref vect is all zeroes */
00541                 {
00542                         return "Reference vector is all zeros"; 
00543                 }
00544                 
00545         ud->refname = tsim->name;
00546         ud->ignore = PLUTO_get_number(plint) ;    /* get number item */
00547         
00548         str = PLUTO_get_string(plint) ;
00549    ud->Dsamp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00550    
00551    /*--------- go to 3rd input line, sampling frequency, and stimulus period ---------*/
00552         
00553    PLUTO_next_option(plint) ;
00554    
00555    ud->fs = PLUTO_get_number(plint) ;    /* get number item */
00556    ud->T = PLUTO_get_number(plint) ;    /* get number item */
00557    
00558    ud->co = PLUTO_get_number(plint) ;    /* get number item */
00559    str = PLUTO_get_string(plint) ;
00560    ud->biasrem = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00561    
00562    /*--------- go to 4th input line, delay units and wrp option---------*/
00563                 
00564    PLUTO_next_option(plint) ;
00565 
00566    ud->Nseg = (int)PLUTO_get_number(plint) ;    /* get number item */
00567    ud->Pover = (int)PLUTO_get_number(plint) ;    /* get number item */
00568    
00569    str = PLUTO_get_string(plint) ;                                                      /* get string item (the method) */
00570    ud->unt = (int)PLUTO_string_index( str ,                                     /* find it in list it is from */
00571                                  NUM_METHOD_STRINGS ,
00572                                  method_strings ) ;
00573         
00574         str = PLUTO_get_string(plint) ;  
00575         ud->wrp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00576         
00577    /*--------- go to 5th input line Output prefix ---------*/
00578                 
00579    PLUTO_next_option(plint) ;
00580                 
00581    ud->new_prefix = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
00582         
00583         /* check to see if the field is empty */
00584         if (ud->new_prefix == NULL)
00585                         nprfx = 0;
00586                 else
00587                         nprfx = 1;
00588         /* check if the size is larger than 0. I did not want to check for this unless it's allocated */
00589         if (nprfx == 1 && (int)strlen (ud->new_prefix) == 0)
00590                 nprfx = 0;
00591                 
00592         if (nprfx == 0)         /* now create the new name and make new_prefix point to it */
00593                 {
00594                         sprintf (nprfxstr,"%s.DEL",DSET_PREFIX (old_dset));
00595                         ud->new_prefix = nprfxstr;
00596                         /*printf ("New prefix is set to be : %s\n\a",ud->new_prefix);*/
00597                 }
00598         
00599    if( ! PLUTO_prefix_ok(ud->new_prefix) )      /* check if it is OK */
00600       return "************************\n"
00601              "Output Prefix is illegal\n"
00602              "************************"  ;
00603         
00604         str = PLUTO_get_string(plint) ;                                 /* write delays to file ? */
00605         ud->out = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
00606         
00607         ud->strout = PLUTO_get_string(plint) ;                          /* strout is for the outiflename, which will be used after the debugging section */
00608         if (ud->strout == NULL)                                         /* if no output name is given, use the new_prefix */
00609                 {ud->strout = ud->new_prefix;}
00610                 else 
00611                         {       
00612                                 if((int)strlen (ud->strout) == 0) ud->strout = ud->new_prefix;
00613                         }
00614                         
00615         str = PLUTO_get_string(plint) ; 
00616         ud->outts = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
00617         
00618         /* ------------------Done with user parameters ---------------------------- */
00619         
00620         ud->nxx = (int)old_dset->daxes->nxx;                            /* get data set dimensions */
00621         ud->nyy = (int)old_dset->daxes->nyy;
00622         ud->nzz = (int)old_dset->daxes->nzz;
00623         
00624         /* No need for users to set these options ...*/
00625         
00626         ud->dtrnd = 0;
00627         
00628         if (ud->ln != (ud->nsamp - ud->ignore))
00629                 {
00630                         ud->errcode = ERROR_BADLENGTH;
00631                         return "***************************\n"
00632                                          "Bad time series length \n"
00633                                          "Check reference time series\n"
00634                                          " or the ignore parameter   \n"
00635                                          "***************************\n";
00636                 }
00637         
00638         if ((ud->unt < 0) || (ud->unt > 2))                                                                             /* unt error Check */
00639                 {
00640          ud->errcode = ERROR_WRONGUNIT;
00641          return "***********************\n"
00642                          " internal error: (ziad)\n"
00643                                          "unt values out of bound\n"
00644                                          "***********************\n";                   /*unt must be between 0 and 2 */
00645                 }
00646           
00647           if ((ud->wrp < 0) || (ud->wrp > 1))                                                                           /* wrp error Check */
00648                 {
00649          ud->errcode = ERROR_WARPVALUES;
00650          return "***********************\n"
00651                          " internal error: (ziad)\n"
00652                                          "wrp values out of bound\n"
00653                                          "***********************\n";                   /* wrp must be between 0 and 1*/
00654                 }
00655           
00656           if (ud->fs < 0.0) {                                         /* fs error Check */
00657          ud->errcode = ERROR_FSVALUES;
00658          return "***********************\n"
00659                          " internal error: (ziad)\n"
00660                                          "fs value is negative !\n"
00661                                          "***********************\n";                   /* fs must be >= 0*/
00662         }
00663           
00664           if (ud->T < 0.0) {                                          /* T error Check */
00665          ud->errcode = ERROR_TVALUES;
00666          return "***********************\n"
00667                          " internal error: (ziad)\n"
00668                                          "T value is negative !\n"
00669                                          "***********************\n";                                   /*T must be >= 0  */
00670         }
00671         
00672                 
00673      if ((ud->T == 0.0) && (ud->unt > 0))                                 /* unt error Check */
00674         {
00675          ud->errcode = ERROR_TaUNITVALUES;
00676          return "***********************\n"
00677                          " internal error: (ziad)\n"
00678                                          "T and unt val. mismatch\n"
00679                                          "***********************\n";                   /*T must be specified, and > 0 in order to use polar units*/
00680         }
00681 
00682     
00683     if ((ud->wrp == 1) && (ud->T == 0.0))                                   /* wrp error Check */
00684         {
00685          ud->errcode = ERROR_TaWRAPVALUES;
00686          return "***********************\n"
00687                          " internal error: (ziad)\n"
00688                                          "wrp and T val. mismatch\n"
00689                                          "***********************\n";                   /*T must be specified, and > 0 in order to use polar warp*/
00690         }
00691          if ((ud->out == NOPE) && (ud->outts == YUP))
00692                         {
00693                          ud->errcode = ERROR_OUTCONFLICT;
00694                          return"***********************\n"
00695                          "error: \n"
00696                                          "Write flag must be on\n"
00697                                          "to use Write ts\n"
00698                                          "***********************\n";   
00699                         
00700                         }
00701         
00702 
00703         /* Open the logfile, regardless of the ascii output files */
00704         sprintf ( tmpstr , "%s.log" , ud->strout);
00705         ud->outlogfile = fopen (tmpstr,"w");
00706 
00707 
00708         if (ud->out == YUP)                                                                     /* open outfile */
00709                                 {                                       
00710                                         ud->outwrite = fopen (ud->strout,"w");
00711                                         
00712                                         if (ud->outts == YUP)
00713                                                 {
00714                                                         sprintf ( tmpstr , "%s.ts" , ud->strout);
00715                                                         ud->outwritets = fopen (tmpstr,"w");
00716                                                         
00717                                                 }
00718                                         
00719                                         if ((ud->outwrite == NULL) || (ud->outlogfile == NULL) ||\
00720                                             (ud->outwritets == NULL && ud->outts == YUP) )
00721                                                 {
00722                                                         ud->errcode = ERROR_FILEOPEN; 
00723                                                         
00724                                                         return "***********************\n"
00725                                                                          "Could Not Write Outfile\n"
00726                                                                          "***********************\n";
00727                                                 }
00728         
00729                                 }
00730         
00731         /* Write out user variables to Logfile */
00732         write_ud (ud);                  /* writes user data to a file */
00733         
00734         /*show_ud (ud,0);       */                      /* For some debugging */
00735 
00736    
00737    /*------------- ready to compute new dataset -----------*/
00738 
00739    new_dset = MAKER_4D_to_typed_fbuc ( old_dset ,             /* input dataset */
00740                                ud->new_prefix ,           /* output prefix */
00741                                -1,                                                      /* negative value indicating data type is like original brick */
00742                                ud->ignore ,               /* ignore count */
00743                                1 ,                    /* detrend = ON Let BOB do it*/
00744                                NBUCKETS,                                        /*Number of values at each voxel*/
00745                                                                                  DELAY_tsfuncV2 ,         /* timeseries processor (bucket version)*/
00746                                                                                  (void *)ud          /* data for tsfunc */
00747                                                                                 ) ; 
00748                                                                                  
00749    /* Setup the label, keywords and types of subbricks */
00750         i = 0;
00751         while (i < NBUCKETS)
00752                 {
00753                         switch (i)
00754                                 {
00755                                         case DELINDX:                                   /* delay value in results vector */
00756                                                 EDIT_BRICK_LABEL (new_dset,i,"Delay");
00757                                                 EDIT_BRICK_ADDKEY (new_dset,i,"D");
00758                                                 ++i;
00759                                                 break;
00760                                         case COVINDX:                                   /* covariance value in results vector */
00761                                                 EDIT_BRICK_LABEL (new_dset,i,"Covariance");
00762                                                 EDIT_BRICK_ADDKEY (new_dset,i,"I");
00763                                                 ++i;
00764                                                 break;
00765                                         case COFINDX:                                   /* cross correlation coefficient value in results vector */
00766                                                 EDIT_BRICK_LABEL (new_dset,i,"Corr. Coef.");
00767                                                 EDIT_BRICK_ADDKEY (new_dset,i,"r");
00768                                                 /* Here you must modify either ud->Nfit or ud->Nort or most likely ud->nsamp based on ud->Navg */
00769                                                 EDIT_BRICK_TO_FICO (new_dset,i,ud->nsamp - ud->ignore,ud->Nfit,ud->Nort);
00770                                                 ++i;
00771                                                 break;
00772                                         case VARINDX:                                   /* FMRI time course variance value in results vector */
00773                                                 EDIT_BRICK_LABEL (new_dset,i,"Variance");
00774                                                 EDIT_BRICK_ADDKEY (new_dset,i,"S2");
00775                                                 ++i;
00776                                                 break;
00777                                         default :
00778                                                 return "*********************\n"
00779                                                                  "Internal Error (ziad)\n"
00780                                                                  " Bad i value \n"
00781                                                                  "*********************\n";
00782                                                 break;
00783                                 }
00784                                 
00785                 }
00786         
00787         PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
00788         
00789         
00790 
00791    if (ud->out == YUP)                                                                  /* close outfile and outlogfile*/
00792                                 {
00793                                         fclose (ud->outlogfile);
00794                                         fclose (ud->outwrite);
00795                                         if (ud->outts  == YUP) fclose (ud->outwritets);
00796                                 }
00797                                 else
00798                                 {
00799                                         if (ud->outlogfile != NULL)     fclose (ud->outlogfile);                /* close outlogfile */
00800                                 }
00801         
00802         free (tmpstr);          
00803         free (nprfxstr);
00804    return NULL ;  /* null string returned means all was OK */
00805 }
00806 
00807 /**********************************************************************
00808    Function that does the real work
00809 ***********************************************************************/
00810 
00811 static void DELAY_tsfuncV2( double T0 , double TR ,
00812                    int npts , float ts[] , double ts_mean , double ts_slope ,
00813                    void * udp , int nbrick , float * buckar)
00814 {
00815    static int nvox , ncall ;
00816         hilbert_data_V2 uda,*ud;
00817         float del, xcorCoef, buckara[4];
00818         float xcor=0.0 ,  tmp=0.0 , tmp2 = 0.0 ,  dtx = 0.0 ,\
00819                          delu = 0.0 , slp = 0.0 , vts = 0.0 , vrvec = 0.0 ;
00820         int i , is_ts_null , status , opt , actv , zpos , ypos , xpos ;
00821         
00822         ud = &uda;
00823         ud = (hilbert_data_V2 *) udp;
00824         
00825    /** is this a "notification"? **/
00826 
00827    if( buckar == NULL ){
00828 
00829       if( npts > 0 ){  /* the "start notification" */
00830 
00831          PLUTO_popup_meter( global_plint ) ;  /* progress meter  */
00832          nvox  = npts ;                       /* keep track of   */
00833          ncall = 0 ;                          /* number of calls */
00834                         
00835       } else {  /* the "end notification" */
00836                         
00837                         opt = 0;                                        /* cleanup in hdelay */
00838                 status = hilbertdelay_V2 (ts,ud->rvec,ud->ln,ud->Nseg,ud->Pover,opt,ud->dtrnd,dtx,ud->biasrem,&tmp,&slp,&xcor,&tmp2,&vts,&vrvec);                                       /* cleanup time */
00839         
00840          PLUTO_set_meter( global_plint , 100 ) ; /* set meter to 100% */
00841 
00842       }
00843       return ;
00844    }
00845 
00846         /* In the old version, I had to check for a match between the lengths of the reference time series and FMRI time series
00847         This is now done before the function is called. */
00848    
00849    if (is_vect_null (ts,npts) == 1) /* check for null vectors */
00850         {
00851                 ud->errcode = ERROR_NULLTIMESERIES;
00852                         error_report (ud , ncall );     /* report the error */
00853                 
00854                 del = 0.0;                                                              /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
00855                 xcorCoef = 0.0;                                         /*  because the data might still be OK */
00856                 xcor = 0.0;
00857         }
00858    
00859    if (ud->errcode == 0)                                /* if there are no errors, proceed */   
00860                 {/* ud->errcode == 0 outer loop */
00861                         opt = 1;                                        /* activate hdelay */
00862                 
00863                 /* transform dtx from seconds to sampling units and correct for the number of points ignored*/
00864                 if (ud->Dsamp == YUP) 
00865                         dtx = (float) (T0 / TR) - ud->ignore;
00866                 else
00867                         dtx = 0.0;
00868                         
00869                 ud->errcode = hilbertdelay_V2 (ts,ud->rvec,ud->ln,ud->Nseg,ud->Pover,opt,ud->dtrnd,dtx,ud->biasrem,&delu,&slp,&xcor,&xcorCoef,&vts,&vrvec);                                     /* cleanup time */
00870                         
00871                         if (ud->errcode == 0) /* If there are no errors, proceed */
00872                                 { /*ud->errcode == 0 inner loop */
00873                                         hunwrap (delu, (float)(1/TR), ud->T, slp, ud->wrp, ud->unt, &del );
00874                                         
00875                                         actv = 1;                                               /* assume voxel is active */
00876         
00877                                         if (xcorCoef < ud->co) actv = 0;                        /* determine if voxel is activated using xcorCoef  */
00878         
00879                                         if ((actv == 1) && (ud->out == YUP))            /* if voxel is truly activated, write results to file without modifying return value */
00880                                                 {
00881                                                         indexTOxyz ( ud , ncall, &xpos , &ypos , &zpos);
00882                                                         fprintf (ud->outwrite,"%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\n", ncall , xpos , ypos , zpos ,  delu , del , xcor , xcorCoef , vts);
00883                                                         if (ud->outts == YUP)
00884                                                                 {
00885                                                                         writets (ud,ts);        
00886                                                                 }
00887                                                 }
00888 
00889                                 }/*ud->errcode == 0 inner loop */
00890                         
00891                         else if (ud->errcode == ERROR_LONGDELAY)
00892                                                 {
00893                                                         error_report ( ud , ncall);     
00894                                                         
00895                                                         del = 0.0;                                                              /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
00896                                                 xcorCoef = 0.0;                                         /*  because the data might still be OK */
00897                                                 xcor = 0.0;
00898                                                         
00899                                                 }
00900                                         else if (ud->errcode != 0)
00901                                                                 {
00902                                                                         error_report ( ud , ncall);     
00903                                                                         
00904                                                                         del = 0.0;                                                              /* Set all the variables to Null and set xcorCoef to an impossible value*/
00905                                                                 xcorCoef = NOWAYXCORCOEF;                                               
00906                                                                 xcor = 0.0;
00907                                                                 }
00908 
00909    }/* ud->errcode == 0 outer loop */
00910    
00911         /* Now fill up the bucket array */
00912         
00913         buckar[DELINDX] = del;
00914         buckar[COVINDX] = xcor;
00915         buckar[COFINDX] = xcorCoef;
00916         buckar[VARINDX] = vts;
00917 
00918 
00919    /** set the progress meter to the % of completion **/
00920    ncall++ ;
00921    
00922    PLUTO_set_meter( global_plint , (100*ncall)/nvox ) ;
00923    
00924    ud->errcode = 0;     /* Rest error to nothing */
00925    
00926    return ;
00927 }
00928 
00929 /* ************************************************************ */
00930 /* function to display user data input (debugging stuff)        */
00931 /* ************************************************************ */
00932 
00933 static void show_ud (hilbert_data_V2* ud,int loc)
00934         {
00935                 printf ("\n\nUser Data Values at location :%d\n",loc);
00936                 printf ("ud->dsetname= %s\n",ud->dsetname);
00937                 printf ("ud->refname= %s\n",ud->refname);
00938                 printf ("ud->rvec= (1st three elements only)");
00939                 disp_vect (ud->rvec,3);
00940                 printf ("ud->nxx= %d\n",ud->nxx);
00941                 printf ("ud->nyy= %d\n",ud->nyy);
00942                 printf ("ud->nzz= %d\n",ud->nzz);
00943                 printf ("ud->fs= %f\n",ud->fs);
00944                 printf ("ud->T= %f\n",ud->T);
00945                 printf ("ud->co= %f\n",ud->co);
00946                 printf ("ud->unt= %d\n",ud->unt);
00947                 printf ("ud->wrp= %d\n",ud->wrp);
00948                 printf ("ud->Navg= %d\n",ud->Navg);
00949                 printf ("ud->Nort= %d\n",ud->Nort);
00950                 printf ("ud->Nfit= %d\n",ud->Nfit);
00951                 printf ("ud->Nseg= %d\n",ud->Nseg);
00952                 printf ("ud->Pover= %d\n",ud->Pover);
00953                 printf ("ud->dtrnd= %d\n",ud->dtrnd);
00954                 printf ("ud->biasrem= %d\n",ud->biasrem);
00955                 printf ("ud->Dsamp= %d\n",ud->Dsamp);
00956                 printf ("ud->ln= %d\n",ud->ln);
00957                 printf ("ud->nsamp= %d\n",ud->nsamp);
00958                 printf ("ud->ignore= %d\n",ud->ignore);
00959                 printf ("ud->errcode= %d\n",ud->errcode);
00960                 printf ("ud->new_prefix= %s\n",ud->new_prefix);
00961                 printf ("ud->out= %d\n",ud->out);
00962                 printf ("ud->strout= %s\n",ud->strout);
00963                 printf ("ud->outts= %d\n",ud->outts);
00964                 printf ("Hit enter to continue\a\n\n");
00965                 getchar ();
00966                 return;
00967         }
00968 
00969 /* ************************************************************ */
00970 /* function to write user data input to log file        */
00971 /* ************************************************************ */
00972 
00973 static void write_ud (hilbert_data_V2* ud)
00974         {
00975                 fprintf (ud->outlogfile,"\nLogfile output by Hilbert Delay98 plugin\n");
00976                 fprintf (ud->outlogfile,"\n\nUser Data Values \n");
00977                 fprintf (ud->outlogfile,"Input data set = %s\n",ud->dsetname);
00978                 fprintf (ud->outlogfile,"Reference file name = %s\n",ud->refname);
00979                 fprintf (ud->outlogfile,"Number of voxels in X dimension = %d\n",ud->nxx);
00980                 fprintf (ud->outlogfile,"Number of voxels in Y dimension = %d\n",ud->nyy);
00981                 fprintf (ud->outlogfile,"Number of voxels in Z dimension = %d\n",ud->nzz);
00982                 fprintf (ud->outlogfile,"Sampling Frequency = %f\n",ud->fs);
00983                 fprintf (ud->outlogfile,"Stimulus Period = %f\n",ud->T);
00984                 fprintf (ud->outlogfile,"Threshold Cut Off value = %f\n",ud->co);
00985                 fprintf (ud->outlogfile,"Delay units = %d\n",ud->unt);
00986                 fprintf (ud->outlogfile,"Delay wrap = %d\n",ud->wrp);
00987                 fprintf (ud->outlogfile,"Number of segments = %d\n",ud->Nseg);
00988                 fprintf (ud->outlogfile,"Number of samples in time series = %d\n",ud->nsamp);
00989                 fprintf (ud->outlogfile,"Ignore = %d\n",ud->ignore);
00990                 fprintf (ud->outlogfile,"Length of reference time series = %d\n",ud->ln);
00991                 fprintf (ud->outlogfile,"Number of fit parameters = %d\n",ud->Nfit);
00992                 fprintf (ud->outlogfile,"Number of nuisance parameters (orts)= %d\n",ud->Nort);
00993                 fprintf (ud->outlogfile,"Percent overlap = %d\n",ud->Pover);
00994                 fprintf (ud->outlogfile,"Plugin detrending = %d (Always 0, mandatory detrending is performed)\n",ud->dtrnd);
00995                 fprintf (ud->outlogfile,"Bias correction = %d\n",ud->biasrem);
00996                 fprintf (ud->outlogfile,"Acquisition time correction = %d\n",ud->Dsamp);
00997                 fprintf (ud->outlogfile,"Prefix for birck output = %s\n",ud->new_prefix);
00998                 fprintf (ud->outlogfile,"Flag for Ascii output file  = %d\n",ud->out);
00999                 fprintf (ud->outlogfile,"Ascii output file name = %s\n",ud->strout);
01000                 fprintf (ud->outlogfile,"Flag for Ascii time series file output = %d\n",ud->outts);
01001                 fprintf (ud->outlogfile,"\nud->errcode (debugging only)= %d\n\n",ud->errcode);
01002                 fprintf (ud->outlogfile,"\nThe format for the output file is the following:\n");
01003                 fprintf (ud->outlogfile,"VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
01004                 fprintf (ud->outlogfile,"\nError Log <message> <index> <x> <y> <z>\n\n");
01005                 
01006                 return;
01007         }
01008 
01009 /* ************************************************************ */ 
01010 /* function to compute x, y, z coordinates from the index       */
01011 /* ************************************************************ */ 
01012 
01013 static void indexTOxyz (hilbert_data_V2* ud, int ncall, int *xpos , int *ypos , int *zpos)      
01014         {
01015                 *zpos = (int)ncall / (int)(ud->nxx*ud->nyy);
01016                 *ypos = (int)(ncall - *zpos * ud->nxx * ud->nyy) / (int)ud->nxx;
01017                 *xpos = ncall - ( *ypos * ud->nxx ) - ( *zpos * ud->nxx * ud->nyy ) ;
01018                 return;
01019         }
01020         
01021 /* ************************************************************ */
01022 /* function to report errors encountered to the logfile         */
01023 /* Only errors that happen during runtime (while delays are      */
01024 /* computed, are handeled here, the others are handeled                  */
01025 /* instantaneously, and need not be logged                                                       */
01026 /* ************************************************************ */
01027 
01028 static void error_report (hilbert_data_V2* ud, int ncall ) 
01029         {
01030                 int xpos,ypos,zpos;
01031                 indexTOxyz (ud, ncall, &xpos , &ypos , &zpos); 
01032 
01033                 switch (ud->errcode)
01034                         {
01035                                 case ERROR_NOTHINGTODO:
01036                                         fprintf (ud->outlogfile,"Nothing to do hilbertdelay_V2 call ");
01037                                         break;
01038                                 case ERROR_LARGENSEG:
01039                                         fprintf (ud->outlogfile,"Number of segments Too Large ");
01040                                         break;
01041                                 case ERROR_LONGDELAY:
01042                                         fprintf (ud->outlogfile,"Could not find zero crossing before maxdel limit ");
01043                                         break;
01044                                 case ERROR_SERIESLENGTH:
01045                                         fprintf (ud->outlogfile,"Vectors have different length ");
01046                                         break;
01047                                 case ERROR_NULLTIMESERIES:
01048                                         fprintf (ud->outlogfile,"Null time series vector ");
01049                                         break;
01050                                 default:
01051                                         fprintf (ud->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",ud->errcode);
01052                                         break;
01053                         }       
01054                 fprintf (ud->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos  );
01055                 return;
01056         }
01057         
01058 /* *************************************************************** */
01059 /* function to write the time course into a line in the given file */
01060 /* *************************************************************** */
01061 
01062 static void writets (hilbert_data_V2 * ud,float * ts)
01063 
01064         {       
01065                 int i;
01066                 
01067                 for (i=0;i<ud->ln;++i)
01068                         {
01069                                 fprintf (ud->outwritets, "%f\t",ts[i]);
01070                         }
01071                 fprintf (ud->outwritets,"\n");
01072         }