00001
00002
00003
00004
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
00015
00016
00017 #include <stdlib.h>
00018 #include <stdio.h>
00019 #include <math.h>
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 typedef struct {
00034 float real, imag;
00035 } COMPLEX;
00036
00037
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
00055
00056
00057 typedef struct
00058 {
00059 int nxx;
00060 int nyy;
00061 int nzz;
00062 char *dsetname;
00063 char *refname;
00064 float *rvec;
00065 float fs;
00066
00067
00068 float T;
00069 float co;
00070 float dmask;
00071 int unt;
00072 int wrp;
00073 int Navg;
00074 int Nort;
00075 int Nfit;
00076 int Nseg;
00077 int nsamp;
00078 int ignore;
00079 int Pover;
00080 int ln;
00081 int dtrnd;
00082 int biasrem;
00083 int Dsamp;
00084 int errcode;
00085 int out;
00086 int outts;
00087 char * new_prefix;
00088 char * strout;
00089 FILE * outwrite;
00090 FILE * outwritets;
00091 FILE * outlogfile;
00092 char outname[PLUGIN_MAX_STRING_RANGE];
00093 }hilbert_data_V2;
00094
00095
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
00214
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
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
00232 #endif
00233
00234 #define NBUCKETS 4
00235 #define DELINDX 0
00236 #define COVINDX 1
00237 #define COFINDX 2
00238 #define VARINDX 3
00239
00240 #define YUP 1
00241 #define NOPE 0
00242
00243 #define ERROR_NOTHINGTODO 1
00244 #define ERROR_LARGENSEG 2
00245 #define ERROR_LONGDELAY 3
00246 #define ERROR_WRONGUNIT 8
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
00260
00261 static char * DELAY_main( PLUGIN_interface * ) ;
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
00278
00279 static PLUGIN_interface * global_plint = NULL ;
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 DEFINE_PLUGIN_PROTOTYPE
00296
00297 PLUGIN_interface * PLUGIN_init( int ncall )
00298 {
00299 PLUGIN_interface * plint ;
00300
00301 if( ncall > 0 ) return NULL ;
00302
00303
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 ;
00311
00312
00313
00314 PLUTO_add_option( plint ,
00315 "Data" ,
00316 "Data" ,
00317 TRUE
00318 ) ;
00319
00320 PLUTO_add_dataset( plint ,
00321 "3D+time" ,
00322 ANAT_ALL_MASK ,
00323 FUNC_ALL_MASK ,
00324 DIMEN_4D_MASK |
00325 BRICK_ALLREAL_MASK
00326 ) ;
00327
00328 PLUTO_add_number( plint ,
00329 "Nort" ,
00330 1 ,
00331 100 ,
00332 0 ,
00333 2 ,
00334 FALSE
00335 ) ;
00336
00337
00338
00339 PLUTO_add_option( plint ,
00340 "Ref." ,
00341 "Ref." ,
00342 TRUE
00343 ) ;
00344
00345 PLUTO_add_timeseries(plint,"Ref. Vect.");
00346
00347 PLUTO_add_number( plint ,
00348 "Ignore" ,
00349 0 ,
00350 50 ,
00351 0 ,
00352 0 ,
00353 FALSE
00354 ) ;
00355
00356 PLUTO_add_string( plint ,
00357 "Dsamp" ,
00358 2,yn_strings,
00359 1
00360 ) ;
00361
00362
00363
00364 PLUTO_add_option( plint ,
00365 "Sig." ,
00366 "Sig." ,
00367 TRUE
00368 ) ;
00369
00370 PLUTO_add_number( plint ,
00371 "fs in Hz" ,
00372 0 ,
00373 2000 ,
00374 1 ,
00375 5 ,
00376 TRUE
00377 ) ;
00378
00379 PLUTO_add_number( plint ,
00380 "Tstim sec" ,
00381 0.0 ,
00382 500 ,
00383 0 ,
00384 40 ,
00385 TRUE
00386 ) ;
00387
00388 PLUTO_add_number( plint ,
00389 "C-Off" ,
00390 -10 ,
00391 10 ,
00392 1 ,
00393 5 ,
00394 TRUE
00395 ) ;
00396
00397
00398 PLUTO_add_string( plint ,
00399 "No-bias" ,
00400 2,yn_strings,
00401 1
00402 ) ;
00403
00404
00405
00406
00407
00408 PLUTO_add_option( plint ,
00409 "Alg." ,
00410 "Alg." ,
00411 TRUE
00412 ) ;
00413
00414 PLUTO_add_number( plint ,
00415 "N seg." ,
00416 1 ,
00417 1 ,
00418 0 ,
00419 1 ,
00420 FALSE
00421 ) ;
00422
00423 PLUTO_add_number( plint ,
00424 "% ovrlp" ,
00425 0 ,
00426 0 ,
00427 0 ,
00428 0 ,
00429 FALSE
00430 ) ;
00431
00432
00433 PLUTO_add_string( plint ,
00434 "Units" ,
00435 3,method_strings,
00436 0
00437 ) ;
00438
00439 PLUTO_add_string( plint ,
00440 "Phz Wrp" ,
00441 2,yn_strings,
00442 0
00443 ) ;
00444
00445
00446
00447
00448 PLUTO_add_option( plint ,
00449 "Output" ,
00450 "Output" ,
00451 TRUE
00452 ) ;
00453
00454 PLUTO_add_string( plint ,
00455 "AFNI Prfx" ,
00456 0,NULL ,
00457 19
00458 ) ;
00459
00460 PLUTO_add_string( plint ,
00461 "Write" ,
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" ,
00470 2,yn_strings ,
00471 1
00472 ) ;
00473
00474
00475 return plint ;
00476 }
00477
00478
00479
00480
00481
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 ;
00489 THD_3dim_dataset * old_dset , * new_dset ;
00490 char *tmpstr , * str , *nprfxstr;
00491 int ntime, nvec ,nprfx, i;
00492 float * vec , fs , T ;
00493
00494
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;
00507 ud->errcode = 0;
00508
00509
00510
00511
00512
00513
00514 PLUTO_next_option(plint) ;
00515
00516 idc = PLUTO_get_idcode(plint) ;
00517 old_dset = PLUTO_find_dset(idc) ;
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 ;
00526 ud->Nort = PLUTO_get_number(plint) ;
00527 ud->Nfit = 2 ;
00528
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;
00536 nvec = tsim -> ny;
00537 ud->rvec = (float *) MRI_FLOAT_PTR(tsim);
00538
00539
00540 if (is_vect_null (ud->rvec,ud->ln) == 1)
00541 {
00542 return "Reference vector is all zeros";
00543 }
00544
00545 ud->refname = tsim->name;
00546 ud->ignore = PLUTO_get_number(plint) ;
00547
00548 str = PLUTO_get_string(plint) ;
00549 ud->Dsamp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00550
00551
00552
00553 PLUTO_next_option(plint) ;
00554
00555 ud->fs = PLUTO_get_number(plint) ;
00556 ud->T = PLUTO_get_number(plint) ;
00557
00558 ud->co = PLUTO_get_number(plint) ;
00559 str = PLUTO_get_string(plint) ;
00560 ud->biasrem = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
00561
00562
00563
00564 PLUTO_next_option(plint) ;
00565
00566 ud->Nseg = (int)PLUTO_get_number(plint) ;
00567 ud->Pover = (int)PLUTO_get_number(plint) ;
00568
00569 str = PLUTO_get_string(plint) ;
00570 ud->unt = (int)PLUTO_string_index( str ,
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
00578
00579 PLUTO_next_option(plint) ;
00580
00581 ud->new_prefix = PLUTO_get_string(plint) ;
00582
00583
00584 if (ud->new_prefix == NULL)
00585 nprfx = 0;
00586 else
00587 nprfx = 1;
00588
00589 if (nprfx == 1 && (int)strlen (ud->new_prefix) == 0)
00590 nprfx = 0;
00591
00592 if (nprfx == 0)
00593 {
00594 sprintf (nprfxstr,"%s.DEL",DSET_PREFIX (old_dset));
00595 ud->new_prefix = nprfxstr;
00596
00597 }
00598
00599 if( ! PLUTO_prefix_ok(ud->new_prefix) )
00600 return "************************\n"
00601 "Output Prefix is illegal\n"
00602 "************************" ;
00603
00604 str = PLUTO_get_string(plint) ;
00605 ud->out = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
00606
00607 ud->strout = PLUTO_get_string(plint) ;
00608 if (ud->strout == NULL)
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
00619
00620 ud->nxx = (int)old_dset->daxes->nxx;
00621 ud->nyy = (int)old_dset->daxes->nyy;
00622 ud->nzz = (int)old_dset->daxes->nzz;
00623
00624
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))
00639 {
00640 ud->errcode = ERROR_WRONGUNIT;
00641 return "***********************\n"
00642 " internal error: (ziad)\n"
00643 "unt values out of bound\n"
00644 "***********************\n";
00645 }
00646
00647 if ((ud->wrp < 0) || (ud->wrp > 1))
00648 {
00649 ud->errcode = ERROR_WARPVALUES;
00650 return "***********************\n"
00651 " internal error: (ziad)\n"
00652 "wrp values out of bound\n"
00653 "***********************\n";
00654 }
00655
00656 if (ud->fs < 0.0) {
00657 ud->errcode = ERROR_FSVALUES;
00658 return "***********************\n"
00659 " internal error: (ziad)\n"
00660 "fs value is negative !\n"
00661 "***********************\n";
00662 }
00663
00664 if (ud->T < 0.0) {
00665 ud->errcode = ERROR_TVALUES;
00666 return "***********************\n"
00667 " internal error: (ziad)\n"
00668 "T value is negative !\n"
00669 "***********************\n";
00670 }
00671
00672
00673 if ((ud->T == 0.0) && (ud->unt > 0))
00674 {
00675 ud->errcode = ERROR_TaUNITVALUES;
00676 return "***********************\n"
00677 " internal error: (ziad)\n"
00678 "T and unt val. mismatch\n"
00679 "***********************\n";
00680 }
00681
00682
00683 if ((ud->wrp == 1) && (ud->T == 0.0))
00684 {
00685 ud->errcode = ERROR_TaWRAPVALUES;
00686 return "***********************\n"
00687 " internal error: (ziad)\n"
00688 "wrp and T val. mismatch\n"
00689 "***********************\n";
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
00704 sprintf ( tmpstr , "%s.log" , ud->strout);
00705 ud->outlogfile = fopen (tmpstr,"w");
00706
00707
00708 if (ud->out == YUP)
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
00732 write_ud (ud);
00733
00734
00735
00736
00737
00738
00739 new_dset = MAKER_4D_to_typed_fbuc ( old_dset ,
00740 ud->new_prefix ,
00741 -1,
00742 ud->ignore ,
00743 1 ,
00744 NBUCKETS,
00745 DELAY_tsfuncV2 ,
00746 (void *)ud
00747 ) ;
00748
00749
00750 i = 0;
00751 while (i < NBUCKETS)
00752 {
00753 switch (i)
00754 {
00755 case DELINDX:
00756 EDIT_BRICK_LABEL (new_dset,i,"Delay");
00757 EDIT_BRICK_ADDKEY (new_dset,i,"D");
00758 ++i;
00759 break;
00760 case COVINDX:
00761 EDIT_BRICK_LABEL (new_dset,i,"Covariance");
00762 EDIT_BRICK_ADDKEY (new_dset,i,"I");
00763 ++i;
00764 break;
00765 case COFINDX:
00766 EDIT_BRICK_LABEL (new_dset,i,"Corr. Coef.");
00767 EDIT_BRICK_ADDKEY (new_dset,i,"r");
00768
00769 EDIT_BRICK_TO_FICO (new_dset,i,ud->nsamp - ud->ignore,ud->Nfit,ud->Nort);
00770 ++i;
00771 break;
00772 case VARINDX:
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)
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);
00800 }
00801
00802 free (tmpstr);
00803 free (nprfxstr);
00804 return NULL ;
00805 }
00806
00807
00808
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
00826
00827 if( buckar == NULL ){
00828
00829 if( npts > 0 ){
00830
00831 PLUTO_popup_meter( global_plint ) ;
00832 nvox = npts ;
00833 ncall = 0 ;
00834
00835 } else {
00836
00837 opt = 0;
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);
00839
00840 PLUTO_set_meter( global_plint , 100 ) ;
00841
00842 }
00843 return ;
00844 }
00845
00846
00847
00848
00849 if (is_vect_null (ts,npts) == 1)
00850 {
00851 ud->errcode = ERROR_NULLTIMESERIES;
00852 error_report (ud , ncall );
00853
00854 del = 0.0;
00855 xcorCoef = 0.0;
00856 xcor = 0.0;
00857 }
00858
00859 if (ud->errcode == 0)
00860 {
00861 opt = 1;
00862
00863
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);
00870
00871 if (ud->errcode == 0)
00872 {
00873 hunwrap (delu, (float)(1/TR), ud->T, slp, ud->wrp, ud->unt, &del );
00874
00875 actv = 1;
00876
00877 if (xcorCoef < ud->co) actv = 0;
00878
00879 if ((actv == 1) && (ud->out == YUP))
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 }
00890
00891 else if (ud->errcode == ERROR_LONGDELAY)
00892 {
00893 error_report ( ud , ncall);
00894
00895 del = 0.0;
00896 xcorCoef = 0.0;
00897 xcor = 0.0;
00898
00899 }
00900 else if (ud->errcode != 0)
00901 {
00902 error_report ( ud , ncall);
00903
00904 del = 0.0;
00905 xcorCoef = NOWAYXCORCOEF;
00906 xcor = 0.0;
00907 }
00908
00909 }
00910
00911
00912
00913 buckar[DELINDX] = del;
00914 buckar[COVINDX] = xcor;
00915 buckar[COFINDX] = xcorCoef;
00916 buckar[VARINDX] = vts;
00917
00918
00919
00920 ncall++ ;
00921
00922 PLUTO_set_meter( global_plint , (100*ncall)/nvox ) ;
00923
00924 ud->errcode = 0;
00925
00926 return ;
00927 }
00928
00929
00930
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
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
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
01023
01024
01025
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
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 }