00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #define PROGRAM_NAME "3ddelay"
00012 #define PROGRAM_AUTHOR "Ziad Saad (using B. Douglas Ward's 3dfim+ to read and write bricks)"
00013 #define PROGRAM_DATE "Jul 22 2005"
00014
00015
00016
00017 #define MAX_FILES 20
00018
00019 #define RA_error FIM_error
00020
00021
00022
00023 #include "mrilib.h"
00024 #include "matrix.h"
00025
00026 #include "delay.c"
00027
00028
00029
00030
00031 static char * method_strings[] = { "Seconds" , "Degrees" , "Radians"} ;
00032 static char * yn_strings[] = { "n" , "y" };
00033
00034
00035 #define CHECK_NULL_STR(str) (str) ? (str) : "(NULL)"
00036
00037
00038 #ifdef ZDBG
00039 #define IPOSx 8
00040 #define IPOSy 38
00041 #define IPOSz 7
00042 #endif
00043
00044 #define NUM_METHOD_STRINGS (sizeof(method_strings)/sizeof(char *))
00045 #define NUM_YN_STRINGS (sizeof(yn_strings)/sizeof(char *))
00046
00047
00048 #define METH_SECONDS 0
00049 #define METH_DEGREES 1
00050 #define METH_RADIANS 2
00051
00052 #undef DELAY
00053 #define DELAY 0
00054 #define XCOR 1
00055 #define XCORCOEF 2
00056 #ifndef NOWAYXCORCOEF
00057 #define NOWAYXCORCOEF 0
00058 #endif
00059
00060
00061 #define NBUCKETS 4
00062 #define DELINDX 0
00063 #define COVINDX 1
00064 #define COFINDX 2
00065 #define VARINDX 3
00066
00067 static char * DELAY_OUTPUT_TYPE_name[NBUCKETS] =
00068 { "Delay", "Covariance", "Corr. Coef.", "Variance" } ;
00069
00070 #define YUP 1
00071 #define NOPE 0
00072
00073 #define ERROR_NOTHINGTODO 1
00074 #define ERROR_LARGENSEG 2
00075 #define ERROR_LONGDELAY 3
00076 #define ERROR_WRONGUNIT 8
00077 #define ERROR_WARPVALUES 9
00078 #define ERROR_FSVALUES 10
00079 #define ERROR_TVALUES 11
00080 #define ERROR_TaUNITVALUES 12
00081 #define ERROR_TaWRAPVALUES 13
00082 #define ERROR_FILEOPEN 15
00083 #define ERROR_SERIESLENGTH 16
00084 #define ERROR_OPTIONS 17
00085 #define ERROR_NULLTIMESERIES 18
00086 #define ERROR_OUTCONFLICT 19
00087 #define ERROR_BADLENGTH 20
00088
00089
00090 typedef struct DELAY_options
00091 {
00092 float fs;
00093
00094
00095 float T;
00096 float co;
00097 int unt;
00098 int wrp;
00099 int Navg;
00100 int Nort;
00101 int Nfit;
00102 int Nseg;
00103 int ignore;
00104 int Pover;
00105 int ln;
00106 int dtrnd;
00107 int biasrem;
00108 int Dsamp;
00109 int errcode;
00110 int out;
00111 int outts;
00112 float *rvec;
00113 int nxx;
00114 int nyy;
00115 int nzz;
00116 FILE * outwrite;
00117 FILE * outwritets;
00118 FILE * outlogfile;
00119
00120 int NFirst;
00121 int NLast;
00122 int N;
00123 int polort;
00124 int num_ortts;
00125 int num_idealts;
00126 int q;
00127 int p;
00128
00129
00130 float fim_thr;
00131 float cdisp;
00132
00133 char * outname;
00134 char * outnamets;
00135 char * outnamelog;
00136
00137 char * input_filename;
00138 char * mask_filename;
00139 char * input1D_filename;
00140
00141 int num_ort_files;
00142 char * ort_filename[MAX_FILES];
00143 int num_ideal_files;
00144 char * ideal_filename[MAX_FILES];
00145 char * bucket_filename;
00146
00147 int output_type[NBUCKETS];
00148
00149 } DELAY_options;
00150
00151
00152
00153
00154
00155
00156
00157 void extract_ts_array
00158 (
00159 THD_3dim_dataset * dset_time,
00160 int iv,
00161 float * ts_array
00162 );
00163
00164
00165
00166
00167
00168
00169 void FIM_error (char * message)
00170 {
00171 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00172 exit(1);
00173 }
00174
00175
00176
00177
00178
00179
00180
00181 void display_help_menu()
00182 {
00183 printf (
00184 "The program estimates the time delay between each voxel time series \n"
00185 "in a 3D+time dataset and a reference time series[1][2]. \n"
00186 "The estimated delays are relative to the reference time series.\n"
00187 "For example, a delay of 4 seconds means that the voxel time series \n"
00188 "is delayed by 4 seconds with respectto the reference time series.\n\n"
00189 " \n"
00190 "Usage: \n"
00191 "3ddelay \n"
00192 "-input fname fname = filename of input 3d+time dataset \n"
00193 "-ideal_file rname rname = input ideal time series file name \n"
00194 " The length of the reference time series should be equal to \n"
00195 " that of the 3d+time data set. \n"
00196 " The reference time series vector is stored in an ascii file. \n"
00197 " The programs assumes that there is one value per line and that all \n"
00198 " values in the file are part of the reference vector. \n"
00199 " PS: Unlike with 3dfim, and FIM in AFNI, values over 33333 are treated\n"
00200 " as part of the time series. \n"
00201 "-fs fs Sampling frequency in Hz. of data time series (1/TR). \n"
00202 "-T Tstim Stimulus period in seconds. \n"
00203 " If the stimulus is not periodic, you can set Tstim to 0.\n"
00204 "[-prefix bucket] The prefix for the results Brick.\n"
00205 " The first subbrick is for Delay.\n"
00206 " The second subbrick is for Covariance, which is an estimate\n"
00207 " of the power in voxel time series at the frequencies present \n"
00208 " in the reference time series.\n"
00209 " The third subbrick is for the Cross Correlation Coefficients between\n"
00210 " FMRI time series and reference time series.\n"
00211 " The fourth subbrick contains estimates of the Variance of voxel time series.\n"
00212 " The default prefix is the prefix of the input 3D+time brick \n"
00213 " with a '.DEL' extension appended to it.\n"
00214 "[-uS/-uD/-uR] Units for delay estimates. (Seconds/Degrees/Radians)\n"
00215 " You can't use Degrees or Radians as units unless \n"
00216 " you specify a value for Tstim > 0.\n"
00217 "[-phzwrp] Delay (or phase) wrap.\n"
00218 " This switch maps delays from: \n"
00219 " (Seconds) 0->T/2 to 0->T/2 and T/2->T to -T/2->0\n"
00220 " (Degrees) 0->180 to 0->180 and 180->360 to -180->0\n"
00221 " (Radians) 0->pi to 0->pi and pi->2pi to -pi->0\n"
00222 " You can't use this option unless you specify a \n"
00223 " value for Tstim > 0.\n\n"
00224 "[-bias] Do not correct for the bias in the estimates [1][2]\n"
00225 "[-mask mname] mname = filename of 3d mask dataset \n"
00226 " only voxels with non-zero values in the mask would be \n"
00227 " considered. \n"
00228 "[-nfirst fnum] fnum = number of first dataset image to use in \n"
00229 " the delay estimate. (default = 0) \n"
00230 "[-nlast lnum] lnum = number of last dataset image to use in \n"
00231 " the delay estimate. (default = last) \n"
00232 "[-nodsamp ] Do not correct a voxel's estimated delay by the time \n"
00233 " at which the slice containing that voxel was acquired.\n\n"
00234 "[-co CCT] Cross Correlation Coefficient threshold value.\n"
00235 " This is only used to limit the ascii output (see below).\n"
00236 "[-nodtrnd] Do not remove the linear trend from the data time series.\n"
00237 " Only the mean is removed. Regardless of this option, \n"
00238 " No detrending is done to the reference time series.\n"
00239 "[-asc [out]] Write the results to an ascii file for voxels with \n"
00240 "[-ascts [out]] cross correlation coefficients larger than CCT.\n"
00241 " If 'out' is not specified, a default name similar \n"
00242 " to the default output prefix is used.\n"
00243 " -asc, only files 'out' and 'out.log' are written to disk (see ahead)\n"
00244 " -ascts, an additional file, 'out.ts', is written to disk (see ahead)\n"
00245 " There a 9 columns in 'out' which hold the following values:\n"
00246 " 1- Voxel Index (VI) : Each voxel in an AFNI brick has a unique index.\n"
00247 " Indices map directly to XYZ coordinates.\n"
00248 " See AFNI plugin documentations for more info.\n"
00249 " 2..4- Voxel coordinates (X Y Z): Those are the voxel slice coordinates.\n"
00250 " You can see these coordinates in the upper left side of the \n"
00251 " AFNI window.To do so, you must first switch the voxel \n"
00252 " coordinate units from mm to slice coordinates. \n"
00253 " Define Datamode -> Misc -> Voxel Coords ?\n"
00254 " PS: The coords that show up in the graph window\n"
00255 " could be different from those in the upper left side \n"
00256 " of AFNI's main window.\n"
00257 " 5- Duff : A value of no interest to you. It is preserved for backward \n"
00258 " compatibility.\n"
00259 " 6- Delay (Del) : The estimated voxel delay.\n"
00260 " 7- Covariance (Cov) : Covariance estimate.\n"
00261 " 8- Cross Correlation Coefficient (xCorCoef) : Cross Correlation Coefficient.\n"
00262 " 9- Variance (VTS) : Variance of voxel's time series.\n\n"
00263 " The file 'out' can be used as an input to two plugins:\n"
00264 " '4Ddump' and '3D+t Extract'\n\n"
00265 " The log file 'out.log' contains all parameter settings used for generating \n"
00266 " the output brick. It also holds any warnings generated by the plugin.\n"
00267 " Some warnings, such as 'null time series ...' , or \n"
00268 " 'Could not find zero crossing ...' are harmless. '\n"
00269 " I might remove them in future versions.\n\n"
00270 " A line (L) in the file 'out.ts' contains the time series of the voxel whose\n"
00271 " results are written on line (L) in the file 'out'.\n"
00272 " The time series written to 'out.ts' do not contain the ignored samples,\n"
00273 " they are detrended and have zero mean.\n\n"
00274 " \n"
00275 "Random Comments/Advice:\n"
00276 " The longer you time series, the better. It is generally recomended that\n"
00277 " the largest delay be less than N/10, N being the length of the time series.\n"
00278 " The algorithm does go all the way to N/2.\n\n"
00279 " If you have/find questions/comments/bugs about the plugin, \n"
00280 " send me an E-mail: ziad@nih.gov\n\n"
00281 " Ziad Saad Dec 8 00.\n\n"
00282 " [1] : Bendat, J. S. (1985). The Hilbert transform and applications to correlation measurements,\n"
00283 " Bruel and Kjaer Instruments Inc.\n"
00284 " [2] : Bendat, J. S. and G. A. Piersol (1986). Random Data analysis and measurement procedures, \n"
00285 " John Wiley & Sons.\n"
00286 " Author's publications on delay estimation using the Hilbert Transform:\n"
00287 " [3] : Saad, Z.S., et al., Analysis and use of FMRI response delays. \n"
00288 " Hum Brain Mapp, 2001. 13(2): p. 74-93.\n"
00289 " [4] : Saad, Z.S., E.A. DeYoe, and K.M. Ropella, Estimation of FMRI Response Delays. \n"
00290 " Neuroimage, 2003. 18(2): p. 494-504.\n\n"
00291 );
00292
00293 exit(0);
00294 }
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 typedef struct {
00308 float real, imag;
00309 } COMPLEX;
00310
00311
00312
00313
00314
00315
00316
00317
00318 extern void dft(COMPLEX *,COMPLEX *,int);
00319 extern void idft(COMPLEX *,COMPLEX *,int);
00320 extern void ham(COMPLEX *,int);
00321 extern void han(COMPLEX *,int);
00322 extern void triang(COMPLEX *,int);
00323 extern void black(COMPLEX *,int);
00324 extern void harris(COMPLEX *,int);
00325
00326 #include "plug_delay_V2.h"
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 void DELAY_tsfuncV2() ;
00340
00341 void show_ud ( DELAY_options* ud,int loc);
00342
00343 void write_ud ( DELAY_options* ud);
00344
00345 void indexTOxyz ( DELAY_options* ud, int ncall, int *xpos , int *ypos , int *zpos);
00346
00347 void xyzTOindex (struct DELAY_options* option_data, int *ncall, int xpos , int ypos , int zpos);
00348
00349 void error_report ( DELAY_options* ud, int ncall );
00350
00351 void writets ( DELAY_options* ud,float * ts);
00352
00353
00354
00355
00356
00357
00358 void initialize_options
00359 (
00360 DELAY_options * option_data
00361 )
00362
00363 {
00364 int is;
00365
00366
00367
00368 option_data->fs = 0;
00369 option_data->co = 0.5;
00370 option_data->T = 0;
00371 option_data->unt = METH_SECONDS;
00372 option_data->wrp = 0;
00373 option_data->Navg = 1;
00374 option_data->Nort = 2;
00375 option_data->Nfit = 2;
00376 option_data->Nseg = 1;
00377 option_data->Pover = 0;
00378 option_data->dtrnd = 1;
00379 option_data->biasrem = 1;
00380 option_data->Dsamp = 1;
00381 option_data->outwrite = NULL;
00382 option_data->outwritets = NULL;
00383 option_data->outlogfile = NULL;
00384 option_data->nxx = 64;
00385 option_data->nyy = 64;
00386 option_data->nzz = 20;
00387 option_data->NFirst = 0;
00388 option_data->NLast = 32767;
00389 option_data->out = 0;
00390 option_data->outts = 0;
00391
00392
00393 option_data->polort = 1;
00394 option_data->num_ortts = 0;
00395 option_data->num_idealts = 0;
00396 option_data->N = 0;
00397 option_data->fim_thr = 0.0999;
00398 option_data->cdisp = -1.0;
00399 option_data->q = 0;
00400 option_data->p = 0;
00401 option_data->num_ort_files = 0;
00402 option_data->num_ideal_files = 0;
00403
00404
00405
00406
00407 option_data->input_filename = NULL;
00408 option_data->mask_filename = NULL;
00409 option_data->input1D_filename = NULL;
00410 option_data->bucket_filename = NULL;
00411 option_data->outname = NULL;
00412
00413 for (is = 0; is < MAX_FILES; is++)
00414 {
00415 option_data->ort_filename[is] = NULL;
00416 option_data->ideal_filename[is] = NULL;
00417 }
00418
00419 }
00420
00421
00422
00423
00424
00425
00426
00427 void get_options
00428 (
00429 int argc,
00430 char ** argv,
00431 DELAY_options * option_data
00432 )
00433
00434 {
00435 int nopt = 1;
00436 int ival, index;
00437 float fval;
00438 char message[THD_MAX_NAME];
00439 int k;
00440
00441
00442
00443 if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu();
00444
00445
00446
00447 initialize_options (option_data);
00448
00449
00450
00451 while (nopt < argc )
00452 {
00453
00454
00455 if (strcmp(argv[nopt], "-input") == 0)
00456 {
00457 nopt++;
00458 if (nopt >= argc) FIM_error ("need argument after -input ");
00459 option_data->input_filename = malloc (sizeof(char)*THD_MAX_NAME);
00460 MTEST (option_data->input_filename);
00461 strcpy (option_data->input_filename, argv[nopt]);
00462 nopt++;
00463 continue;
00464 }
00465
00466
00467 if (strcmp(argv[nopt], "-mask") == 0)
00468 {
00469 nopt++;
00470 if (nopt >= argc) FIM_error ("need argument after -mask ");
00471 option_data->mask_filename = malloc (sizeof(char)*THD_MAX_NAME);
00472 MTEST (option_data->mask_filename);
00473 strcpy (option_data->mask_filename, argv[nopt]);
00474 nopt++;
00475 continue;
00476 }
00477
00478
00479
00480 if (strcmp(argv[nopt], "-nfirst") == 0)
00481 {
00482 nopt++;
00483 if (nopt >= argc) FIM_error ("need argument after -nfirst ");
00484 sscanf (argv[nopt], "%d", &ival);
00485 if (ival < 0)
00486 FIM_error ("illegal argument after -nfirst ");
00487 option_data->NFirst = ival;
00488 nopt++;
00489 continue;
00490 }
00491
00492
00493
00494 if (strcmp(argv[nopt], "-nlast") == 0)
00495 {
00496 nopt++;
00497 if (nopt >= argc) FIM_error ("need argument after -nlast ");
00498 sscanf (argv[nopt], "%d", &ival);
00499 if (ival < 0)
00500 FIM_error ("illegal argument after -nlast ");
00501 option_data->NLast = ival;
00502 nopt++;
00503 continue;
00504 }
00505
00506
00507 if (strcmp(argv[nopt], "-fs") == 0)
00508 {
00509 nopt++;
00510 if (nopt >= argc) FIM_error ("need argument after -fs ");
00511 sscanf (argv[nopt], "%f", &fval);
00512 if (fval < 0)
00513 FIM_error ("illegal argument after -fs ");
00514 option_data->fs = fval;
00515 nopt++;
00516 continue;
00517 }
00518
00519
00520 if (strcmp(argv[nopt], "-T") == 0)
00521 {
00522 nopt++;
00523 if (nopt >= argc) FIM_error ("need argument after -T ");
00524 sscanf (argv[nopt], "%f", &fval);
00525 if (fval < 0)
00526 FIM_error ("illegal argument after -T ");
00527 option_data->T = fval;
00528 nopt++;
00529 continue;
00530 }
00531
00532
00533 if (strcmp(argv[nopt], "-ideal_file") == 0)
00534 {
00535 nopt++;
00536 if (nopt >= argc) FIM_error ("need argument after -ideal_file");
00537
00538 k = option_data->num_ideal_files;
00539 if (k+1 > MAX_FILES)
00540 {
00541 sprintf (message, "Too many ( > %d ) ideal time series files. ",
00542 MAX_FILES);
00543 FIM_error (message);
00544 }
00545
00546 option_data->ideal_filename[k]
00547 = malloc (sizeof(char)*THD_MAX_NAME);
00548 MTEST (option_data->ideal_filename[k]);
00549 strcpy (option_data->ideal_filename[k], argv[nopt]);
00550 option_data->num_ideal_files++;
00551 nopt++;
00552 continue;
00553 }
00554
00555
00556
00557 if (strcmp(argv[nopt], "-prefix") == 0)
00558 {
00559 nopt++;
00560 if (nopt >= argc) FIM_error ("need file prefixname after -bucket ");
00561 option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
00562 MTEST (option_data->bucket_filename);
00563 strcpy (option_data->bucket_filename, argv[nopt]);
00564 nopt++;
00565 continue;
00566 }
00567
00568
00569
00570 if (strcmp(argv[nopt], "-uS") == 0)
00571 {
00572 option_data->unt = METH_SECONDS;
00573 nopt++;
00574 continue;
00575 }
00576
00577
00578 if (strcmp(argv[nopt], "-uR") == 0)
00579 {
00580 option_data->unt = METH_RADIANS;
00581 nopt++;
00582 continue;
00583 }
00584
00585
00586 if (strcmp(argv[nopt], "-uD") == 0)
00587 {
00588 option_data->unt = METH_DEGREES;
00589 nopt++;
00590 continue;
00591 }
00592
00593
00594 if (strcmp(argv[nopt], "-phzwrp") == 0)
00595 {
00596 option_data->wrp = 1;
00597 nopt++;
00598 continue;
00599 }
00600
00601
00602 if (strcmp(argv[nopt], "-bias") == 0)
00603 {
00604 option_data->biasrem = 0;
00605 nopt++;
00606 continue;
00607 }
00608
00609
00610 if (strcmp(argv[nopt], "-nodsamp") == 0)
00611 {
00612 option_data->Dsamp = 0;
00613 nopt++;
00614 continue;
00615 }
00616
00617
00618 if (strcmp(argv[nopt], "-nodtrnd") == 0)
00619 {
00620 option_data->dtrnd = 0;
00621 nopt++;
00622 continue;
00623 }
00624
00625
00626 if (strcmp(argv[nopt], "-co") == 0)
00627 {
00628 nopt++;
00629 if (nopt >= argc) FIM_error ("need argument after -co");
00630 sscanf (argv[nopt], "%f", &fval);
00631 if (fval < 0)
00632 FIM_error ("illegal argument after -co ");
00633 option_data->co = fval;
00634 nopt++;
00635 continue;
00636 }
00637
00638
00639 if (strcmp(argv[nopt], "-asc") == 0)
00640 {
00641 nopt++;
00642 option_data->out = 1;
00643 if (nopt >= argc) {
00644 option_data->outname = NULL;
00645 option_data->outnamelog = NULL;
00646
00647 continue; }
00648 if (strncmp(argv[nopt], "-", 1) == 0) {
00649 option_data->outname = NULL;
00650 option_data->outnamelog = NULL;
00651 continue; }
00652
00653 option_data->outname = malloc (sizeof(char)*THD_MAX_NAME);
00654 option_data->outnamelog = malloc (sizeof(char)*(THD_MAX_NAME+4));
00655
00656 MTEST (option_data->outname);
00657 MTEST (option_data->outnamelog);
00658 strcpy (option_data->outname, argv[nopt]);
00659 sprintf (option_data->outnamelog, "%s.log", option_data->outname);
00660
00661 nopt++;
00662 continue;
00663 }
00664
00665
00666 if (strcmp(argv[nopt], "-ascts") == 0)
00667 {
00668 nopt++;
00669 option_data->out = 1;
00670 option_data->outts = 1;
00671 if (nopt >= argc) {
00672 option_data->outname = NULL;
00673 option_data->outnamelog = NULL;
00674 option_data->outnamets = NULL;
00675 continue; }
00676 if (strncmp(argv[nopt], "-", 1) == 0) {
00677 option_data->outnamelog = NULL;
00678 option_data->outnamets = NULL;
00679 option_data->outname = NULL;
00680 continue; }
00681
00682 option_data->outname = malloc (sizeof(char)*THD_MAX_NAME);
00683 option_data->outnamelog = malloc (sizeof(char)*(THD_MAX_NAME+4));
00684 option_data->outnamets = malloc (sizeof(char)*(THD_MAX_NAME+3));
00685
00686 MTEST (option_data->outname);
00687 MTEST (option_data->outnamets);
00688 MTEST (option_data->outnamelog);
00689
00690 strcpy (option_data->outname, argv[nopt]);
00691 sprintf (option_data->outnamets, "%s.ts", option_data->outname);
00692 sprintf (option_data->outnamelog, "%s.log", option_data->outname);
00693
00694 nopt++;
00695 continue;
00696 }
00697
00698
00699
00700
00701
00702 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00703 FIM_error (message);
00704
00705 }
00706
00707 }
00708
00709
00710
00711
00712
00713
00714
00715
00716 float * read_one_time_series
00717 (
00718 char * ts_filename,
00719 int * ts_length
00720 )
00721
00722 {
00723 char message[THD_MAX_NAME];
00724 char * cpt;
00725 char filename[THD_MAX_NAME];
00726 char subv[THD_MAX_NAME];
00727 MRI_IMAGE * im, * flim;
00728
00729 float * far;
00730 int nx;
00731 int ny;
00732 int iy;
00733 int ipt;
00734 float * ts_data;
00735
00736
00737
00738 flim = mri_read_1D( ts_filename ) ;
00739 if (flim == NULL)
00740 {
00741 sprintf (message, "Unable to read time series file: %s", ts_filename);
00742 FIM_error (message);
00743 }
00744
00745
00746
00747 nx = flim->nx;
00748 ny = flim->ny; iy = 0 ;
00749 if( ny > 1 ){
00750 fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny);
00751 }
00752
00753
00754
00755
00756 *ts_length = nx;
00757 ts_data = (float *) malloc (sizeof(float) * nx);
00758 MTEST (ts_data);
00759 for (ipt = 0; ipt < nx; ipt++)
00760 ts_data[ipt] = far[ipt + iy*nx];
00761
00762
00763 mri_free (flim); flim = NULL;
00764
00765 return (ts_data);
00766 }
00767
00768
00769
00770
00771
00772
00773
00774
00775 MRI_IMAGE * read_time_series
00776 (
00777 char * ts_filename,
00778 int ** column_list
00779 )
00780
00781 {
00782 char message[THD_MAX_NAME];
00783 char * cpt;
00784 char filename[THD_MAX_NAME];
00785 char subv[THD_MAX_NAME];
00786 MRI_IMAGE * im, * flim;
00787
00788 float * far;
00789 int nx;
00790 int ny;
00791
00792
00793
00794
00795 flim = mri_read_1D(ts_filename) ;
00796 if (flim == NULL)
00797 {
00798 sprintf (message, "Unable to read time series file: %s", ts_filename);
00799 FIM_error (message);
00800 }
00801
00802
00803 far = MRI_FLOAT_PTR(flim);
00804 nx = flim->nx;
00805 ny = flim->ny;
00806 *column_list = NULL;
00807
00808 return (flim);
00809 }
00810
00811
00812
00813
00814
00815
00816
00817 void read_input_data
00818 (
00819 DELAY_options * option_data,
00820 THD_3dim_dataset ** dset_time,
00821 THD_3dim_dataset ** mask_dset,
00822 float ** fmri_data,
00823 int * fmri_length,
00824 MRI_IMAGE ** ort_array,
00825 int ** ort_list,
00826 MRI_IMAGE ** ideal_array,
00827 int ** ideal_list
00828 )
00829
00830 {
00831 char message[THD_MAX_NAME];
00832
00833 int num_ort_files;
00834 int num_ideal_files;
00835 int is;
00836 int polort;
00837 int num_ortts;
00838 int num_idealts;
00839 int q;
00840 int p;
00841
00842
00843
00844
00845 polort = option_data->polort;
00846 num_ort_files = option_data->num_ort_files;
00847 num_ideal_files = option_data->num_ideal_files;
00848
00849
00850
00851 if (option_data->input1D_filename != NULL)
00852 {
00853
00854 *fmri_data = read_one_time_series (option_data->input1D_filename,
00855 fmri_length);
00856 if (*fmri_data == NULL)
00857 {
00858 sprintf (message, "Unable to read time series file: %s",
00859 option_data->input1D_filename);
00860 FIM_error (message);
00861 }
00862 *dset_time = NULL;
00863 }
00864
00865 else if (option_data->input_filename != NULL)
00866 {
00867
00868 *dset_time = THD_open_one_dataset (option_data->input_filename);
00869 if (!ISVALID_3DIM_DATASET(*dset_time))
00870 {
00871 sprintf (message, "Unable to open data file: %s",
00872 option_data->input_filename);
00873 FIM_error (message);
00874 }
00875 THD_load_datablock ((*dset_time)->dblk);
00876
00877 if (option_data->mask_filename != NULL)
00878 {
00879
00880 *mask_dset = THD_open_dataset (option_data->mask_filename);
00881 if (!ISVALID_3DIM_DATASET(*mask_dset))
00882 {
00883 sprintf (message, "Unable to open mask file: %s",
00884 option_data->mask_filename);
00885 FIM_error (message);
00886 }
00887 THD_load_datablock ((*mask_dset)->dblk);
00888 }
00889 }
00890 else
00891 FIM_error ("Must specify input measurement data");
00892
00893
00894
00895 for (is = 0; is < num_ideal_files; is++)
00896 {
00897 ideal_array[is] = read_time_series (option_data->ideal_filename[is],
00898 &(ideal_list[is]));
00899
00900 if (ideal_array[is] == NULL)
00901 {
00902 sprintf (message, "Unable to read ideal time series file: %s",
00903 option_data->ideal_filename[is]);
00904 FIM_error (message);
00905 }
00906 }
00907
00908
00909
00910 num_ortts = 0;
00911 for (is = 0; is < num_ort_files; is++)
00912 {
00913 if (ort_list[is] == NULL)
00914 num_ortts += ort_array[is]->ny;
00915 else
00916 num_ortts += ort_list[is][0];
00917 }
00918 q = polort + 1 + num_ortts;
00919
00920 num_idealts = 0;
00921 for (is = 0; is < num_ideal_files; is++)
00922 {
00923 if (ideal_list[is] == NULL)
00924 num_idealts += ideal_array[is]->ny;
00925 else
00926 num_idealts += ideal_list[is][0];
00927 }
00928 p = q + num_idealts;
00929
00930 option_data->num_ortts = num_ortts;
00931 option_data->num_idealts = num_idealts;
00932 option_data->q = q;
00933 option_data->p = p;
00934
00935
00936 }
00937
00938
00939
00940
00941
00942
00943
00944 void check_one_output_file
00945 (
00946 THD_3dim_dataset * dset_time,
00947 char * filename
00948 )
00949
00950 {
00951 char message[THD_MAX_NAME];
00952 THD_3dim_dataset * new_dset=NULL;
00953 int ierror;
00954
00955
00956
00957 new_dset = EDIT_empty_copy( dset_time ) ;
00958
00959
00960 ierror = EDIT_dset_items( new_dset ,
00961 ADN_prefix , filename ,
00962 ADN_label1 , filename ,
00963 ADN_self_name , filename ,
00964 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
00965 GEN_FUNC_TYPE ,
00966 ADN_none ) ;
00967
00968 if( ierror > 0 )
00969 {
00970 sprintf (message,
00971 "*** %d errors in attempting to create output dataset!\n",
00972 ierror);
00973 FIM_error (message);
00974 }
00975
00976 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00977 {
00978 sprintf (message,
00979 "Output dataset file %s already exists "
00980 " -- cannot continue!\a\n",
00981 new_dset->dblk->diskptr->header_name);
00982 FIM_error (message);
00983 }
00984
00985
00986 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00987
00988 }
00989
00990
00991
00992
00993
00994
00995
00996 void check_output_files
00997 (
00998 DELAY_options * option_data,
00999 THD_3dim_dataset * dset_time
01000 )
01001
01002 {
01003
01004 if ((option_data->bucket_filename != NULL)
01005 && (option_data->input1D_filename == NULL))
01006 check_one_output_file (dset_time, option_data->bucket_filename);
01007
01008 }
01009
01010
01011
01012
01013
01014
01015
01016 void check_for_valid_inputs
01017 (
01018 DELAY_options * option_data,
01019 THD_3dim_dataset * dset_time,
01020 THD_3dim_dataset * mask_dset,
01021 int fmri_length,
01022 MRI_IMAGE ** ort_array,
01023 MRI_IMAGE ** ideal_array
01024 )
01025
01026 {
01027 char message[THD_MAX_NAME];
01028 int is;
01029 int num_ort_files;
01030 int num_ideal_files;
01031 int num_idealts;
01032 int nt;
01033 int NFirst;
01034 int NLast;
01035 int N;
01036 int lncheck;
01037
01038
01039 if (option_data->input1D_filename != NULL)
01040 nt = fmri_length;
01041 else
01042 nt = DSET_NUM_TIMES (dset_time);
01043
01044 num_ort_files = option_data->num_ort_files;
01045 num_ideal_files = option_data->num_ideal_files;
01046 num_idealts = option_data->num_idealts;
01047
01048
01049 NFirst = option_data->NFirst;
01050
01051 NLast = option_data->NLast;
01052 if (NLast > nt-1) NLast = nt-1;
01053 option_data->NLast = NLast;
01054
01055 N = NLast - NFirst + 1;
01056 option_data->N = N;
01057
01058
01059
01060 if (num_idealts < 1) FIM_error ("No ideal time series?");
01061
01062
01063
01064 if (mask_dset != NULL)
01065 {
01066 if ( (DSET_NX(dset_time) != DSET_NX(mask_dset))
01067 || (DSET_NY(dset_time) != DSET_NY(mask_dset))
01068 || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) )
01069 {
01070 sprintf (message, "%s and %s have incompatible dimensions",
01071 option_data->input_filename, option_data->mask_filename);
01072 FIM_error (message);
01073 }
01074
01075 if (DSET_NVALS(mask_dset) != 1 )
01076 FIM_error ("Must specify 1 sub-brick from mask dataset");
01077 }
01078
01079
01080
01081
01082
01083 check_output_files (option_data, dset_time);
01084
01085
01086 option_data->ln = option_data->NLast - option_data->NFirst + 1;
01087 option_data->rvec = (float *) malloc (sizeof(float) * option_data->ln+1); MTEST (option_data->rvec);
01088
01089
01090 lncheck = float_file_size (option_data->ideal_filename[0]);
01091 if (lncheck != nt)
01092 {
01093 printf("Error: Reference filename contains %d values.\n %d values were expected.\n", lncheck, nt);
01094 exit (0);
01095 }
01096
01097 Read_part_file_delay (option_data->rvec, option_data->ideal_filename[0], option_data->NFirst + 1,option_data->NLast + 1);
01098
01099
01100
01101 if (option_data->bucket_filename == NULL)
01102 {
01103 option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
01104 MTEST (option_data->bucket_filename);
01105 sprintf (option_data->bucket_filename, "%s.DEL", DSET_PREFIX (dset_time));
01106
01107 check_output_files (option_data, dset_time);
01108 }
01109
01110
01111
01112 if (option_data->outname == NULL)
01113 {
01114 option_data->outnamelog = malloc (sizeof(char)*(THD_MAX_NAME+4));
01115 MTEST (option_data->outnamelog);
01116 sprintf (option_data->outnamelog, "%s.log", option_data->bucket_filename);
01117 }
01118 if (option_data->out || option_data->outts)
01119 {
01120 if (option_data->outname == NULL)
01121 {
01122 option_data->outname = malloc (sizeof(char)*THD_MAX_NAME);
01123 MTEST (option_data->outname);
01124 sprintf (option_data->outname, "%s", option_data->bucket_filename);
01125 }
01126 if (option_data->outts)
01127 {
01128 option_data->outnamets = malloc (sizeof(char)*(THD_MAX_NAME+3));
01129 MTEST (option_data->outnamets);
01130 sprintf (option_data->outnamets, "%s.ts", option_data->outname);
01131 }
01132 }
01133
01134
01135
01136 option_data->outlogfile = fopen (option_data->outnamelog,"w");
01137
01138 if (option_data->out == YUP)
01139 {
01140 option_data->outwrite = fopen (option_data->outname,"w");
01141
01142 if (option_data->outts == YUP)
01143 {
01144 option_data->outwritets = fopen (option_data->outnamets,"w");
01145
01146 }
01147
01148 if ((option_data->outwrite == NULL) || (option_data->outlogfile == NULL) ||\
01149 (option_data->outwritets == NULL && option_data->outts == YUP) )
01150 {
01151 printf ("\nCould not open ascii output files for writing\n");
01152 exit (1);
01153 }
01154
01155 }
01156
01157
01158 write_ud (option_data);
01159
01160 #ifdef ZDBG
01161 show_ud(option_data, 1);
01162 #endif
01163
01164 }
01165
01166
01167
01168
01169
01170
01171
01172 void allocate_memory
01173 (
01174 DELAY_options * option_data,
01175 THD_3dim_dataset * dset,
01176 float *** fim_params_vol
01177 )
01178
01179 {
01180 int ip;
01181 int nxyz;
01182 int ixyz;
01183
01184
01185
01186 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01187
01188
01189
01190 *fim_params_vol = (float **) malloc (sizeof(float *) * NBUCKETS);
01191 MTEST(*fim_params_vol);
01192
01193 for (ip = 0; ip < NBUCKETS; ip++)
01194 {
01195 (*fim_params_vol)[ip] = (float *) malloc (sizeof(float) * nxyz);
01196 MTEST((*fim_params_vol)[ip]);
01197 for (ixyz = 0; ixyz < nxyz; ixyz++)
01198 {
01199 (*fim_params_vol)[ip][ixyz] = 0.0;
01200 }
01201 }
01202
01203 }
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213 void initialize_program
01214 (
01215 int argc,
01216 char ** argv,
01217 DELAY_options ** option_data,
01218 THD_3dim_dataset ** dset_time,
01219 THD_3dim_dataset ** mask_dset,
01220 float ** fmri_data,
01221 int * fmri_length,
01222 MRI_IMAGE ** ort_array,
01223 int ** ort_list,
01224 MRI_IMAGE ** ideal_array,
01225 int ** ideal_list,
01226 float *** fim_params_vol
01227 )
01228
01229 {
01230
01231 *option_data = (DELAY_options *) malloc (sizeof(DELAY_options));
01232
01233
01234
01235 get_options (argc, argv, *option_data);
01236
01237
01238
01239 read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length,
01240 ort_array, ort_list, ideal_array, ideal_list);
01241
01242
01243
01244 check_for_valid_inputs (*option_data, *dset_time, *mask_dset,
01245 *fmri_length, ort_array, ideal_array);
01246
01247
01248
01249 if ((*option_data)->input1D_filename == NULL)
01250 allocate_memory (*option_data, *dset_time, fim_params_vol);
01251
01252 }
01253
01254
01255
01256
01257
01258
01259
01260 void extract_ts_array
01261 (
01262 THD_3dim_dataset * dset_time,
01263 int iv,
01264 float * ts_array
01265 )
01266
01267 {
01268 MRI_IMAGE * im;
01269 float * ar;
01270 int ts_length;
01271 int it;
01272
01273
01274
01275 im = THD_extract_series (iv, dset_time, 0);
01276
01277
01278
01279 if (im == NULL) FIM_error ("Unable to extract data from 3d+time dataset");
01280
01281
01282
01283 ts_length = DSET_NUM_TIMES (dset_time);
01284 ar = MRI_FLOAT_PTR (im);
01285 for (it = 0; it < ts_length; it++)
01286 {
01287 ts_array[it] = ar[it];
01288 }
01289
01290
01291
01292 mri_free (im); im = NULL;
01293
01294 }
01295
01296
01297
01298
01299
01300
01301
01302 void save_voxel
01303 (
01304 int iv,
01305 float * fim_params,
01306 float ** fim_params_vol
01307 )
01308
01309 {
01310 int ip;
01311
01312
01313
01314 for (ip = 0; ip < NBUCKETS; ip++)
01315 {
01316 if (fim_params_vol[ip] != NULL)
01317 fim_params_vol[ip][iv] = fim_params[ip];
01318 }
01319
01320 }
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330 void calculate_results
01331 (
01332 DELAY_options * option_data,
01333 THD_3dim_dataset * dset,
01334 THD_3dim_dataset * mask,
01335 float * fmri_data,
01336 int fmri_length,
01337 MRI_IMAGE ** ort_array,
01338 int ** ort_list,
01339 MRI_IMAGE ** ideal_array,
01340 int ** ideal_list,
01341 float ** fim_params_vol
01342 )
01343
01344 {
01345 float * ts_array = NULL;
01346 float mask_val[1];
01347
01348 int q;
01349 int p;
01350
01351 int m;
01352 int n;
01353
01354
01355 matrix xdata;
01356 matrix x_base;
01357 matrix xtxinvxt_base;
01358 matrix x_ideal[MAX_FILES];
01359 matrix xtxinvxt_ideal[MAX_FILES];
01360
01361 vector y;
01362
01363 int ixyz;
01364 int nxyz;
01365
01366 int nt;
01367 int NFirst;
01368 int NLast;
01369 int N;
01370
01371 int num_ort_files;
01372 int num_ideal_files;
01373 int polort;
01374 int num_ortts;
01375 int num_idealts;
01376
01377 int i;
01378 int is;
01379 int ilag;
01380 float stddev;
01381 char * label;
01382 int xpos, ypos, zpos;
01383 double tzero , tdelta , ts_mean , ts_slope;
01384 int ii , iz,izold, nxy , nuse, use_fac, kk;
01385 float * x_bot = NULL;
01386 float * x_ave = NULL;
01387 float * x_top = NULL;
01388 int * good_list = NULL;
01389 float ** rarray = NULL;
01390 float FimParams[NBUCKETS];
01391
01392 float * dtr = NULL ;
01393 float * fac = NULL ;
01394 float * vox_vect;
01395 float *ref_ts;
01396 float slp, delu, del, xcor, xcorCoef,vts, vrvec, dtx, d0fac , d1fac , x0,x1;
01397 int actv, opt, iposdbg;
01398
01399 #ifdef ZDBG
01400 xyzTOindex (option_data, &iposdbg, IPOSx, IPOSy , IPOSz);
01401 printf ("Debug for %d: %d, %d, %d\n\n", iposdbg, IPOSx, IPOSy , IPOSz);
01402 #endif
01403
01404
01405 matrix_initialize (&xdata);
01406 matrix_initialize (&x_base);
01407 matrix_initialize (&xtxinvxt_base);
01408 for (is =0; is < MAX_FILES; is++)
01409 {
01410 matrix_initialize (&x_ideal[is]);
01411 matrix_initialize (&xtxinvxt_ideal[is]);
01412 }
01413 vector_initialize (&y);
01414
01415
01416
01417 if (option_data->input1D_filename != NULL)
01418 {
01419 nxyz = 1;
01420 nt = fmri_length;
01421 }
01422 else
01423 {
01424 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01425 nt = DSET_NUM_TIMES (dset);
01426 }
01427
01428 NFirst = option_data->NFirst;
01429 NLast = option_data->NLast;
01430 N = option_data->N;
01431 p = option_data->p;
01432 q = option_data->q;
01433
01434 polort = option_data->polort;
01435 num_idealts = option_data->num_idealts;
01436 num_ort_files = option_data->num_ort_files;
01437 num_ideal_files = option_data->num_ideal_files;
01438
01439
01440
01441 ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array);
01442 x_bot = (float *) malloc (sizeof(float) * p); MTEST (x_bot);
01443 x_ave = (float *) malloc (sizeof(float) * p); MTEST (x_ave);
01444 x_top = (float *) malloc (sizeof(float) * p); MTEST (x_top);
01445 good_list = (int *) malloc (sizeof(int) * N); MTEST (good_list);
01446 rarray = (float **) malloc (sizeof(float *) * num_idealts); MTEST (rarray);
01447 vox_vect = (float *) malloc (sizeof(float) * nt); MTEST (vox_vect);
01448
01449
01450 N = init_indep_var_matrix (q, p, NFirst, N, polort,
01451 num_ort_files, num_ideal_files,
01452 ort_array, ort_list, ideal_array, ideal_list,
01453 x_bot, x_ave, x_top, good_list, &xdata);
01454 option_data->N = N;
01455
01456
01457
01458
01459
01460 init_delay (q, p, N, num_idealts, xdata, &x_base,
01461 &xtxinvxt_base, x_ideal, xtxinvxt_ideal, rarray);
01462
01463
01464 vector_create (N, &y);
01465
01466
01467
01468 tdelta = dset->taxis->ttdel ;
01469 if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
01470 if( tdelta == 0.0 ) tdelta = 1.0 ;
01471 izold = -666 ;
01472 nxy = dset->daxes->nxx * dset->daxes->nyy ;
01473
01474
01475
01476 nuse = DSET_NUM_TIMES(dset) - option_data->NFirst ;
01477 fac = (float *) malloc( sizeof(float) * nuse ) ; MTEST (fac);
01478
01479
01480 use_fac = 0 ;
01481 for( kk=0 ; kk < nuse ; kk++ ){
01482 fac[kk] = DSET_BRICK_FACTOR(dset,kk+option_data->NFirst) ;
01483 if( fac[kk] != 0.0 ) use_fac++ ;
01484 else fac[kk] = 1.0 ;
01485 }
01486 if( !use_fac ) FREEUP(fac) ;
01487
01488
01489
01490 dtr = (float *) malloc( sizeof(float) * nuse ) ;MTEST (dtr);
01491
01492
01493 d0fac = 1.0 / nuse ;
01494 d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
01495 for( kk=0 ; kk < nuse ; kk++ )
01496 dtr[kk] = kk - 0.5 * (nuse-1) ;
01497
01498
01499
01500 for (ixyz = 0; ixyz < nxyz; ixyz++)
01501 {
01502 #ifdef ZDBG
01503 if (ixyz == iposdbg)
01504 {
01505 printf ("\nAt voxel, checking for mask\n");
01506 }
01507 #endif
01508
01509 indexTOxyz ( option_data , ixyz, &xpos , &ypos , &zpos);
01510
01511 if (mask != NULL)
01512 {
01513 extract_ts_array (mask, ixyz, mask_val);
01514 if (mask_val[0] == 0.0) continue;
01515 }
01516
01517 #ifdef ZDBG
01518 if (ixyz == iposdbg)
01519 {
01520 printf ("Past the mask, extracting TS\n");
01521 }
01522 #endif
01523
01524
01525
01526 if (option_data->input1D_filename != NULL)
01527 {
01528 for (i = 0; i < N; i++)
01529 vox_vect[i] = (float)fmri_data[good_list[i]+NFirst];
01530 }
01531 else
01532 {
01533 extract_ts_array (dset, ixyz, ts_array);
01534 for (i = 0; i < N; i++)
01535 vox_vect[i] = (float)ts_array[good_list[i]+NFirst];
01536 }
01537
01538 #ifdef ZDBG
01539 if (ixyz == iposdbg)
01540 {
01541 printf ("TS extracted\n");
01542 }
01543 #endif
01544
01545
01546
01547 opt = 1;
01548
01549
01550
01551 #ifdef ZDBG
01552 if (ixyz == iposdbg)
01553 {
01554 printf("Before Scale\n");
01555 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1], vox_vect[2], vox_vect[3], vox_vect[4]);
01556
01557 }
01558 #endif
01559
01560 if( use_fac )
01561 for( kk=0 ; kk < nuse ; kk++ ) vox_vect[kk] *= fac[kk] ;
01562
01563 #ifdef ZDBG
01564 if (ixyz == iposdbg)
01565 {
01566 printf("Before Detrending\n");
01567 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1], vox_vect[2], vox_vect[3], vox_vect[4]);
01568
01569 }
01570 #endif
01571
01572
01573
01574 x0 = x1 = 0.0 ;
01575 for( kk=0 ; kk < nuse ; kk++ ){
01576 x0 += vox_vect[kk] ; x1 += vox_vect[kk] * dtr[kk] ;
01577 }
01578
01579 x0 *= d0fac ; x1 *= d1fac ;
01580
01581 ts_mean = x0 ;
01582 ts_slope = x1 / tdelta ;
01583
01584
01585
01586 if( option_data->dtrnd )
01587 for( kk=0 ; kk < nuse ; kk++ ) vox_vect[kk] -= (x0 + x1 * dtr[kk]) ;
01588 else
01589 for( kk=0 ; kk < nuse ; kk++ ) vox_vect[kk] -= x0;
01590
01591 #ifdef ZDBG
01592 if (ixyz == iposdbg)
01593 {
01594 printf("After Detrending (or just zero-meaning)\n");
01595 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1], vox_vect[2], vox_vect[3], vox_vect[4]);
01596
01597 }
01598 #endif
01599
01600
01601
01602
01603 iz = ixyz / nxy ;
01604
01605 if( iz != izold ){
01606 tzero = THD_timeof( option_data->NFirst ,
01607 dset->daxes->zzorg
01608 + iz*dset->daxes->zzdel , dset->taxis ) ;
01609 izold = iz ;
01610
01611 if( DSET_TIMEUNITS(dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
01612
01613 if (option_data->Dsamp == YUP)
01614 dtx = (float) (tzero / tdelta) - option_data->NFirst;
01615 else
01616 dtx = 0.0;
01617 }
01618
01619 #ifdef ZDBG
01620 if (ixyz == iposdbg)
01621 {
01622 printf("dtx = %f\n", dtx);
01623 }
01624 #endif
01625
01626 option_data->errcode = hilbertdelay_V2 (vox_vect, option_data->rvec,option_data->ln,option_data->Nseg,option_data->Pover,opt,0,dtx,option_data->biasrem,&delu,&slp,&xcor,&xcorCoef,&vts,&vrvec);
01627 #ifdef ZDBG
01628 if (ixyz == iposdbg)
01629 {
01630 printf ("%d err code @ixyz = %d\n", option_data->errcode, ixyz);
01631 }
01632 #endif
01633 if (option_data->errcode == 0)
01634 {
01635 hunwrap (delu, (float)(option_data->fs), option_data->T, slp, option_data->wrp, option_data->unt, &del );
01636
01637 actv = 1;
01638
01639 if (xcorCoef < option_data->co) actv = 0;
01640
01641 if ((actv == 1) && (option_data->out == YUP))
01642 {
01643 indexTOxyz ( option_data , ixyz, &xpos , &ypos , &zpos);
01644 fprintf (option_data->outwrite,"%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\n", ixyz , xpos , ypos , zpos , delu , del , xcor , xcorCoef , vts);
01645 if (option_data->outts == YUP)
01646 {
01647 writets (option_data, vox_vect);
01648 }
01649 }
01650
01651 }
01652
01653 else if (option_data->errcode == ERROR_LONGDELAY)
01654 {
01655 error_report ( option_data , ixyz);
01656
01657 del = 0.0;
01658 xcorCoef = 0.0;
01659 xcor = 0.0;
01660
01661 }
01662 else if (option_data->errcode != 0)
01663 {
01664 error_report ( option_data , ixyz);
01665
01666 del = 0.0;
01667 xcorCoef = NOWAYXCORCOEF;
01668 xcor = 0.0;
01669 }
01670
01671
01672 FimParams[DELINDX] = del;
01673 FimParams[COVINDX] = xcor;
01674 FimParams[COFINDX] = xcorCoef;
01675 FimParams[VARINDX] = vts;
01676
01677 #ifdef ZDBG
01678 if (ixyz == iposdbg)
01679 {
01680 printf("VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
01681 printf("%d\t%d\t%d\t%d\t", ixyz, xpos, ypos, zpos);
01682 printf("%f\t%f\t%f\t%f\t%f\n", delu, del, xcor, xcorCoef, vts);
01683 printf("TS: %f\t%f\t%f\t%f\t%f\n", vox_vect[0], vox_vect[1], vox_vect[2], vox_vect[3], vox_vect[4]);
01684 printf("%f\t%f\t%f\t%f\t%f\n", dtx, delu, del, xcor, xcorCoef);
01685
01686 }
01687 #endif
01688
01689
01690 if (option_data->input1D_filename == NULL)
01691 save_voxel (ixyz, FimParams, fim_params_vol);
01692
01693
01694
01695 }
01696
01697
01698 if (option_data->out == YUP)
01699 {
01700 fclose (option_data->outlogfile);
01701 fclose (option_data->outwrite);
01702 if (option_data->outts == YUP) fclose (option_data->outwritets);
01703 }
01704 else
01705 {
01706 if (option_data->outlogfile != NULL) fclose (option_data->outlogfile);
01707 }
01708
01709
01710
01711 vector_destroy (&y);
01712 for (is = 0; is < MAX_FILES; is++)
01713 {
01714 matrix_destroy (&xtxinvxt_ideal[is]);
01715 matrix_destroy (&x_ideal[is]);
01716 }
01717 matrix_destroy (&xtxinvxt_base);
01718 matrix_destroy (&x_base);
01719 matrix_destroy (&xdata);
01720
01721
01722
01723 free (ts_array); ts_array = NULL;
01724 free (x_bot); x_bot = NULL;
01725 free (x_ave); x_ave = NULL;
01726 free (x_top); x_top = NULL;
01727 free (good_list); good_list = NULL;
01728 free (vox_vect); vox_vect = NULL;
01729
01730 for (is = 0; is < num_idealts; is++)
01731 {
01732 free (rarray[is]); rarray[is] = NULL;
01733 }
01734 free (rarray); rarray = NULL;
01735
01736 }
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750 float EDIT_coerce_autoscale_new( int nxyz ,
01751 int itype,void *ivol , int otype,void *ovol )
01752 {
01753 float fac=0.0 , top ;
01754
01755 if( MRI_IS_INT_TYPE(otype) ){
01756 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01757 if (top == 0.0) fac = 0.0;
01758 else fac = MRI_TYPE_maxval[otype]/top;
01759 }
01760
01761 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01762 return ( fac );
01763 }
01764
01765
01766
01767
01768
01769
01770
01771 void attach_sub_brick
01772 (
01773 THD_3dim_dataset * new_dset,
01774 int ibrick,
01775 float * volume,
01776 int nxyz,
01777 int brick_type,
01778 char * brick_label,
01779 int nsam,
01780 int nfit,
01781 int nort,
01782 short ** bar
01783 )
01784
01785 {
01786 const float EPSILON = 1.0e-10;
01787 float factor;
01788
01789
01790
01791 bar[ibrick] = (short *) malloc (sizeof(short) * nxyz);
01792 MTEST (bar[ibrick]);
01793 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01794 MRI_short, bar[ibrick]);
01795
01796 if (factor < EPSILON) factor = 0.0;
01797 else factor = 1.0 / factor;
01798
01799
01800
01801 EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
01802 EDIT_BRICK_FACTOR (new_dset, ibrick, factor);
01803
01804 if (brick_type == FUNC_COR_TYPE)
01805 EDIT_BRICK_TO_FICO (new_dset, ibrick, nsam, nfit, nort);
01806
01807
01808
01809 EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]);
01810
01811 }
01812
01813
01814
01815
01816
01817
01818 void write_bucket_data
01819 (
01820 int argc,
01821 char ** argv,
01822 DELAY_options * option_data,
01823 float ** fim_params_vol
01824 )
01825
01826 {
01827 THD_3dim_dataset * old_dset = NULL;
01828 THD_3dim_dataset * new_dset = NULL;
01829 char output_prefix[THD_MAX_NAME];
01830 char output_session[THD_MAX_NAME];
01831 int nbricks;
01832 short ** bar = NULL;
01833
01834 int brick_type;
01835 int brick_coef;
01836 char brick_label[THD_MAX_NAME];
01837
01838 int ierror;
01839 float * volume;
01840
01841 int N;
01842 int q;
01843 int num_idealts;
01844 int ip;
01845 int nxyz;
01846 int ibrick;
01847 int nsam;
01848 int nfit;
01849 int nort;
01850 char label[THD_MAX_NAME];
01851
01852
01853
01854 old_dset = THD_open_one_dataset (option_data->input_filename);
01855
01856
01857
01858 nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
01859 num_idealts = option_data->num_idealts;
01860 q = option_data->q;
01861 N = option_data->N;
01862
01863
01864
01865 nbricks = NBUCKETS;
01866
01867 strcpy (output_prefix, option_data->bucket_filename);
01868 strcpy (output_session, "./");
01869
01870
01871
01872 bar = (short **) malloc (sizeof(short *) * nbricks);
01873 MTEST (bar);
01874
01875
01876
01877 new_dset = EDIT_empty_copy (old_dset);
01878
01879
01880
01881 tross_Copy_History( old_dset , new_dset ) ;
01882 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
01883 sprintf (label, "Output prefix: %s", output_prefix);
01884 tross_Append_History ( new_dset, label);
01885
01886
01887
01888 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
01889
01890
01891
01892
01893 ierror = EDIT_dset_items (new_dset,
01894 ADN_prefix, output_prefix,
01895 ADN_directory_name, output_session,
01896 ADN_type, HEAD_FUNC_TYPE,
01897 ADN_func_type, FUNC_BUCK_TYPE,
01898 ADN_datum_all, MRI_short ,
01899 ADN_ntt, 0,
01900 ADN_nvals, nbricks,
01901 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
01902 ADN_none ) ;
01903
01904 if( ierror > 0 )
01905 {
01906 fprintf(stderr,
01907 "*** %d errors in attempting to create bucket dataset!\n",
01908 ierror);
01909 exit(1);
01910 }
01911
01912 if (THD_is_file(DSET_HEADNAME(new_dset)))
01913 {
01914 fprintf(stderr,
01915 "*** Output dataset file %s already exists--cannot continue!\n",
01916 DSET_HEADNAME(new_dset));
01917 exit(1);
01918 }
01919
01920
01921
01922 ibrick = 0;
01923 for (ip = 0; ip < NBUCKETS; ip++)
01924 {
01925 strcpy (brick_label, DELAY_OUTPUT_TYPE_name[ip]);
01926
01927 if (ip == COFINDX)
01928 {
01929 brick_type = FUNC_COR_TYPE;
01930 nsam = option_data->ln; nort = option_data->Nort;
01931 nfit = option_data->Nfit;
01932 }
01933 else
01934 {
01935 brick_type = FUNC_FIM_TYPE;
01936 nsam = 0; nfit = 0; nort = 0;
01937 }
01938
01939 volume = fim_params_vol[ip];
01940 attach_sub_brick (new_dset, ibrick, volume, nxyz,
01941 brick_type, brick_label, nsam, nfit, nort, bar);
01942
01943 ibrick++;
01944 }
01945
01946
01947
01948 THD_load_statistics (new_dset);
01949 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01950 fprintf(stderr,"Wrote bucket dataset ");
01951 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01952
01953
01954
01955 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01956
01957 }
01958
01959
01960
01961
01962
01963
01964
01965 void output_results
01966 (
01967 int argc,
01968 char ** argv,
01969 DELAY_options * option_data,
01970 float ** fim_params_vol
01971 )
01972
01973 {
01974 int q;
01975 int num_idealts;
01976 int ib;
01977 int is;
01978 int ts_length;
01979 int N;
01980
01981
01982
01983 q = option_data->polort + 1;
01984 num_idealts = option_data->num_idealts;
01985 N = option_data->N;
01986
01987
01988
01989 if (option_data->bucket_filename != NULL)
01990 {
01991 write_bucket_data (argc, argv, option_data, fim_params_vol);
01992 }
01993
01994 }
01995
01996
01997
01998
01999
02000 void show_ud (struct DELAY_options* option_data,int loc)
02001 {
02002 printf ("\n\nUser Data Values at location :%d\n",loc);
02003 printf ("option_data->rvec= (1st five elements only)");
02004 disp_vect (option_data->rvec,5);
02005 printf ("Input Data Set: %s\n", option_data->input_filename);
02006 printf ("ASCII output file: %s\n", option_data->outname);
02007 printf ("ASCII output file (log): %s\n", option_data->outnamelog);
02008 printf ("ASCII output file (ts): %s\n", option_data->outnamets);
02009 printf ("Mask Data Set: %s\n", option_data->mask_filename);
02010 printf ("Reference File Name: %s\n", option_data->ideal_filename[0]);
02011 printf ("Number of voxels in X direction: %d\n", option_data->nxx);
02012 printf ("Number of voxels in Y direction: %d\n", option_data->nyy);
02013 printf ("Number of voxels in Z direction: %d\n", option_data->nzz);
02014 printf ("option_data->fs= %f\n",option_data->fs);
02015 printf ("option_data->T= %f\n",option_data->T);
02016 printf ("option_data->unt= %d\n",option_data->unt);
02017 printf ("option_data->wrp= %d\n",option_data->wrp);
02018 printf ("option_data->Navg= %d\n",option_data->Navg);
02019 printf ("option_data->Nort= %d\n",option_data->Nort);
02020 printf ("option_data->Nfit= %d\n",option_data->Nfit);
02021 printf ("option_data->Nseg= %d\n",option_data->Nseg);
02022 printf ("option_data->Pover= %d\n",option_data->Pover);
02023 printf ("option_data->dtrnd= %d\n",option_data->dtrnd);
02024 printf ("option_data->biasrem= %d\n",option_data->biasrem);
02025 printf ("option_data->Dsamp= %d\n",option_data->Dsamp);
02026 printf ("option_data->co= %f\n",option_data->co);
02027 printf ("option_data->ln= %d\n",option_data->ln);
02028 printf ("option_data->errcode= %d\n",option_data->errcode);
02029 printf ("option_data->out= %d\n",option_data->out);
02030 printf ("option_data->outts= %d\n",option_data->outts);
02031 printf ("Hit enter to continue\a\n\n");
02032 getchar ();
02033 return;
02034 }
02035
02036
02037
02038
02039 void write_ud (struct DELAY_options* option_data)
02040 {
02041 fprintf (option_data->outlogfile,"\nLogfile output by Hilbert Delay98 plugin\n");
02042 fprintf (option_data->outlogfile,"\n\nUser Data Values \n");
02043
02044
02045 fprintf (option_data->outlogfile,"Input Data Set: %s\n",
02046 CHECK_NULL_STR(option_data->input_filename));
02047 fprintf (option_data->outlogfile,"Mask Data Set: %s\n",
02048 CHECK_NULL_STR(option_data->mask_filename));
02049 fprintf (option_data->outlogfile,"Ascii output file name: %s\n",
02050 CHECK_NULL_STR(option_data->outname));
02051 fprintf (option_data->outlogfile,"Reference File Name: %s\n",
02052 CHECK_NULL_STR(option_data->ideal_filename[0]));
02053
02054 fprintf (option_data->outlogfile,"Number of voxels in X direction: %d\n", option_data->nxx);
02055 fprintf (option_data->outlogfile,"Number of voxels in Y direction: %d\n", option_data->nyy);
02056 fprintf (option_data->outlogfile,"Number of voxels in Z direction: %d\n", option_data->nzz);
02057 fprintf (option_data->outlogfile,"Sampling Frequency = %f\n",option_data->fs);
02058 fprintf (option_data->outlogfile,"Stimulus Period = %f\n",option_data->T);
02059 fprintf (option_data->outlogfile,"Delay units = %d\n",option_data->unt);
02060 fprintf (option_data->outlogfile,"Delay wrap = %d\n",option_data->wrp);
02061 fprintf (option_data->outlogfile,"Number of segments = %d\n",option_data->Nseg);
02062 fprintf (option_data->outlogfile,"Length of reference time series = %d\n",option_data->ln);
02063 fprintf (option_data->outlogfile,"Number of fit parameters = %d\n",option_data->Nfit);
02064 fprintf (option_data->outlogfile,"Number of nuisance parameters (orts)= %d\n",option_data->Nort);
02065 fprintf (option_data->outlogfile,"Percent overlap = %d\n",option_data->Pover);
02066 fprintf (option_data->outlogfile,"Plugin detrending = %d (Always 0, mandatory detrending is performed)\n",option_data->dtrnd);
02067 fprintf (option_data->outlogfile,"Bias correction = %d\n",option_data->biasrem);
02068 fprintf (option_data->outlogfile,"Acquisition time correction = %d\n",option_data->Dsamp);
02069 fprintf (option_data->outlogfile,"Correlation Coefficient Threshold= %f\n",option_data->co);
02070 fprintf (option_data->outlogfile,"Flag for Ascii output file = %d\n",option_data->out);
02071 fprintf (option_data->outlogfile,"Flag for Ascii time series file output = %d\n",option_data->outts);
02072 fprintf (option_data->outlogfile,"\noption_data->errcode (debugging only)= %d\n\n",option_data->errcode);
02073 fprintf (option_data->outlogfile,"\nThe format for the output file is the following:\n");
02074 fprintf (option_data->outlogfile,"VI\tX\tY\tZ\tDuff\tDel\tCov\txCorCoef\tVTS\n");
02075 fprintf (option_data->outlogfile,"\nError Log <message> <index> <x> <y> <z>\n\n");
02076
02077 return;
02078 }
02079
02080
02081
02082
02083
02084 void indexTOxyz (struct DELAY_options* option_data, int ncall, int *xpos , int *ypos , int *zpos)
02085 {
02086 *zpos = (int)ncall / (int)(option_data->nxx*option_data->nyy);
02087 *ypos = (int)(ncall - *zpos * option_data->nxx * option_data->nyy) / (int)option_data->nxx;
02088 *xpos = ncall - ( *ypos * option_data->nxx ) - ( *zpos * option_data->nxx * option_data->nyy ) ;
02089 return;
02090 }
02091
02092 void xyzTOindex (struct DELAY_options* option_data, int *ncall, int xpos , int ypos , int zpos)
02093 {
02094 *ncall = zpos * ( option_data->nxx*option_data->nyy ) + ypos * option_data->nxx + xpos;
02095 return;
02096 }
02097
02098
02099
02100
02101
02102
02103
02104
02105 void error_report (struct DELAY_options* option_data, int ncall )
02106 {
02107 int xpos,ypos,zpos;
02108 indexTOxyz (option_data, ncall, &xpos , &ypos , &zpos);
02109
02110 switch (option_data->errcode)
02111 {
02112 case ERROR_NOTHINGTODO:
02113 fprintf (option_data->outlogfile,"Nothing to do hilbertdelay_V2 call ");
02114 break;
02115 case ERROR_LARGENSEG:
02116 fprintf (option_data->outlogfile,"Number of segments Too Large ");
02117 break;
02118 case ERROR_LONGDELAY:
02119 fprintf (option_data->outlogfile,"Could not find zero crossing before maxdel limit ");
02120 break;
02121 case ERROR_SERIESLENGTH:
02122 fprintf (option_data->outlogfile,"Vectors have different length ");
02123 break;
02124 case ERROR_NULLTIMESERIES:
02125 fprintf (option_data->outlogfile,"Null time series vector ");
02126 break;
02127 default:
02128 fprintf (option_data->outlogfile,"De Fault, De Fault (%d), the two sweetest words in the english langage ! ",option_data->errcode);
02129 break;
02130 }
02131 fprintf (option_data->outlogfile,"%d\t%d\t%d\t%d\t\n", ncall , xpos , ypos , zpos );
02132 return;
02133 }
02134
02135
02136
02137
02138
02139 void writets (struct DELAY_options * option_data,float * ts)
02140
02141 {
02142 int i;
02143
02144 for (i=0;i<option_data->ln;++i)
02145 {
02146 fprintf (option_data->outwritets, "%f\t",ts[i]);
02147 }
02148 fprintf (option_data->outwritets,"\n");
02149 }
02150
02151
02152
02153 void terminate_program
02154 (
02155 DELAY_options ** option_data,
02156 MRI_IMAGE ** ort_array,
02157 int ** ort_list,
02158 MRI_IMAGE ** ideal_array,
02159 int ** ideal_list,
02160 float *** fim_params_vol
02161 )
02162
02163 {
02164 int num_idealts;
02165 int ip;
02166 int is;
02167
02168
02169
02170 num_idealts = (*option_data)->num_idealts;
02171
02172
02173
02174 free (*option_data); *option_data = NULL;
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185 if (*fim_params_vol != NULL)
02186 {
02187 for (ip = 0; ip < NBUCKETS; ip++)
02188 {
02189 if ((*fim_params_vol)[ip] != NULL)
02190 { free ((*fim_params_vol)[ip]); (*fim_params_vol)[ip] = NULL; }
02191 }
02192
02193 free (*fim_params_vol); *fim_params_vol = NULL;
02194 }
02195
02196
02197
02198 }
02199
02200
02201
02202
02203 int main
02204 (
02205 int argc,
02206 char ** argv
02207 )
02208
02209 {
02210 DELAY_options * option_data;
02211 THD_3dim_dataset * dset_time = NULL;
02212 THD_3dim_dataset * mask_dset = NULL;
02213 float * fmri_data = NULL;
02214 int fmri_length;
02215 MRI_IMAGE * ort_array[MAX_FILES];
02216 int * ort_list[MAX_FILES];
02217 MRI_IMAGE * ideal_array[MAX_FILES];
02218 int * ideal_list[MAX_FILES];
02219
02220 float ** fim_params_vol = NULL;
02221
02222
02223
02224
02225 printf ("\n\n");
02226 printf ("Program: %s \n", PROGRAM_NAME);
02227 printf ("Author: %s \n", PROGRAM_AUTHOR);
02228 printf ("Date: %s \n", PROGRAM_DATE);
02229 printf ("\n");
02230
02231
02232 initialize_program (argc, argv, &option_data, &dset_time, &mask_dset,
02233 &fmri_data, &fmri_length,
02234 ort_array, ort_list, ideal_array, ideal_list,
02235 &fim_params_vol);
02236
02237
02238
02239 calculate_results (option_data, dset_time, mask_dset,
02240 fmri_data, fmri_length,
02241 ort_array, ort_list, ideal_array, ideal_list,
02242 fim_params_vol);
02243
02244
02245
02246 if (dset_time != NULL)
02247 { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; }
02248 if (mask_dset != NULL)
02249 { THD_delete_3dim_dataset (mask_dset, False); mask_dset = NULL; }
02250
02251
02252
02253 if (option_data->input1D_filename == NULL)
02254 output_results (argc, argv, option_data, fim_params_vol);
02255
02256
02257
02258 terminate_program (&option_data,
02259 ort_array, ort_list, ideal_array, ideal_list,
02260 &fim_params_vol);
02261
02262 }
