1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "afni.h"
8 #include "afni_plugin.h"
9
10 #ifndef ALLOW_PLUGINS
11 # error "Plugins not properly set up -- see machdep.h"
12 #endif
13
14 /* -----------------------START---------------------------------*/
15 /* definition and decleration part to suport the main algorithm */
16
17 #include <stdlib.h>
18 #include <stdio.h>
19 #include <math.h>
20
21 /*-------------------------------------------------------------------*/
22 /* taken from #include "/usr/people/ziad/Programs/C/DSP_in_C/dft.h" */
23 /* dft.h - function prototypes and structures for dft and fft functions */
24
25 #include "plug_delay_V2.h"
26
27
28 /***********************************************************************
29 Plugin to compute a 3D+time dataset voxelwise delay with respect to
30 a reference waveform
31 ************************************************************************/
32 typedef struct
33 {
34 int nxx; /* number of voxels in the x direction */
35 int nyy; /* number of voxels in the y direction */
36 int nzz; /* number of voxels in the z direction */
37 char *dsetname; /* prefix of data set analyzed */
38 char *refname; /* reference vector filename */
39 float *rvec; /* reference vetor */
40 float fs; /* Sampling frequency */
41 /* it is only used for the log file in this version*/
42 /* the ts_func, gives TR automatically */
43 float T; /* Stimulus period */
44 float co; /* Correlation Coefficient Threshold*/
45 float dmask; /* delays to be masked will take this value ( Not used anymore )*/
46 int unt; /* Delay units */
47 int wrp; /* flag for Polar Wrap */
48 int Navg; /* number of data sets averaged to obtain the brick (for statistical stuff) */
49 int Nort; /* Number of nuisance parameters (orts) (for statistical stuff) */
50 int Nfit; /* Number of fit parameters (for statistical stuff) */
51 int Nseg; /* Number of segments */
52 int nsamp; /* number of points in time series */
53 int ignore; /* number ofpoints to ignore from time courses */
54 int Pover; /* Percent overlap */
55 int ln; /* length of FMRI vector */
56 int dtrnd; /* remove linear trend or just the mean */
57 int biasrem; /* flag for removing delay bias */
58 int Dsamp; /* flag for correction of non uniform sampling start time */
59 int errcode; /* error code number returned from hdelay */
60 int out; /* flag for writing delays to a file */
61 int outts; /* flag for writing time series to a file */
62 char * new_prefix; /* new prefix for data set */
63 char * strout;
64 FILE * outwrite;
65 FILE * outwritets;
66 FILE * outlogfile;
67 char outname[PLUGIN_MAX_STRING_RANGE]; /* output data file name */
68 }hilbert_data_V2;
69
70 /*--------------------- string to 'help' the user --------------------*/
71
72 static char helpstring[] =
73 " Hilbert Delay98 Plugin \n\n"
74 "The Plugin estimates the time delay between each voxel time series in a 3D+time dataset \n"
75 "and a reference time series[1][2]. The estimated delays are relative to the reference time series.\n"
76 "For example, a delay of 4 seconds means that the voxel time series is delayed by 4 seconds with respect\n"
77 "to the reference time series.\n\n"
78 "Plugin Inputs:\n\n"
79 " 1- Data : \n"
80 " 3d+time -> 3D+time dataset to analyze.\n"
81 " Nort -> Number of orts or nuisance parameters. \n"
82 " It is set to 2 by default because the mean \n"
83 " and linear trends are always removed from the time series.\n"
84 " This value is used for estimating the p-value at the cross\n"
85 " correlation threshold selected with AFNI's threshold slider.\n"
86 " The p-value is estimated only when the cross correlation coefficient\n"
87 " subbrick is used for thresholding.\n\n"
88 " 2- Ref. :\n"
89 " Ref. Vect. -> 1D Reference time series vector. \n"
90 " The reference time series vector is stored in an ascii file.\n"
91 " The plugin assumes that there is one value per line and that all\n"
92 " values in the file are part of the reference vector.\n"
93 " PS: Unlike with 3dfim, and FIM in AFNI, values over 33333 are treated \n"
94 " as part of the time series.\n"
95 " The reference vectors must be placed in a directory specified in \n"
96 " AFNI_TSPATH environment variable. Read AFNI documentation for more info.\n"
97 " Ignore -> Number of samples to ignore from the beginning of each voxel's time series.\n"
98 " Ignore is NOT applied to the reference time series.\n"
99 " Dsamp -> (Y/N) Correct a voxel's estimated delay by the time at which the slice\n"
100 " containing that voxel was acquired.\n\n"
101 " 3- Sig. :\n"
102 " fs in Hz -> Sampling frequency in Hz. of data time series. \n"
103 " Tstim sec -> Stimulus period in seconds. \n"
104 " If the stimulus is not periodic, you can set Tstim to 0.\n"
105 " C-Off -> Cross Correlation Coefficient threshold value.\n"
106 " This is only used to censor the ascii output (see below).\n"
107 " No-bias -> (y/n) Correct for the bias in the cross correlation coefficient estimate [1][2].\n\n"
108 " 4- Alg. :\n"
109 " N seg. -> Number of segments used to estimate the periodogram.\n"
110 " (you can't modify this parameter in this version)\n"
111 " % ovrlp -> Percent segment overlap when estimating periodogram.\n"
112 " (you can't modify this parameter in this version)\n"
113 " Units -> (Seconds/Degrees/Radians) Units for delay estimates.\n"
114 " You can't use Degrees or Radians as units unless you specify\n"
115 " a value for Tstim > 0.\n"
116 " Phz Wrp -> (Y/N) Delay (or phase) wrap.\n"
117 " This switch maps delays from: \n"
118 " (Seconds) 0->T/2 to 0->T/2 and T/2->T to -T/2->0\n"
119 " (Degrees) 0->180 to 0->180 and 180->360 to -180->0\n"
120 " (Radians) 0->pi to 0->pi and pi->2pi to -pi->0\n"
121 " You can't use this option unless you specify a value for Tstim > 0.\n\n"
122 " 5- Output :\n"
123 " AFNI Prfx -> Prefix for output brick of the bucket type (fbuc).\n"
124 " The first subbrick is for Delay.\n"
125 " The second subbrick is for Covariance, which is an estimate\n"
126 " of the power in voxel time series at the frequencies present \n"
127 " in the reference time series.\n"
128 " The third subbrick is for the Cross Correlation Coefficients between\n"
129 " FMRI time series and reference time series.\n"
130 " The fourth subbrick contains estimates of the Variance of voxel time series.\n"
131 " If you leave the field empty, a default prefix is used.\n"
132 " The default prefix is the prefix of the input 3D+time brick \n"
133 " with a '.DEL' extension appended to it.\n"
134 " Write -> (Y/N) Write the results to an ascii file for voxels with \n"
135 " Cross Correlation Coefficients larger than C-Off.\n"
136 " Each line in the output file contains information for one voxel.\n"
137 " There a 9 columns in the file which hold the following values:\n"
138 " 1- Voxel Index (VI) : Each voxel in an AFNI brick has a unique index.\n"
139 " Indices map directly to XYZ coordinates.\n"
140 " See AFNI plugin documentations for more info.\n"
141 " 2..4- Voxel coordinates (X Y Z): Those are the voxel slice coordinates.\n"
142 " You can see these coordinates in the upper left side\n"
143 " of the AFNI window. To do so, you must first switch the\n"
144 " voxel coordinate units from mm to slice coordinates.\n"
145 " Define Datamode -> Misc -> Voxel Coords ?\n"
146 " PS: The coords that show up in the graph window\n"
147 " could be different from those in the upper left side \n"
148 " of AFNI's main window.\n"
149 " 5- Duff : A value of no interest to you. It is preserved for it's historical value.\n"
150 " 6- Delay (Del) : The estimated voxel delay.\n"
151 " 7- Covariance (Cov) : Covariance estimate.\n"
152 " 8- Cross Correlation Coefficient (xCorCoef) : Cross Correlation Coefficient estimate.\n"
153 " 9- Variance (VTS) : Variance of voxel's time series.\n"
154 " This output file can be used as an input to two other plugins:\n"
155 " '4Ddump' and '3D+t Extract'\n\n"
156 " If Write is used, a log file is written too.\n"
157 " The log file contains all parameter settings used for generating the output brick.\n"
158 " The name of the log file is the same as 'Filename' (see below) with a '.log' extension\n"
159 " appended at the end. The log file also holds any warnings generated by the plugin.\n"
160 " Some warnings, such as 'null time series ...' , or 'Could not find zero crossing ...'\n"
161 " are harmless, I might remove them in future versions.\n"
162 " Filename -> Name of the ascii file to write results to.\n"
163 " If the field is left empty, a default name similar \n"
164 " to the default output prefix is used.\n"
165 " Write ts -> (Y/N) Write the time series to an ascii file for voxels with \n"
166 " Cross Correlation Coefficients larger than C-Off.\n"
167 " The file name is that used in 'Filename' with a '.ts' extension appended at the end.\n"
168 " A line (L) in the file 'Filename.ts' contains the time series of the voxel whose\n"
169 " results are written on line (L) in the file 'Filename'.\n"
170 " The time series written to 'Filename.ts' do not contain the ignored samples,\n"
171 " they are detrended and have zero mean.\n\n"
172 "Random Comments/Advice:\n"
173 "Of course, the longer you time series, the better. It is generally recomended that\n"
174 "the largest delay be less than N/10, N being the length of the time series.\n"
175 "The algorithm does go all the way to N/2, but that's not too good.\n\n"
176 "Disclaimer: \n"
177 "(Blaaa bla bla bla bla .... --> legal terms you probably wouldn't read) \n"
178 " I am not responsible for anything bad.\n\n"
179 "If you have/find questions/comments/bugs about the plugin, \n"
180 "send me an E-mail: ziad@image.bien.mu.edu\n\n"
181 " Ziad Saad June 20 97, lastest update Aug 21 00.\n\n"
182 "[1] : Bendat, J. S. (1985). The Hilbert transform and applications to correlation measurements,\n"
183 " Bruel and Kjaer Instruments Inc.\n"
184 "[2] : Bendat, J. S. and G. A. Piersol (1986). Random Data analysis and measurement procedures, \n"
185 " John Wiley & Sons.\n"
186 ;
187
188 /*--------------------- strings for output format --------------------*/
189 /* do not change the order in this string*/
190 static char * method_strings[] = { "Seconds" , "Degrees" , "Radians"} ;
191 static char * yn_strings[] = { "n" , "y" };
192
193 #define NUM_METHOD_STRINGS (sizeof(method_strings)/sizeof(char *))
194 #define NUM_YN_STRINGS (sizeof(yn_strings)/sizeof(char *))
195
196 /* do not change these three*/
197 #define METH_SECONDS 0
198 #define METH_DEGREES 1
199 #define METH_RADIANS 2
200
201 #undef DELAY
202 #define DELAY 0
203 #define XCOR 1
204 #define XCORCOEF 2
205 #ifndef NOWAYXCORCOEF
206 #define NOWAYXCORCOEF 0 /* A flag value indicating that something lethal went on */
207 #endif
208
209 #define NBUCKETS 4 /* Number of values per voxel in Buket data set */
210 #define DELINDX 0 /* index of delay value in results vector */
211 #define COVINDX 1 /* index of covariance value in results vector */
212 #define COFINDX 2 /* index of cross correlation coefficient value in results vector */
213 #define VARINDX 3 /* index of FMRI time course variance value in results vector */
214
215 #define YUP 1
216 #define NOPE 0
217
218
219 /*----------------- prototypes for internal routines -----------------*/
220
221 static char * DELAY_main( PLUGIN_interface * ) ; /* the entry point */
222
223 static void DELAY_tsfuncV2( double T0 , double TR ,
224 int npts , float ts[] , double ts_mean , double ts_slope ,
225 void * udp , int nbrick , float * buckar) ;
226
227 static void show_ud (hilbert_data_V2* ud,int loc);
228
229 static void write_ud (hilbert_data_V2* ud);
230
231 static void indexTOxyz (hilbert_data_V2* ud, int ncall, int *xpos , int *ypos , int *zpos);
232
233 static void error_report (hilbert_data_V2* ud, int ncall );
234
235 static void writets (hilbert_data_V2* ud,float * ts);
236
237 /*---------------------------- global data ---------------------------*/
238
239 static PLUGIN_interface * global_plint = NULL ;
240
241 /***********************************************************************
242 Set up the interface to the user:
243 1) Create a new interface using "PLUTO_new_interface";
244
245 2) For each line of inputs, create the line with "PLUTO_add_option"
246 (this line of inputs can be optional or mandatory);
247
248 3) For each item on the line, create the item with
249 "PLUTO_add_dataset" for a dataset chooser,
250 "PLUTO_add_string" for a string chooser,
251 "PLUTO_add_number" for a number chooser.
252 ************************************************************************/
253
254
255 DEFINE_PLUGIN_PROTOTYPE
256
PLUGIN_init(int ncall)257 PLUGIN_interface * PLUGIN_init( int ncall )
258 {
259 PLUGIN_interface * plint ; /* will be the output of this routine */
260
261 if( ncall > 0 ) return NULL ; /* only one interface */
262 CHECK_IF_ALLOWED("HILBERTDELAY98","Hilbert Delay98") ; /* 30 Sep 2016 */
263
264 /*---------------- set titles and call point ----------------*/
265
266 plint = PLUTO_new_interface( "Hilbert Delay98" ,
267 "Time delay between FMRI and reference time series" ,
268 helpstring ,
269 PLUGIN_CALL_VIA_MENU , DELAY_main ) ;
270
271 global_plint = plint ; /* make global copy */
272
273 /*--------- 1st line: Input dataset ---------*/
274
275 PLUTO_add_option( plint ,
276 "Data" , /* label at left of input line */
277 "Data" , /* tag to return to plugin */
278 TRUE /* is this mandatory? */
279 ) ;
280
281 PLUTO_add_dataset( plint ,
282 "3D+time" , /* label next to button */
283 ANAT_ALL_MASK , /* take only EPI datasets */
284 FUNC_ALL_MASK , /* No fim funcs */
285 DIMEN_4D_MASK | /* need 3D+time datasets */
286 BRICK_ALLREAL_MASK /* need real-valued datasets */
287 ) ;
288
289 PLUTO_add_number( plint ,
290 "Nort" , /* label next to chooser */
291 1 , /* smallest possible value */
292 100 , /* largest possible value (inactivated for now)*/
293 0 , /* decimal shift (none in this case) */
294 2 , /* default value */
295 FALSE /* allow user to edit value? */
296 ) ;
297
298 /*---------- 2nd line: Input time series ----------*/
299
300 PLUTO_add_option( plint ,
301 "Ref." , /* label at left of input line */
302 "Ref." , /* tag to return to plugin */
303 TRUE /* is this mandatory? */
304 ) ;
305
306 PLUTO_add_timeseries(plint,"Ref. Vect.");
307
308 PLUTO_add_number( plint ,
309 "Ignore" , /* label next to chooser */
310 0 , /* smallest possible value */
311 50 , /* largest possible value (inactivated for now)*/
312 0 , /* decimal shift (none in this case) */
313 0 , /* default value */
314 FALSE /* allow user to edit value? */
315 ) ;
316
317 PLUTO_add_string( plint ,
318 "Dsamp" , /*label next to textfield */
319 2,yn_strings, /* strings to choose among */
320 1 /* Default option */
321 ) ;
322
323 /*---------- 3rd line: sampling frequency ----------*/
324
325 PLUTO_add_option( plint ,
326 "Sig." , /* label at left of input line */
327 "Sig." , /* tag to return to plugin */
328 TRUE /* is this mandatory? */
329 ) ;
330
331 PLUTO_add_number( plint ,
332 "fs in Hz" , /* label next to chooser */
333 0 , /* smallest possible value */
334 2000 , /* largest possible value */
335 1 , /* decimal shift (none in this case) */
336 5 , /* default value */
337 TRUE /* allow user to edit value? */
338 ) ;
339
340 PLUTO_add_number( plint ,
341 "Tstim sec" , /* label next to chooser */
342 0.0 , /* smallest possible value */
343 500 , /* largest possible value */
344 0 , /* decimal shift (none in this case) */
345 40 , /* default value */
346 TRUE /* allow user to edit value? */
347 ) ;
348
349 PLUTO_add_number( plint ,
350 "C-Off" , /* label next to chooser */
351 -10 , /* smallest possible value */
352 10 , /* largest possible value */
353 1 , /* decimal shift */
354 5 , /* default value */
355 TRUE /* allow user to edit value? */
356 ) ;
357
358
359 PLUTO_add_string( plint ,
360 "No-bias" , /*label next to textfield */
361 2,yn_strings, /* strings to choose among */
362 1 /* Default option */
363 ) ;
364
365
366
367 /*---------- 4th line: Delay Units ----------*/
368
369 PLUTO_add_option( plint ,
370 "Alg." , /* label at left of input line */
371 "Alg." , /* tag to return to plugin */
372 TRUE /* is this mandatory? */
373 ) ;
374
375 PLUTO_add_number( plint ,
376 "N seg." , /* label next to chooser */
377 1 , /* smallest possible value */
378 1 , /* largest possible value (turned Off for the moment, supporting code is present)*/
379 0 , /* decimal shift (none in this case) */
380 1 , /* default value */
381 FALSE /* allow user to edit value? */
382 ) ;
383
384 PLUTO_add_number( plint ,
385 "% ovrlp" , /* label next to chooser */
386 0 , /* smallest possible value */
387 0 , /* largest possible value (not implemented)*/
388 0 , /* decimal shift (none in this case) */
389 0 , /* default value */
390 FALSE /* allow user to edit value? */
391 ) ;
392
393
394 PLUTO_add_string( plint ,
395 "Units" , /* label next to textfield */
396 3,method_strings, /* strings to choose among */
397 0 /* Default option */
398 ) ;
399
400 PLUTO_add_string( plint ,
401 "Phz Wrp" , /* label next to textfield */
402 2,yn_strings, /* strings to choose among */
403 0 /* Default option */
404 ) ;
405
406
407 /*---------- 5th line: Output dataset ----------*/
408
409 PLUTO_add_option( plint ,
410 "Output" , /* label at left of input line */
411 "Output" , /* tag to return to plugin */
412 TRUE /* is this mandatory? */
413 ) ;
414
415 PLUTO_add_string( plint ,
416 "AFNI Prfx" , /* label next to textfield */
417 0,NULL , /* no fixed strings to choose among */
418 19 /* 19 spaces for typing in value */
419 ) ;
420
421 PLUTO_add_string( plint ,
422 "Write" , /* label next to textfield */
423 2,yn_strings ,
424 1
425 ) ;
426
427 PLUTO_add_string( plint , "Filename" , 0 , NULL , 19 ) ;
428
429 PLUTO_add_string( plint ,
430 "Write ts" , /* label next to textfield */
431 2,yn_strings ,
432 1
433 ) ;
434
435 /*--------- done with interface setup ---------*/
436 return plint ;
437 }
438
439 /***************************************************************************
440 Main routine for this plugin (will be called from AFNI).
441 If the return string is not NULL, some error transpired, and
442 AFNI will popup the return string in a message box.
443 ****************************************************************************/
444
DELAY_main(PLUGIN_interface * plint)445 static char * DELAY_main( PLUGIN_interface * plint )
446 {
447 hilbert_data_V2 uda,*ud;
448 MRI_IMAGE * tsim;
449 MCW_idcode * idc ; /* input dataset idcode */
450 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
451 char *tmpstr , * str , *nprfxstr; /* strings from user */
452 int ntime, nvec ,nprfx, i;
453 float * vec , fs , T ;
454
455 /* Allocate as much character space as Bob specifies in afni.h + a bit more */
456
457 tmpstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
458 nprfxstr = (char *) calloc (PLUGIN_MAX_STRING_RANGE+10,sizeof(char));
459
460 if (tmpstr == NULL || nprfxstr == NULL)
461 return "********************\n"
462 "Could not Allocate\n"
463 "a teeni weeni bit of\n"
464 "Memory ! \n"
465 "********************\n";
466
467 ud = &uda; /* ud now points to an allocated space */
468 ud->errcode = 0; /*reset error flag */
469
470 /*--------------------------------------------------------------------*/
471 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
472
473 /*--------- go to first input line ---------*/
474
475 PLUTO_next_option(plint) ;
476
477 idc = PLUTO_get_idcode(plint) ; /* get dataset item */
478 old_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
479 if( old_dset == NULL )
480 return "*************************\n"
481 "Cannot find Input Dataset\n"
482 "*************************" ;
483
484 ud->dsetname = DSET_FILECODE (old_dset);
485 ud->nsamp = DSET_NUM_TIMES (old_dset);
486 ud->Navg = 1 ; /* Navg does not play a role for the p value, averaging increases sensitivity */
487 ud->Nort = PLUTO_get_number(plint) ; /* Should be two by default, for mean and linear trend */
488 ud->Nfit = 2 ; /* Always 2 for phase and amplitude for this plugin */
489 /*--------- go to 2nd input line, input time series ---------*/
490
491 PLUTO_next_option(plint) ;
492
493 tsim = PLUTO_get_timeseries(plint);
494 if (tsim == NULL) return "No Timeseries Input";
495
496 ud->ln = (int)tsim -> nx; /* number of points in each vector */
497 nvec = tsim -> ny; /* number of vectors */
498 ud->rvec = (float *) MRI_FLOAT_PTR(tsim); /* vec[i+j*nx] = ith point of jth vector */
499 /* for i=0 .. ntime-1 and j=0 .. nvec-1 */
500
501 if (is_vect_null (ud->rvec,ud->ln) == 1) /* check if ref vect is all zeroes */
502 {
503 return "Reference vector is all zeros";
504 }
505
506 ud->refname = tsim->name;
507 ud->ignore = PLUTO_get_number(plint) ; /* get number item */
508
509 str = PLUTO_get_string(plint) ;
510 ud->Dsamp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
511
512 /*--------- go to 3rd input line, sampling frequency, and stimulus period ---------*/
513
514 PLUTO_next_option(plint) ;
515
516 ud->fs = PLUTO_get_number(plint) ; /* get number item */
517 ud->T = PLUTO_get_number(plint) ; /* get number item */
518
519 ud->co = PLUTO_get_number(plint) ; /* get number item */
520 str = PLUTO_get_string(plint) ;
521 ud->biasrem = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
522
523 /*--------- go to 4th input line, delay units and wrp option---------*/
524
525 PLUTO_next_option(plint) ;
526
527 ud->Nseg = (int)PLUTO_get_number(plint) ; /* get number item */
528 ud->Pover = (int)PLUTO_get_number(plint) ; /* get number item */
529
530 str = PLUTO_get_string(plint) ; /* get string item (the method) */
531 ud->unt = (int)PLUTO_string_index( str , /* find it in list it is from */
532 NUM_METHOD_STRINGS ,
533 method_strings ) ;
534
535 str = PLUTO_get_string(plint) ;
536 ud->wrp = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings ) ;
537
538 /*--------- go to 5th input line Output prefix ---------*/
539
540 PLUTO_next_option(plint) ;
541
542 ud->new_prefix = PLUTO_get_string(plint) ; /* get string item (the output prefix) */
543
544 /* check to see if the field is empty */
545 if (ud->new_prefix == NULL)
546 nprfx = 0;
547 else
548 nprfx = 1;
549 /* check if the size is larger than 0. I did not want to check for this unless it's allocated */
550 if (nprfx == 1 && (int)strlen (ud->new_prefix) == 0)
551 nprfx = 0;
552
553 if (nprfx == 0) /* now create the new name and make new_prefix point to it */
554 {
555 sprintf (nprfxstr,"%s.DEL",DSET_PREFIX (old_dset));
556 ud->new_prefix = nprfxstr;
557 /*printf ("New prefix is set to be : %s\n\a",ud->new_prefix);*/
558 }
559
560 if( ! PLUTO_prefix_ok(ud->new_prefix) ) /* check if it is OK */
561 return "************************\n"
562 "Output Prefix is illegal\n"
563 "************************" ;
564
565 str = PLUTO_get_string(plint) ; /* write delays to file ? */
566 ud->out = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
567
568 ud->strout = PLUTO_get_string(plint) ; /* strout is for the outiflename, which will be used after the debugging section */
569 if (ud->strout == NULL) /* if no output name is given, use the new_prefix */
570 {ud->strout = ud->new_prefix;}
571 else
572 {
573 if((int)strlen (ud->strout) == 0) ud->strout = ud->new_prefix;
574 }
575
576 str = PLUTO_get_string(plint) ;
577 ud->outts = (int)PLUTO_string_index( str , NUM_YN_STRINGS , yn_strings );
578
579 /* ------------------Done with user parameters ---------------------------- */
580
581 ud->nxx = (int)old_dset->daxes->nxx; /* get data set dimensions */
582 ud->nyy = (int)old_dset->daxes->nyy;
583 ud->nzz = (int)old_dset->daxes->nzz;
584
585 /* No need for users to set these options ...*/
586
587 ud->dtrnd = 0;
588
589 if (ud->ln != (ud->nsamp - ud->ignore))
590 {
591 ud->errcode = ERROR_BADLENGTH;
592 return "***************************\n"
593 "Bad time series length \n"
594 "Check reference time series\n"
595 " or the ignore parameter \n"
596 "***************************\n";
597 }
598
599 if ((ud->unt < 0) || (ud->unt > 2)) /* unt error Check */
600 {
601 ud->errcode = ERROR_WRONGUNIT;
602 return "***********************\n"
603 " internal error: (ziad)\n"
604 "unt values out of bound\n"
605 "***********************\n"; /*unt must be between 0 and 2 */
606 }
607
608 if ((ud->wrp < 0) || (ud->wrp > 1)) /* wrp error Check */
609 {
610 ud->errcode = ERROR_WARPVALUES;
611 return "***********************\n"
612 " internal error: (ziad)\n"
613 "wrp values out of bound\n"
614 "***********************\n"; /* wrp must be between 0 and 1*/
615 }
616
617 if (ud->fs < 0.0) { /* fs error Check */
618 ud->errcode = ERROR_FSVALUES;
619 return "***********************\n"
620 " internal error: (ziad)\n"
621 "fs value is negative !\n"
622 "***********************\n"; /* fs must be >= 0*/
623 }
624
625 if (ud->T < 0.0) { /* T error Check */
626 ud->errcode = ERROR_TVALUES;
627 return "***********************\n"
628 " internal error: (ziad)\n"
629 "T value is negative !\n"
630 "***********************\n"; /*T must be >= 0 */
631 }
632
633
634 if ((ud->T == 0.0) && (ud->unt > 0)) /* unt error Check */
635 {
636 ud->errcode = ERROR_TaUNITVALUES;
637 return "***********************\n"
638 " internal error: (ziad)\n"
639 "T and unt val. mismatch\n"
640 "***********************\n"; /*T must be specified, and > 0 in order to use polar units*/
641 }
642
643
644 if ((ud->wrp == 1) && (ud->T == 0.0)) /* wrp error Check */
645 {
646 ud->errcode = ERROR_TaWRAPVALUES;
647 return "***********************\n"
648 " internal error: (ziad)\n"
649 "wrp and T val. mismatch\n"
650 "***********************\n"; /*T must be specified, and > 0 in order to use polar warp*/
651 }
652 if ((ud->out == NOPE) && (ud->outts == YUP))
653 {
654 ud->errcode = ERROR_OUTCONFLICT;
655 return"***********************\n"
656 "error: \n"
657 "Write flag must be on\n"
658 "to use Write ts\n"
659 "***********************\n";
660
661 }
662
663
664 /* Open the logfile, regardless of the ascii output files */
665 sprintf ( tmpstr , "%s.log" , ud->strout);
666 ud->outlogfile = fopen (tmpstr,"w");
667
668
669 if (ud->out == YUP) /* open outfile */
670 {
671 ud->outwrite = fopen (ud->strout,"w");
672
673 if (ud->outts == YUP)
674 {
675 sprintf ( tmpstr , "%s.ts" , ud->strout);
676 ud->outwritets = fopen (tmpstr,"w");
677
678 }
679
680 if ((ud->outwrite == NULL) || (ud->outlogfile == NULL) ||\
681 (ud->outwritets == NULL && ud->outts == YUP) )
682 {
683 ud->errcode = ERROR_FILEOPEN;
684
685 return "***********************\n"
686 "Could Not Write Outfile\n"
687 "***********************\n";
688 }
689
690 }
691
692 /* Write out user variables to Logfile */
693 write_ud (ud); /* writes user data to a file */
694
695 /*show_ud (ud,0); */ /* For some debugging */
696
697
698 /*------------- ready to compute new dataset -----------*/
699
700 new_dset = MAKER_4D_to_typed_fbuc ( old_dset , /* input dataset */
701 ud->new_prefix , /* output prefix */
702 -1, /* negative value indicating data type is like original brick */
703 ud->ignore , /* ignore count */
704 1 , /* detrend = ON Let BOB do it*/
705 NBUCKETS, /*Number of values at each voxel*/
706 DELAY_tsfuncV2 , /* timeseries processor (bucket version)*/
707 (void *)ud, /* data for tsfunc */
708 NULL, 0 ) ;
709
710 /* Setup the label, keywords and types of subbricks */
711 i = 0;
712 while (i < NBUCKETS)
713 {
714 switch (i)
715 {
716 case DELINDX: /* delay value in results vector */
717 EDIT_BRICK_LABEL (new_dset,i,"Delay");
718 EDIT_BRICK_ADDKEY (new_dset,i,"D");
719 ++i;
720 break;
721 case COVINDX: /* covariance value in results vector */
722 EDIT_BRICK_LABEL (new_dset,i,"Covariance");
723 EDIT_BRICK_ADDKEY (new_dset,i,"I");
724 ++i;
725 break;
726 case COFINDX: /* cross correlation coefficient value in results vector */
727 EDIT_BRICK_LABEL (new_dset,i,"Corr. Coef.");
728 EDIT_BRICK_ADDKEY (new_dset,i,"r");
729 /* Here you must modify either ud->Nfit or ud->Nort or most likely ud->nsamp based on ud->Navg */
730 EDIT_BRICK_TO_FICO (new_dset,i,ud->nsamp - ud->ignore,ud->Nfit,ud->Nort);
731 ++i;
732 break;
733 case VARINDX: /* FMRI time course variance value in results vector */
734 EDIT_BRICK_LABEL (new_dset,i,"Variance");
735 EDIT_BRICK_ADDKEY (new_dset,i,"S2");
736 ++i;
737 break;
738 default :
739 return "*********************\n"
740 "Internal Error (ziad)\n"
741 " Bad i value \n"
742 "*********************\n";
743 break;
744 }
745
746 }
747
748
749 if (!AFNI_noenv("AFNI_AUTOMATIC_FDR")) {
750 THD_create_all_fdrcurves( new_dset );
751 }
752 PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
753
754
755
756 if (ud->out == YUP) /* close outfile and outlogfile*/
757 {
758 fclose (ud->outlogfile);
759 fclose (ud->outwrite);
760 if (ud->outts == YUP) fclose (ud->outwritets);
761 }
762 else
763 {
764 if (ud->outlogfile != NULL) fclose (ud->outlogfile); /* close outlogfile */
765 }
766
767 free (tmpstr);
768 free (nprfxstr);
769 return NULL ; /* null string returned means all was OK */
770 }
771
772 /**********************************************************************
773 Function that does the real work
774 ***********************************************************************/
775
DELAY_tsfuncV2(double T0,double TR,int npts,float ts[],double ts_mean,double ts_slope,void * udp,int nbrick,float * buckar)776 static void DELAY_tsfuncV2( double T0 , double TR ,
777 int npts , float ts[] , double ts_mean , double ts_slope ,
778 void * udp , int nbrick , float * buckar)
779 {
780 static int nvox , ncall ;
781 hilbert_data_V2 uda,*ud;
782 float del, xcorCoef, buckara[4];
783 float xcor=0.0 , tmp=0.0 , tmp2 = 0.0 , dtx = 0.0 ,\
784 delu = 0.0 , slp = 0.0 , vts = 0.0 , vrvec = 0.0 ;
785 int i , is_ts_null , status , opt , actv , zpos , ypos , xpos ;
786
787 ud = &uda;
788 ud = (hilbert_data_V2 *) udp;
789
790 /** is this a "notification"? **/
791
792 if( buckar == NULL ){
793
794 if( npts > 0 ){ /* the "start notification" */
795
796 PLUTO_popup_meter( global_plint ) ; /* progress meter */
797 nvox = npts ; /* keep track of */
798 ncall = 0 ; /* number of calls */
799
800 } else { /* the "end notification" */
801 opt = 0; /* cleanup in hdelay */
802 #if 0
803 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 */
804 #else
805 hilbertdelay_V2reset();
806 #endif
807 PLUTO_set_meter( global_plint , 100 ) ; /* set meter to 100% */
808
809 }
810 return ;
811 }
812
813 /* In the old version, I had to check for a match between the lengths of the reference time series and FMRI time series
814 This is now done before the function is called. */
815
816 if (is_vect_null (ts,npts) == 1) /* check for null vectors */
817 {
818 ud->errcode = ERROR_NULLTIMESERIES;
819 error_report (ud , ncall ); /* report the error */
820
821 del = 0.0; /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
822 xcorCoef = 0.0; /* because the data might still be OK */
823 xcor = 0.0;
824 }
825
826 if (ud->errcode == 0) /* if there are no errors, proceed */
827 {/* ud->errcode == 0 outer loop */
828 opt = 1; /* activate hdelay */
829
830 /* transform dtx from seconds to sampling units and correct for the number of points ignored*/
831 if (ud->Dsamp == YUP)
832 dtx = (float) (T0 / TR) - ud->ignore;
833 else
834 dtx = 0.0;
835
836 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 */
837
838 if (ud->errcode == 0) /* If there are no errors, proceed */
839 { /*ud->errcode == 0 inner loop */
840 hunwrap (delu, (float)(1/TR), ud->T, slp, ud->wrp, ud->unt,
841 0, 1.0, &del );
842
843 actv = 1; /* assume voxel is active */
844
845 if (xcorCoef < ud->co) actv = 0; /* determine if voxel is activated using xcorCoef */
846
847 if ((actv == 1) && (ud->out == YUP)) /* if voxel is truly activated, write results to file without modifying return value */
848 {
849 indexTOxyz ( ud , ncall, &xpos , &ypos , &zpos);
850 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);
851 if (ud->outts == YUP)
852 {
853 writets (ud,ts);
854 }
855 }
856
857 }/*ud->errcode == 0 inner loop */
858
859 else if (ud->errcode == ERROR_LONGDELAY)
860 {
861 error_report ( ud , ncall);
862
863 del = 0.0; /* Set all the variables to Null and don't set xcorCoef to an impossible value*/
864 xcorCoef = 0.0; /* because the data might still be OK */
865 xcor = 0.0;
866
867 }
868 else if (ud->errcode != 0)
869 {
870 error_report ( ud , ncall);
871
872 del = 0.0; /* Set all the variables to Null and set xcorCoef to an impossible value*/
873 xcorCoef = NOWAYXCORCOEF;
874 xcor = 0.0;
875 }
876
877 }/* ud->errcode == 0 outer loop */
878
879 /* Now fill up the bucket array */
880
881 buckar[DELINDX] = del;
882 buckar[COVINDX] = xcor;
883 buckar[COFINDX] = xcorCoef;
884 buckar[VARINDX] = vts;
885
886
887 /** set the progress meter to the % of completion **/
888 ncall++ ;
889
890 PLUTO_set_meter( global_plint , (100*ncall)/nvox ) ;
891
892 ud->errcode = 0; /* Rest error to nothing */
893
894 return ;
895 }
896
897 /* ************************************************************ */
898 /* function to display user data input (debugging stuff) */
899 /* ************************************************************ */
900
show_ud(hilbert_data_V2 * ud,int loc)901 static void show_ud (hilbert_data_V2* ud,int loc)
902 {
903 printf ("\n\nUser Data Values at location :%d\n",loc);
904 printf ("ud->dsetname= %s\n",ud->dsetname);
905 printf ("ud->refname= %s\n",ud->refname);
906 printf ("ud->rvec= (1st three elements only)");
907 disp_vect (ud->rvec,3);
908 printf ("ud->nxx= %d\n",ud->nxx);
909 printf ("ud->nyy= %d\n",ud->nyy);
910 printf ("ud->nzz= %d\n",ud->nzz);
911 printf ("ud->fs= %f\n",ud->fs);
912 printf ("ud->T= %f\n",ud->T);
913 printf ("ud->co= %f\n",ud->co);
914 printf ("ud->unt= %d\n",ud->unt);
915 printf ("ud->wrp= %d\n",ud->wrp);
916 printf ("ud->Navg= %d\n",ud->Navg);
917 printf ("ud->Nort= %d\n",ud->Nort);
918 printf ("ud->Nfit= %d\n",ud->Nfit);
919 printf ("ud->Nseg= %d\n",ud->Nseg);
920 printf ("ud->Pover= %d\n",ud->Pover);
921 printf ("ud->dtrnd= %d\n",ud->dtrnd);
922 printf ("ud->biasrem= %d\n",ud->biasrem);
923 printf ("ud->Dsamp= %d\n",ud->Dsamp);
924 printf ("ud->ln= %d\n",ud->ln);
925 printf ("ud->nsamp= %d\n",ud->nsamp);
926 printf ("ud->ignore= %d\n",ud->ignore);
927 printf ("ud->errcode= %d\n",ud->errcode);
928 printf ("ud->new_prefix= %s\n",ud->new_prefix);
929 printf ("ud->out= %d\n",ud->out);
930 printf ("ud->strout= %s\n",ud->strout);
931 printf ("ud->outts= %d\n",ud->outts);
932 printf ("Hit enter to continue\a\n\n");
933 getchar ();
934 return;
935 }
936
937 /* ************************************************************ */
938 /* function to write user data input to log file */
939 /* ************************************************************ */
940
write_ud(hilbert_data_V2 * ud)941 static void write_ud (hilbert_data_V2* ud)
942 {
943 fprintf (ud->outlogfile,"\nLogfile output by Hilbert Delay98 plugin\n");
944 fprintf (ud->outlogfile,"\n\nUser Data Values \n");
945 fprintf (ud->outlogfile,"Input data set = %s\n",ud->dsetname);
946 fprintf (ud->outlogfile,"Reference file name = %s\n",ud->refname);
947 fprintf (ud->outlogfile,"Number of voxels in X dimension = %d\n",ud->nxx);
948 fprintf (ud->outlogfile,"Number of voxels in Y dimension = %d\n",ud->nyy);
949 fprintf (ud->outlogfile,"Number of voxels in Z dimension = %d\n",ud->nzz);
950 fprintf (ud->outlogfile,"Sampling Frequency = %f\n",ud->fs);
951 fprintf (ud->outlogfile,"Stimulus Period = %f\n",ud->T);
952 fprintf (ud->outlogfile,"Threshold Cut Off value = %f\n",ud->co);
953 fprintf (ud->outlogfile,"Delay units = %d\n",ud->unt);
954 fprintf (ud->outlogfile,"Delay wrap = %d\n",ud->wrp);
955 fprintf (ud->outlogfile,"Number of segments = %d\n",ud->Nseg);
956 fprintf (ud->outlogfile,"Number of samples in time series = %d\n",ud->nsamp);
957 fprintf (ud->outlogfile,"Ignore = %d\n",ud->ignore);
958 fprintf (ud->outlogfile,"Length of reference time series = %d\n",ud->ln);
959 fprintf (ud->outlogfile,"Number of fit parameters = %d\n",ud->Nfit);
960 fprintf (ud->outlogfile,"Number of nuisance parameters (orts)= %d\n",ud->Nort);
961 fprintf (ud->outlogfile,"Percent overlap = %d\n",ud->Pover);
962 fprintf (ud->outlogfile,"Plugin detrending = %d (Always 0, mandatory detrending is performed)\n",ud->dtrnd);
963 fprintf (ud->outlogfile,"Bias correction = %d\n",ud->biasrem);
964 fprintf (ud->outlogfile,"Acquisition time correction = %d\n",ud->Dsamp);
965 fprintf (ud->outlogfile,"Prefix for birck output = %s\n",ud->new_prefix);
966 fprintf (ud->outlogfile,"Flag for Ascii output file = %d\n",ud->out);
967 fprintf (ud->outlogfile,"Ascii output file name = %s\n",ud->strout);
968 fprintf (ud->outlogfile,"Flag for Ascii time series file output = %d\n",ud->outts);
969 fprintf (ud->outlogfile,"\nud->errcode (debugging only)= %d\n\n",ud->errcode);
970 fprintf (ud->outlogfile,"\nThe format for the output file is the following:\n");
971 fprintf (ud->outlogfile,"VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
972 fprintf (ud->outlogfile,"\nError Log <message> <index> <x> <y> <z>\n\n");
973
974 return;
975 }
976
977 /* ************************************************************ */
978 /* function to compute x, y, z coordinates from the index */
979 /* ************************************************************ */
980
indexTOxyz(hilbert_data_V2 * ud,int ncall,int * xpos,int * ypos,int * zpos)981 static void indexTOxyz (hilbert_data_V2* ud, int ncall, int *xpos , int *ypos , int *zpos)
982 {
983 *zpos = (int)ncall / (int)(ud->nxx*ud->nyy);
984 *ypos = (int)(ncall - *zpos * ud->nxx * ud->nyy) / (int)ud->nxx;
985 *xpos = ncall - ( *ypos * ud->nxx ) - ( *zpos * ud->nxx * ud->nyy ) ;
986 return;
987 }
988
989 /* ************************************************************ */
990 /* function to report errors encountered to the logfile */
991 /* Only errors that happen during runtime (while delays are */
992 /* computed, are handeled here, the others are handeled */
993 /* instantaneously, and need not be logged */
994 /* ************************************************************ */
995
error_report(hilbert_data_V2 * ud,int ncall)996 static void error_report (hilbert_data_V2* ud, int ncall )
997 {
998 int xpos,ypos,zpos;
999 indexTOxyz (ud, ncall, &xpos , &ypos , &zpos);
1000
1001 switch (ud->errcode)
1002 {
1003 case ERROR_NOTHINGTODO:
1004 fprintf (ud->outlogfile,"Nothing to do hilbertdelay_V2 call ");
1005 break;
1006 case ERROR_LARGENSEG:
1007 fprintf (ud->outlogfile,"Number of segments Too Large ");
1008 break;
1009 case ERROR_LONGDELAY:
1010 fprintf (ud->outlogfile,"Could not find zero crossing before maxdel limit ");
1011 break;
1012 case ERROR_SERIESLENGTH:
1013 fprintf (ud->outlogfile,"Vectors have different length ");
1014 break;
1015 case ERROR_NULLTIMESERIES:
1016 fprintf (ud->outlogfile,"Null time series vector ");
1017 break;
1018 default:
1019 fprintf (ud->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",ud->errcode);
1020 break;
1021 }
1022 fprintf (ud->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos );
1023 return;
1024 }
1025
1026 /* *************************************************************** */
1027 /* function to write the time course into a line in the given file */
1028 /* *************************************************************** */
1029
writets(hilbert_data_V2 * ud,float * ts)1030 static void writets (hilbert_data_V2 * ud,float * ts)
1031
1032 {
1033 int i;
1034
1035 for (i=0;i<ud->ln;++i)
1036 {
1037 fprintf (ud->outwritets, "%f\t",ts[i]);
1038 }
1039 fprintf (ud->outwritets,"\n");
1040 }
1041