1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 Program to calculate the cross-correlation of an ideal reference waveform
10 with the measured FMRI time series for each voxel.
11 This is the stand-alone (batch command) version of the afni fim+ function.
12
13 File: 3dfim+.c
14 Author: B. Douglas Ward
15 Date: 28 April 2000
16
17
18 Mod: Use output_type array in regression_analysis routine to avoid
19 some unnecessary calculations.
20 Date: 18 May 2000
21
22 Mod: Added call to AFNI_logger.
23 Date: 15 August 2001
24
25 Mod: Add FICO-ness of sub-bricks for Spearman and Quadrant correlation.
26 Date 29 Oct 2004 - RWCox
27
28 */
29
30 /*---------------------------------------------------------------------------*/
31
32 #define PROGRAM_NAME "3dfim+" /* name of this program */
33 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
34 #define PROGRAM_INITIAL "28 April 2000" /* date of initial program release */
35 #define PROGRAM_LATEST "29 October 2004" /* date of latest program revision */
36
37 /*---------------------------------------------------------------------------*/
38
39 #define MAX_FILES 200 /* maximum number of ideal files */
40 /* = maximum number of ort files */
41 /* It looks like this is also the
42 limit on the number of regressors
43 (columns) in an ideal file
44 ZSS May 15 08 */
45 #define RA_error FIM_error
46
47 /*---------------------------------------------------------------------------*/
48
49 #include "mrilib.h"
50 #include "matrix.h"
51
52 #include "fim+.c"
53
54 /*---------------------------------------------------------------------------*/
55
56 typedef struct FIM_options
57 {
58 int NFirst; /* first image from input 3d+time dataset to use */
59 int NLast; /* last image from input 3d+time dataset to use */
60 int N; /* number of usable data points from input data */
61 int polort; /* degree of polynomial for baseline model */
62 int num_ortts; /* number of ort time series */
63 int num_idealts; /* number of ideal time series */
64 int q; /* number of parameters in the baseline model */
65 int p; /* number of parameters in the baseline model
66 plus number of ideals */
67
68 float fim_thr; /* threshold for internal fim mask */
69 float cdisp; /* minimum correlation coefficient for display */
70
71 char * input_filename; /* input 3d+time dataset filename */
72 char * mask_filename; /* input mask dataset filename */
73 char * input1D_filename; /* input fMRI measurement time series */
74
75 int num_ort_files; /* number of ort files */
76 char * ort_filename[MAX_FILES]; /* input ort time series file names */
77 int num_ideal_files; /* number of ideal files */
78 char * ideal_filename[MAX_FILES]; /* input ideal time series file names */
79 char * bucket_filename; /* output bucket dataset file name */
80
81 int output_type[MAX_OUTPUT_TYPE]; /* output type options */
82
83 } FIM_options;
84
85
86 /*---------------------------------------------------------------------------*/
87 /*
88 Get the time series for one voxel from the AFNI 3d+time data set.
89 */
90
91 void extract_ts_array
92 (
93 THD_3dim_dataset * dset_time, /* input 3d+time dataset */
94 int iv, /* get time series for this voxel */
95 float * ts_array /* time series data for voxel #iv */
96 );
97
98 /*---------------------------------------------------------------------------*/
99 /*
100 Print error message and stop.
101 */
102
FIM_error(char * message)103 void FIM_error (char * message)
104 {
105 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
106 exit(1);
107 }
108
109
110 /*---------------------------------------------------------------------------*/
111 /*
112 Routine to display 3dfim+ help menu.
113 */
114
display_help_menu()115 void display_help_menu()
116 {
117 printf (
118 "Program to calculate the cross-correlation of an ideal reference waveform \n"
119 "with the measured FMRI time series for each voxel. \n"
120 " \n"
121 "Usage: \n"
122 "3dfim+ \n"
123 "-input fname fname = filename of input 3d+time dataset \n"
124 "[-input1D dname] dname = filename of single (fMRI) .1D time series \n"
125 "[-mask mname] mname = filename of 3d mask dataset \n"
126 "[-nfirst fnum] fnum = number of first dataset image to use in \n"
127 " the cross-correlation procedure. (default = 0) \n"
128 "[-nlast lnum] lnum = number of last dataset image to use in \n"
129 " the cross-correlation procedure. (default = last) \n"
130 "[-polort pnum] pnum = degree of polynomial corresponding to the \n"
131 " baseline model (pnum = 0, 1, etc.) \n"
132 " (default: pnum = 1). Use -1 for no baseline model.\n"
133 "[-fim_thr p] p = fim internal mask threshold value (0 <= p <= 1) \n"
134 " to get rid of low intensity voxels. \n"
135 " (default: p = 0.0999), set p = 0.0 for no masking.\n"
136 "[-cdisp cval] Write (to screen) results for those voxels \n"
137 " whose correlation stat. > cval (0 <= cval <= 1) \n"
138 " (default: disabled) \n"
139 "[-ort_file sname] sname = input ort time series file name \n"
140 "-ideal_file rname rname = input ideal time series file name \n"
141 " \n"
142 " Note: The -ort_file and -ideal_file commands may be used \n"
143 " more than once. \n"
144 " Note: If files sname or rname contain multiple columns, \n"
145 " then ALL columns will be used as ort or ideal \n"
146 " time series. However, individual columns or \n"
147 " a subset of columns may be selected using a file \n"
148 " name specification like 'fred.1D[0,3,5]', which \n"
149 " indicates that only columns #0, #3, and #5 will \n"
150 " be used for input. \n"
151 );
152
153 printf("\n"
154 "[-out param] Flag to output the specified parameter, where \n"
155 " the string 'param' may be any one of the following: \n"
156 " \n"
157 "%12s L.S. fit coefficient for Best Ideal \n"
158 "%12s Index number for Best Ideal (count starts at 1) \n"
159 "%12s P-P amplitude of signal response / Baseline \n"
160 "%12s Average of baseline model response \n"
161 "%12s Best Ideal product-moment correlation coefficient \n"
162 "%12s P-P amplitude of signal response / Average \n"
163 "%12s Baseline + average of signal response \n"
164 "%12s P-P amplitude of signal response / Topline \n"
165 "%12s Baseline + P-P amplitude of signal response \n"
166 "%12s Std. Dev. of residuals from best fit \n"
167 "%9sAll This specifies all of the above parameters \n"
168 "%12s Spearman correlation coefficient \n"
169 "%12s Quadrant correlation coefficient \n"
170 " \n"
171 " Note: Multiple '-out' commands may be used. \n"
172 " Note: If a parameter name contains imbedded spaces, the \n"
173 " entire parameter name must be enclosed by quotes, \n"
174 " e.g., -out '%8s' \n"
175 " \n"
176 "[-bucket bprefix] Create one AFNI 'bucket' dataset containing the \n"
177 " parameters of interest, as specified by the above \n"
178 " '-out' commands. \n"
179 " The output 'bucket' dataset is written to a file \n"
180 " with the prefix name bprefix. \n"
181 ,
182 OUTPUT_TYPE_name[FIM_FitCoef],
183 OUTPUT_TYPE_name[FIM_BestIndex],
184 OUTPUT_TYPE_name[FIM_PrcntChange],
185 OUTPUT_TYPE_name[FIM_Baseline],
186 OUTPUT_TYPE_name[FIM_Correlation],
187 OUTPUT_TYPE_name[FIM_PrcntFromAve],
188 OUTPUT_TYPE_name[FIM_Average],
189 OUTPUT_TYPE_name[FIM_PrcntFromTop],
190 OUTPUT_TYPE_name[FIM_Topline],
191 OUTPUT_TYPE_name[FIM_SigmaResid],
192 "",
193 OUTPUT_TYPE_name[FIM_SpearmanCC],
194 OUTPUT_TYPE_name[FIM_QuadrantCC],
195 OUTPUT_TYPE_name[FIM_FitCoef]
196 );
197
198 PRINT_COMPILE_DATE ; exit(0);
199 }
200
201
202 /*---------------------------------------------------------------------------*/
203 /*
204 Routine to initialize the input options.
205 */
206
initialize_options(FIM_options * option_data)207 void initialize_options
208 (
209 FIM_options * option_data /* fim program options */
210 )
211
212 {
213 int is; /* index */
214
215
216 /*----- Initialize default values -----*/
217 option_data->NFirst = 0;
218 option_data->NLast = 32767;
219 option_data->N = 0;
220 option_data->polort = 1;
221 option_data->num_ortts = 0;
222 option_data->num_idealts = 0;
223 option_data->q = 0;
224 option_data->p = 0;
225
226 option_data->fim_thr = 0.0999;
227 option_data->cdisp = -1.0;
228
229 option_data->num_ort_files = 0;
230 option_data->num_ideal_files = 0;
231
232
233 /*----- Initialize output flags -----*/
234 for (is = 0; is < MAX_OUTPUT_TYPE; is++)
235 option_data->output_type[is] = 0;
236
237
238 /*----- Initialize file names -----*/
239 option_data->input_filename = NULL;
240 option_data->mask_filename = NULL;
241 option_data->input1D_filename = NULL;
242 option_data->bucket_filename = NULL;
243
244 for (is = 0; is < MAX_FILES; is++)
245 {
246 option_data->ort_filename[is] = NULL;
247 option_data->ideal_filename[is] = NULL;
248 }
249
250 }
251
252
253 /*---------------------------------------------------------------------------*/
254 /*
255 Routine to get user specified input options.
256 */
257
get_options(int argc,char ** argv,FIM_options * option_data)258 void get_options
259 (
260 int argc, /* number of input arguments */
261 char ** argv, /* array of input arguments */
262 FIM_options * option_data /* fim program options */
263 )
264
265 {
266 int nopt = 1; /* input option argument counter */
267 int ival, index; /* integer input */
268 float fval; /* float input */
269 char message[THD_MAX_NAME]; /* error message */
270 int k; /* ideal time series index */
271
272
273 /*----- does user request help menu? -----*/
274 if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu();
275
276
277 /*----- add to program log -----*/
278 AFNI_logger (PROGRAM_NAME,argc,argv);
279
280
281 /*----- initialize the input options -----*/
282 initialize_options (option_data);
283
284
285 /*----- main loop over input options -----*/
286 while (nopt < argc )
287 {
288
289 /*----- -input filename -----*/
290 if (strcmp(argv[nopt], "-input") == 0)
291 {
292 nopt++;
293 if (nopt >= argc) FIM_error ("need argument after -input ");
294 option_data->input_filename = malloc (sizeof(char)*THD_MAX_NAME);
295 MTEST (option_data->input_filename);
296 strcpy (option_data->input_filename, argv[nopt]);
297 nopt++;
298 continue;
299 }
300
301
302 /*----- -mask filename -----*/
303 if (strcmp(argv[nopt], "-mask") == 0)
304 {
305 nopt++;
306 if (nopt >= argc) FIM_error ("need argument after -mask ");
307 option_data->mask_filename = malloc (sizeof(char)*THD_MAX_NAME);
308 MTEST (option_data->mask_filename);
309 strcpy (option_data->mask_filename, argv[nopt]);
310 nopt++;
311 continue;
312 }
313
314
315 /*----- -input1D filename -----*/
316 if (strcmp(argv[nopt], "-input1D") == 0)
317 {
318 nopt++;
319 if (nopt >= argc) FIM_error ("need argument after -input1D ");
320 option_data->input1D_filename =
321 malloc (sizeof(char)*THD_MAX_NAME);
322 MTEST (option_data->input1D_filename);
323 strcpy (option_data->input1D_filename, argv[nopt]);
324 nopt++;
325 continue;
326 }
327
328
329 /*----- -nfirst num -----*/
330 if (strcmp(argv[nopt], "-nfirst") == 0)
331 {
332 nopt++;
333 if (nopt >= argc) FIM_error ("need argument after -nfirst ");
334 sscanf (argv[nopt], "%d", &ival);
335 if (ival < 0)
336 FIM_error ("illegal argument after -nfirst ");
337 option_data->NFirst = ival;
338 nopt++;
339 continue;
340 }
341
342
343 /*----- -nlast num -----*/
344 if (strcmp(argv[nopt], "-nlast") == 0)
345 {
346 nopt++;
347 if (nopt >= argc) FIM_error ("need argument after -nlast ");
348 sscanf (argv[nopt], "%d", &ival);
349 if (ival < 0)
350 FIM_error ("illegal argument after -nlast ");
351 option_data->NLast = ival;
352 nopt++;
353 continue;
354 }
355
356
357 /*----- -polort num -----*/
358 if (strcmp(argv[nopt], "-polort") == 0)
359 {
360 char *qpt ;
361 nopt++;
362 if (nopt >= argc) FIM_error ("need argument after -polort ");
363 ival = (int)strtod(argv[nopt],&qpt) ;
364 if( *qpt != '\0' ) WARNING_message("Illegal non-numeric value after -polort") ;
365
366 #undef PMAX
367 #ifdef USE_LEGENDRE
368 # define PMAX 19
369 #else
370 # define PMAX 2
371 #endif
372
373 if ((ival < -1) || (ival > PMAX)) /* ZSS May 08, allowed -1 */
374 FIM_error ("illegal argument after -polort ");
375
376 #ifdef USE_LEGENDRE
377 if( 0 && ival > 2 ) /* ZSS June 2010, No need for fear mongering */
378 fprintf(stderr,
379 "** WARNING: -polort > 2 is a new untested option: 29 Mar 2005\n") ;
380 #endif
381
382 option_data->polort = ival;
383 nopt++;
384 continue;
385 }
386
387
388 /*----- -fim_thr r -----*/
389 if (strcmp(argv[nopt], "-fim_thr") == 0)
390 {
391 nopt++;
392 if (nopt >= argc) FIM_error ("need argument after -fim_thr ");
393 sscanf (argv[nopt], "%f", &fval);
394 if ((fval < 0.0) || (fval > 1.0))
395 FIM_error ("illegal argument after -fim_thr ");
396 option_data->fim_thr = fval;
397 nopt++;
398 continue;
399 }
400
401
402 /*----- -cdisp cval -----*/
403 if (strcmp(argv[nopt], "-cdisp") == 0)
404 {
405 nopt++;
406 if (nopt >= argc) FIM_error ("need argument after -cdisp ");
407 sscanf (argv[nopt], "%f", &fval);
408 if ((fval < 0.0) || (fval > 1.0))
409 FIM_error ("illegal argument after -cdisp ");
410 option_data->cdisp = fval;
411 nopt++;
412 continue;
413 }
414
415
416 /*----- -ort_file sname -----*/
417 if (strcmp(argv[nopt], "-ort_file") == 0)
418 {
419 nopt++;
420 if (nopt >= argc) FIM_error ("need argument after -ort_file");
421
422 k = option_data->num_ort_files;
423 if (k+1 > MAX_FILES)
424 {
425 sprintf (message, "Too many ( > %d ) ort time series files. ",
426 MAX_FILES);
427 FIM_error (message);
428 }
429
430 option_data->ort_filename[k]
431 = malloc (sizeof(char)*THD_MAX_NAME);
432 MTEST (option_data->ort_filename[k]);
433 strcpy (option_data->ort_filename[k], argv[nopt]);
434 option_data->num_ort_files++;
435 nopt++;
436 continue;
437 }
438
439
440 /*----- -ideal_file rname -----*/
441 if (strcmp(argv[nopt], "-ideal_file") == 0)
442 {
443 nopt++;
444 if (nopt >= argc) FIM_error ("need argument after -ideal_file");
445
446 k = option_data->num_ideal_files;
447 if (k+1 > MAX_FILES)
448 {
449 sprintf (message, "Too many ( > %d ) ideal time series files. ",
450 MAX_FILES);
451 FIM_error (message);
452 }
453
454 option_data->ideal_filename[k]
455 = malloc (sizeof(char)*THD_MAX_NAME);
456 MTEST (option_data->ideal_filename[k]);
457 strcpy (option_data->ideal_filename[k], argv[nopt]);
458 option_data->num_ideal_files++;
459 nopt++;
460 continue;
461 }
462
463
464 /*----- -out option -----*/
465 if (strcmp(argv[nopt], "-out") == 0)
466 {
467 nopt++;
468 if (nopt >= argc) FIM_error ("need argument after -out ");
469 for (ival = 0; ival < MAX_OUTPUT_TYPE; ival++)
470 if (strcmp(argv[nopt], OUTPUT_TYPE_name[ival]) == 0) break;
471 if (ival < MAX_OUTPUT_TYPE)
472 {
473 option_data->output_type[ival] = 1;
474 nopt++;
475 continue;
476 }
477 else if (strcmp(argv[nopt], "All") == 0)
478 {
479 for (ival = 0; ival < MAX_OUTPUT_TYPE-2; ival++)
480 option_data->output_type[ival] = 1;
481 nopt++;
482 continue;
483 }
484 else
485 {
486 sprintf(message,"Unrecognized output type: %s\n", argv[nopt]);
487 FIM_error (message);
488 }
489 }
490
491
492 /*----- -bucket filename -----*/
493 if (strcmp(argv[nopt], "-bucket") == 0)
494 {
495 nopt++;
496 if (nopt >= argc) FIM_error ("need file prefixname after -bucket ");
497 option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
498 MTEST (option_data->bucket_filename);
499 strcpy (option_data->bucket_filename, argv[nopt]);
500 nopt++;
501 continue;
502 }
503
504
505 /*----- unknown command -----*/
506 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
507 FIM_error (message);
508
509 }
510
511 }
512
513
514 /*---------------------------------------------------------------------------*/
515 /*
516 Read one time series from specified file name. This file name may have
517 a column selector attached.
518 */
519
read_one_time_series(char * ts_filename,int * ts_length)520 float * read_one_time_series
521 (
522 char * ts_filename, /* time series file name (plus column index) */
523 int * ts_length /* output value for time series length */
524 )
525
526 {
527 char message[THD_MAX_NAME]; /* error message */
528 char * cpt; /* pointer to column suffix */
529 char subv[THD_MAX_NAME]; /* string containing column index */
530 MRI_IMAGE * im, * flim; /* pointers to image structures
531 -- used to read 1D ASCII */
532 float * far; /* pointer to MRI_IMAGE floating point data */
533 int nx; /* number of time points in time series */
534 int ny; /* number of columns in time series file */
535 int iy; /* time series file column index */
536 int ipt; /* time point index */
537 float * ts_data; /* input time series data */
538
539
540 /*----- Read the time series file -----*/
541 flim = mri_read_1D(ts_filename) ;
542 if (flim == NULL)
543 { /* filename -> ts_filename 11 Sep 2003 [rickr] */
544 sprintf (message, "Unable to read time series file: %s", ts_filename);
545 FIM_error (message);
546 }
547
548
549 /*----- Set pointer to data, and set dimensions -----*/
550 far = MRI_FLOAT_PTR(flim);
551 nx = flim->nx;
552 ny = flim->ny; iy = 0 ;
553 if( ny > 1 ){
554 fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny);
555 }
556
557
558 /*----- Save the time series data -----*/
559 *ts_length = nx;
560 ts_data = (float *) malloc (sizeof(float) * nx);
561 MTEST (ts_data);
562 for (ipt = 0; ipt < nx; ipt++)
563 ts_data[ipt] = far[ipt + iy*nx];
564
565
566 mri_free (flim); flim = NULL;
567
568 return (ts_data);
569 }
570
571
572 /*---------------------------------------------------------------------------*/
573 /*
574 Read multiple time series from specified file name. This file name may have
575 a column selector attached.
576 */
577
read_time_series(char * ts_filename,int ** column_list)578 MRI_IMAGE * read_time_series
579 (
580 char * ts_filename, /* time series file name (plus column selectors) */
581 int ** column_list /* list of selected columns */
582 )
583
584 {
585 char message[THD_MAX_NAME]; /* error message */
586 char * cpt; /* pointer to column suffix */
587 char filename[THD_MAX_NAME]; /* time series file name w/o column index */
588 char subv[THD_MAX_NAME]; /* string containing column index */
589 MRI_IMAGE * im, * flim; /* pointers to image structures
590 -- used to read 1D ASCII */
591 float * far; /* pointer to MRI_IMAGE floating point data */
592 int nx; /* number of time points in time series */
593 int ny; /* number of columns in time series file */
594
595
596 /*----- Read the time series file -----*/
597 flim = mri_read_1D(ts_filename) ;
598 if (flim == NULL)
599 {
600 sprintf (message, "Unable to read time series file: %s", ts_filename);
601 FIM_error (message);
602 }
603
604
605 /*----- Set pointer to data, and set dimensions -----*/
606 far = MRI_FLOAT_PTR(flim);
607 nx = flim->nx;
608 ny = flim->ny;
609 *column_list = NULL; /* mri_read_1D does column selection */
610
611 return (flim);
612 }
613
614
615 /*---------------------------------------------------------------------------*/
616 /*
617 Read the input data files.
618 */
619
read_input_data(FIM_options * option_data,THD_3dim_dataset ** dset_time,THD_3dim_dataset ** mask_dset,float ** fmri_data,int * fmri_length,MRI_IMAGE ** ort_array,int ** ort_list,MRI_IMAGE ** ideal_array,int ** ideal_list)620 void read_input_data
621 (
622 FIM_options * option_data, /* fim program options */
623 THD_3dim_dataset ** dset_time, /* input 3d+time data set */
624 THD_3dim_dataset ** mask_dset, /* input mask data set */
625 float ** fmri_data, /* input fMRI time series data */
626 int * fmri_length, /* length of fMRI time series */
627 MRI_IMAGE ** ort_array, /* ort time series arrays */
628 int ** ort_list, /* list of ort time series */
629 MRI_IMAGE ** ideal_array, /* ideal time series arrays */
630 int ** ideal_list /* list of ideal time series */
631 )
632
633 {
634 char message[THD_MAX_NAME]; /* error message */
635
636 int num_ort_files; /* number of ort time series files */
637 int num_ideal_files; /* number of ideal time series files */
638 int is; /* time series file index */
639 int polort; /* degree of polynomial for baseline model */
640 int num_ortts; /* number of ort time series */
641 int num_idealts; /* number of ideal time series */
642 int q; /* number of parameters in the baseline model */
643 int p; /* number of parameters in the baseline model
644 plus number of ideals */
645
646
647 /*----- Initialize local variables -----*/
648 polort = option_data->polort;
649 num_ort_files = option_data->num_ort_files;
650 num_ideal_files = option_data->num_ideal_files;
651
652
653 /*----- Read the input fMRI measurement data -----*/
654 if (option_data->input1D_filename != NULL)
655 {
656 /*----- Read the input fMRI 1D time series -----*/
657 *fmri_data = read_one_time_series (option_data->input1D_filename,
658 fmri_length);
659 if (*fmri_data == NULL)
660 {
661 sprintf (message, "Unable to read time series file: %s",
662 option_data->input1D_filename);
663 FIM_error (message);
664 }
665 *dset_time = NULL;
666 }
667
668 else if (option_data->input_filename != NULL)
669 {
670 /*----- Read the input 3d+time dataset -----*/
671 *dset_time = THD_open_dataset (option_data->input_filename);
672 if (!ISVALID_3DIM_DATASET(*dset_time))
673 {
674 sprintf (message, "Unable to open data file: %s",
675 option_data->input_filename);
676 FIM_error (message);
677 }
678 DSET_load(*dset_time); CHECK_LOAD_ERROR(*dset_time);
679
680 if (option_data->mask_filename != NULL)
681 {
682 /*----- Read the input mask dataset -----*/
683 *mask_dset = THD_open_dataset (option_data->mask_filename);
684 if (!ISVALID_3DIM_DATASET(*mask_dset))
685 {
686 sprintf (message, "Unable to open mask file: %s",
687 option_data->mask_filename);
688 FIM_error (message);
689 }
690 DSET_load(*mask_dset); CHECK_LOAD_ERROR(*mask_dset);
691 }
692 }
693 else
694 FIM_error ("Must specify input measurement data");
695
696
697 /*----- Read the input ort time series files -----*/
698 for (is = 0; is < num_ort_files; is++)
699 {
700 ort_array[is] = read_time_series (option_data->ort_filename[is],
701 &(ort_list[is]));
702
703 if (ort_array[is] == NULL)
704 {
705 sprintf (message, "Unable to read ort time series file: %s",
706 option_data->ort_filename[is]);
707 FIM_error (message);
708 }
709 }
710
711
712 /*----- Read the input ideal time series files -----*/
713 for (is = 0; is < num_ideal_files; is++)
714 {
715 ideal_array[is] = read_time_series (option_data->ideal_filename[is],
716 &(ideal_list[is]));
717
718 if (ideal_array[is] == NULL)
719 {
720 sprintf (message, "Unable to read ideal time series file: %s",
721 option_data->ideal_filename[is]);
722 FIM_error (message);
723 }
724 }
725
726
727 /*----- Count number of ort and number of ideal time series -----*/
728 num_ortts = 0;
729 for (is = 0; is < num_ort_files; is++)
730 {
731 if (ort_list[is] == NULL)
732 num_ortts += ort_array[is]->ny;
733 else
734 num_ortts += ort_list[is][0];
735 }
736 q = polort + 1 + num_ortts;
737
738 num_idealts = 0;
739 for (is = 0; is < num_ideal_files; is++)
740 {
741 if (ideal_list[is] == NULL)
742 num_idealts += ideal_array[is]->ny;
743 else
744 num_idealts += ideal_list[is][0];
745 }
746 p = q + num_idealts;
747
748 option_data->num_ortts = num_ortts;
749 option_data->num_idealts = num_idealts;
750 option_data->q = q;
751 option_data->p = p;
752
753 }
754
755
756 /*---------------------------------------------------------------------------*/
757 /*
758 Routine to check whether one output file already exists.
759 */
760
check_one_output_file(THD_3dim_dataset * dset_time,char * filename)761 void check_one_output_file
762 (
763 THD_3dim_dataset * dset_time, /* input 3d+time data set */
764 char * filename /* name of output file */
765 )
766
767 {
768 char message[THD_MAX_NAME]; /* error message */
769 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
770 int ierror; /* number of errors in editing data */
771
772
773 /*----- make an empty copy of input dataset -----*/
774 new_dset = EDIT_empty_copy( dset_time ) ;
775
776
777 ierror = EDIT_dset_items( new_dset ,
778 ADN_prefix , filename ,
779 ADN_label1 , filename ,
780 ADN_self_name , filename ,
781 ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
782 GEN_FUNC_TYPE ,
783 ADN_none ) ;
784
785 if( ierror > 0 )
786 {
787 sprintf (message,
788 "*** %d errors in attempting to create output dataset!\n",
789 ierror);
790 FIM_error (message);
791 }
792
793 if( THD_is_file(new_dset->dblk->diskptr->header_name) )
794 {
795 sprintf (message,
796 "Output dataset file %s already exists "
797 " -- cannot continue!\a\n",
798 new_dset->dblk->diskptr->header_name);
799 FIM_error (message);
800 }
801
802 /*----- deallocate memory -----*/
803 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
804
805 }
806
807
808 /*---------------------------------------------------------------------------*/
809 /*
810 Routine to check whether output files already exist.
811 */
812
check_output_files(FIM_options * option_data,THD_3dim_dataset * dset_time)813 void check_output_files
814 (
815 FIM_options * option_data, /* fim program options */
816 THD_3dim_dataset * dset_time /* input 3d+time data set */
817 )
818
819 {
820
821 if ((option_data->bucket_filename != NULL)
822 && (option_data->input1D_filename == NULL))
823 check_one_output_file (dset_time, option_data->bucket_filename);
824
825 }
826
827
828 /*---------------------------------------------------------------------------*/
829 /*
830 Routine to check for valid inputs.
831 */
832
check_for_valid_inputs(FIM_options * option_data,THD_3dim_dataset * dset_time,THD_3dim_dataset * mask_dset,int fmri_length,MRI_IMAGE ** ort_array,MRI_IMAGE ** ideal_array)833 void check_for_valid_inputs
834 (
835 FIM_options * option_data, /* fim program options */
836 THD_3dim_dataset * dset_time, /* input 3d+time data set */
837 THD_3dim_dataset * mask_dset, /* input mask data set */
838 int fmri_length, /* length of input fMRI time series */
839 MRI_IMAGE ** ort_array, /* ort time series arrays */
840 MRI_IMAGE ** ideal_array /* ideal time series arrays */
841 )
842
843 {
844 char message[THD_MAX_NAME]; /* error message */
845 int is; /* ideal index */
846 int num_ort_files; /* number of ort time series files */
847 int num_ideal_files; /* number of ideal time series files */
848 int num_idealts; /* number of ideal time series */
849 int nt; /* number of images in input 3d+time dataset */
850 int NFirst; /* first image from input 3d+time dataset to use */
851 int NLast; /* last image from input 3d+time dataset to use */
852 int N; /* number of usable time points */
853
854
855 /*----- Initialize local variables -----*/
856 if (option_data->input1D_filename != NULL)
857 nt = fmri_length;
858 else
859 nt = DSET_NUM_TIMES (dset_time);
860
861 num_ort_files = option_data->num_ort_files;
862 num_ideal_files = option_data->num_ideal_files;
863 num_idealts = option_data->num_idealts;
864
865
866 NFirst = option_data->NFirst;
867
868 NLast = option_data->NLast;
869 if (NLast > nt-1) NLast = nt-1;
870 option_data->NLast = NLast;
871
872 N = NLast - NFirst + 1;
873 option_data->N = N;
874
875
876 /*----- Check number of ideal time series -----*/
877 if (num_idealts < 1) FIM_error ("No ideal time series?");
878 if (num_idealts < 2) option_data->output_type[FIM_BestIndex] = 0;
879
880
881 /*----- If mask is used, check for compatible dimensions -----*/
882 if (mask_dset != NULL)
883 {
884 if ( (DSET_NX(dset_time) != DSET_NX(mask_dset))
885 || (DSET_NY(dset_time) != DSET_NY(mask_dset))
886 || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) )
887 {
888 sprintf (message, "%s and %s have incompatible dimensions",
889 option_data->input_filename, option_data->mask_filename);
890 FIM_error (message);
891 }
892
893 if (DSET_NVALS(mask_dset) != 1 )
894 FIM_error ("Must specify 1 sub-brick from mask dataset");
895 }
896
897
898 /*----- Check lengths of ort time series -----*/
899 for (is = 0; is < num_ort_files; is++)
900 {
901 if (ort_array[is]->nx < NLast+1)
902 {
903 sprintf (message, "Input ort time series file %s is too short",
904 option_data->ort_filename[is]);
905 FIM_error (message);
906 }
907 }
908
909
910 /*----- Check lengths of ideal time series -----*/
911 for (is = 0; is < num_ideal_files; is++)
912 {
913 if (ideal_array[is]->nx < NLast+1)
914 {
915 sprintf (message, "Input ideal time series file %s is too short",
916 option_data->ideal_filename[is]);
917 FIM_error (message);
918 }
919 }
920
921
922 /*----- Check whether any of the output files already exist -----*/
923 if( THD_deathcon() ) check_output_files (option_data, dset_time);
924
925 }
926
927
928 /*---------------------------------------------------------------------------*/
929 /*
930 Allocate memory for output volumes.
931 */
932
allocate_memory(FIM_options * option_data,THD_3dim_dataset * dset,float *** fim_params_vol)933 void allocate_memory
934 (
935 FIM_options * option_data, /* fim program options */
936 THD_3dim_dataset * dset, /* input 3d+time data set */
937 float *** fim_params_vol /* array of volumes containing fim parameters */
938 )
939
940 {
941 int ip; /* parameter index */
942 int nxyz; /* total number of voxels */
943 int ixyz; /* voxel index */
944
945
946 /*----- Initialize local variables -----*/
947 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
948
949
950 /*----- Allocate memory space for fim parameters -----*/
951 *fim_params_vol = (float **) malloc (sizeof(float *) * MAX_OUTPUT_TYPE);
952 MTEST(*fim_params_vol);
953
954 for (ip = 0; ip < MAX_OUTPUT_TYPE; ip++)
955 {
956 if (option_data->output_type[ip])
957 {
958 (*fim_params_vol)[ip] = (float *) malloc (sizeof(float) * nxyz);
959 MTEST((*fim_params_vol)[ip]);
960 for (ixyz = 0; ixyz < nxyz; ixyz++)
961 {
962 (*fim_params_vol)[ip][ixyz] = 0.0;
963 }
964 }
965 else
966 (*fim_params_vol)[ip] = NULL;
967 }
968
969 }
970
971
972
973
974 /*---------------------------------------------------------------------------*/
975 /*
976 Perform all program initialization.
977 */
978
initialize_program(int argc,char ** argv,FIM_options ** option_data,THD_3dim_dataset ** dset_time,THD_3dim_dataset ** mask_dset,float ** fmri_data,int * fmri_length,MRI_IMAGE ** ort_array,int ** ort_list,MRI_IMAGE ** ideal_array,int ** ideal_list,float *** fim_params_vol)979 void initialize_program
980 (
981 int argc, /* number of input arguments */
982 char ** argv, /* array of input arguments */
983 FIM_options ** option_data, /* fim algorithm options */
984 THD_3dim_dataset ** dset_time, /* input 3d+time data set */
985 THD_3dim_dataset ** mask_dset, /* input mask data set */
986 float ** fmri_data, /* input fMRI time series data */
987 int * fmri_length, /* length of fMRI time series */
988 MRI_IMAGE ** ort_array, /* ort time series arrays */
989 int ** ort_list, /* list of ort time series */
990 MRI_IMAGE ** ideal_array, /* ideal time series arrays */
991 int ** ideal_list, /* list of ideal time series */
992 float *** fim_params_vol /* array of volumes containing fim parameters */
993 )
994
995 {
996 /*----- Allocate memory -----*/
997 *option_data = (FIM_options *) malloc (sizeof(FIM_options));
998
999
1000 /*----- Get command line inputs -----*/
1001 get_options (argc, argv, *option_data);
1002
1003
1004 /*----- Read input data -----*/
1005 read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length,
1006 ort_array, ort_list, ideal_array, ideal_list);
1007
1008
1009 /*----- Check for valid inputs -----*/
1010 check_for_valid_inputs (*option_data, *dset_time, *mask_dset,
1011 *fmri_length, ort_array, ideal_array);
1012
1013
1014 /*----- Allocate memory for output volumes -----*/
1015 if ((*option_data)->input1D_filename == NULL)
1016 allocate_memory (*option_data, *dset_time, fim_params_vol);
1017
1018 }
1019
1020
1021 /*---------------------------------------------------------------------------*/
1022 /*
1023 Get the time series for one voxel from the AFNI 3d+time data set.
1024 */
1025
extract_ts_array(THD_3dim_dataset * dset_time,int iv,float * ts_array)1026 void extract_ts_array
1027 (
1028 THD_3dim_dataset * dset_time, /* input 3d+time dataset */
1029 int iv, /* get time series for this voxel */
1030 float * ts_array /* time series data for voxel #iv */
1031 )
1032
1033 {
1034 MRI_IMAGE * im; /* intermediate float data */
1035 float * ar; /* pointer to float data */
1036 int ts_length; /* length of input 3d+time data set */
1037 int it; /* time index */
1038
1039
1040 /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/
1041 im = THD_extract_series (iv, dset_time, 0);
1042
1043
1044 /*----- Verify extraction -----*/
1045 if (im == NULL) FIM_error ("Unable to extract data from 3d+time dataset");
1046
1047
1048 /*----- Now extract time series from MRI_IMAGE -----*/
1049 ts_length = DSET_NUM_TIMES (dset_time);
1050 ar = MRI_FLOAT_PTR (im);
1051 for (it = 0; it < ts_length; it++)
1052 {
1053 ts_array[it] = ar[it];
1054 }
1055
1056
1057 /*----- Release memory -----*/
1058 mri_free (im); im = NULL;
1059
1060 }
1061
1062
1063 /*---------------------------------------------------------------------------*/
1064 /*
1065 Save results for this voxel.
1066 */
1067
save_voxel(int iv,float * fim_params,float ** fim_params_vol)1068 void save_voxel
1069 (
1070 int iv, /* current voxel index */
1071 float * fim_params, /* array of fim parameters */
1072 float ** fim_params_vol /* array of volumes of fim output parameters */
1073 )
1074
1075 {
1076 int ip; /* parameter index */
1077
1078
1079 /*----- Saved user requested fim parameters -----*/
1080 for (ip = 0; ip < MAX_OUTPUT_TYPE; ip++)
1081 {
1082 if (fim_params_vol[ip] != NULL)
1083 fim_params_vol[ip][iv] = fim_params[ip];
1084 }
1085
1086 }
1087
1088
1089 /*---------------------------------------------------------------------------*/
1090 /*
1091 Set fim threshold level.
1092 */
1093
set_fim_thr_level(int NFirst,float fim_thr,THD_3dim_dataset * dset)1094 float set_fim_thr_level
1095 (
1096 int NFirst, /* first usable data point */
1097 float fim_thr, /* fim threshold (as proportion of mean) */
1098 THD_3dim_dataset * dset /* input 3d+time data set */
1099 )
1100
1101 {
1102 int nt; /* number of time points */
1103 int nxyz; /* total number of voxels */
1104 int ixyz; /* voxel index */
1105 double sum; /* sum of voxel intensities */
1106 float fthr; /* fim threshold (as intensity level) */
1107 float * ts_array = NULL; /* time series data for individual voxel */
1108
1109
1110 /*----- Initialize local variables -----*/
1111 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
1112 nt = DSET_NUM_TIMES (dset);
1113
1114 ts_array = (float *) malloc (sizeof(float) * nt);
1115 MTEST (ts_array);
1116
1117 sum = 0.0; /* Ides March 2004 [rickr] */
1118 /*----- Sum values of all voxels at initial time point -----*/
1119 for (ixyz = 0; ixyz < nxyz; ixyz++)
1120 {
1121 extract_ts_array (dset, ixyz, ts_array);
1122 sum += fabs(ts_array[NFirst]);
1123 }
1124
1125
1126 /*----- Set fim intensity level threshold -----*/
1127 fthr = fim_thr * sum / nxyz;
1128
1129
1130 free (ts_array); ts_array = NULL;
1131
1132 return (fthr);
1133 }
1134
1135
1136 /*---------------------------------------------------------------------------*/
1137 /*
1138 Calculate the best cross correlation for each voxel.
1139 */
1140
calculate_results(FIM_options * option_data,THD_3dim_dataset * dset,THD_3dim_dataset * mask,float * fmri_data,int fmri_length,MRI_IMAGE ** ort_array,int ** ort_list,MRI_IMAGE ** ideal_array,int ** ideal_list,float ** fim_params_vol)1141 void calculate_results
1142 (
1143 FIM_options * option_data, /* fim program options */
1144 THD_3dim_dataset * dset, /* input 3d+time data set */
1145 THD_3dim_dataset * mask, /* input mask data set */
1146 float * fmri_data, /* input fMRI time series data */
1147 int fmri_length, /* length of fMRI time series */
1148 MRI_IMAGE ** ort_array, /* ort time series arrays */
1149 int ** ort_list, /* list of ort time series */
1150 MRI_IMAGE ** ideal_array, /* ideal time series arrays */
1151 int ** ideal_list, /* list of ideal time series */
1152 float ** fim_params_vol /* array of volumes of fim output parameters */
1153 )
1154
1155 {
1156 float * ts_array = NULL; /* array of measured data for one voxel */
1157 float mask_val[1]; /* value of mask at current voxel */
1158 float fthr=0.0; /* internal mask threshold level */
1159
1160 int q; /* number of parameters in the baseline model */
1161 int p; /* number of parameters in the baseline model
1162 plus number of ideals */
1163 int m; /* parameter index */
1164 int n; /* data point index */
1165
1166
1167 matrix xdata; /* independent variable matrix */
1168 matrix x_base; /* extracted X matrix for baseline model */
1169 matrix xtxinvxt_base; /* matrix: (1/(X'X))X' for baseline model */
1170 matrix x_ideal[MAX_FILES]; /* extracted X matrices for ideal models */
1171 matrix xtxinvxt_ideal[MAX_FILES];
1172 /* matrix: (1/(X'X))X' for ideal models */
1173 vector y; /* vector of measured data */
1174
1175 int ixyz; /* voxel index */
1176 int nxyz; /* number of voxels per dataset */
1177
1178 int nt; /* number of images in input 3d+time dataset */
1179 int NFirst; /* first image from input 3d+time dataset to use */
1180 int NLast; /* last image from input 3d+time dataset to use */
1181 int N; /* number of usable data points */
1182
1183 int num_ort_files; /* number of ort time series files */
1184 int num_ideal_files; /* number of ideal time series files */
1185 int polort; /* degree of polynomial ort */
1186 int num_ortts; /* number of ort time series */
1187 int num_idealts; /* number of ideal time series */
1188
1189 int i; /* data point index */
1190 int is; /* ideal index */
1191 int ilag; /* time lag index */
1192 float stddev; /* normalized parameter standard deviation */
1193 char * label; /* string containing stat. summary of results */
1194
1195 float * x_bot = NULL; /* minimum of stimulus time series */
1196 float * x_ave = NULL; /* average of stimulus time series */
1197 float * x_top = NULL; /* maximum of stimulus time series */
1198 int * good_list = NULL; /* list of good (usable) time points */
1199 float ** rarray = NULL; /* ranked arrays of ideal time series */
1200 float FimParams[MAX_OUTPUT_TYPE]; /* output fim parameters */
1201
1202
1203 /*----- Initialize matrices and vectors -----*/
1204 matrix_initialize (&xdata);
1205 matrix_initialize (&x_base);
1206 matrix_initialize (&xtxinvxt_base);
1207 for (is =0; is < MAX_FILES; is++)
1208 {
1209 matrix_initialize (&x_ideal[is]);
1210 matrix_initialize (&xtxinvxt_ideal[is]);
1211 }
1212 vector_initialize (&y);
1213
1214
1215 /*----- Initialize local variables -----*/
1216 if (option_data->input1D_filename != NULL)
1217 {
1218 nxyz = 1;
1219 nt = fmri_length;
1220 }
1221 else
1222 {
1223 nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
1224 nt = DSET_NUM_TIMES (dset);
1225 }
1226 NFirst = option_data->NFirst;
1227 NLast = option_data->NLast;
1228 N = option_data->N;
1229 p = option_data->p;
1230 q = option_data->q;
1231
1232 polort = option_data->polort;
1233 num_idealts = option_data->num_idealts;
1234 num_ort_files = option_data->num_ort_files;
1235 num_ideal_files = option_data->num_ideal_files;
1236 if (num_idealts > MAX_FILES ||
1237 num_ort_files > MAX_FILES ||
1238 num_ideal_files > MAX_FILES ) {
1239 /* ZSS: I used one ideal file with 30 regressors when MAX_FILE was 20
1240 and I spent most of the day chasing memory corruption errors.
1241 Tested with 20 regressors and all was well so now limit is 200 ZSS May 08 */
1242 fprintf(stderr,"Error: the number of ideal or ort regressors\n"
1243 "exceeds the hard coded limit of %d\n"
1244 "Contact authors if you need the limit extended.\n",
1245 MAX_FILES);
1246 exit(1);
1247 }
1248
1249 /*----- Allocate memory -----*/
1250 ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array);
1251 x_bot = (float *) malloc (sizeof(float) * p); MTEST (x_bot);
1252 x_ave = (float *) malloc (sizeof(float) * p); MTEST (x_ave);
1253 x_top = (float *) malloc (sizeof(float) * p); MTEST (x_top);
1254 good_list = (int *) malloc (sizeof(int) * N); MTEST (good_list);
1255 rarray = (float **) malloc (sizeof(float *) * num_idealts); MTEST (rarray);
1256
1257
1258 /*----- Initialize the independent variable matrix -----*/
1259 N = init_indep_var_matrix (q, p, NFirst, N, polort,
1260 num_ort_files, num_ideal_files,
1261 ort_array, ort_list, ideal_array, ideal_list,
1262 x_bot, x_ave, x_top, good_list, &xdata);
1263 option_data->N = N;
1264
1265
1266 /*----- Initialize fim threshold level -----*/
1267 if (option_data->input1D_filename == NULL)
1268 fthr = set_fim_thr_level (good_list[0]+NFirst, option_data->fim_thr, dset);
1269
1270
1271 /*----- Initialization for the regression analysis -----*/
1272 init_regression_analysis (q, p, N, num_idealts, xdata, &x_base,
1273 &xtxinvxt_base, x_ideal, xtxinvxt_ideal, rarray);
1274
1275
1276 vector_create (N, &y);
1277
1278
1279 /*----- Loop over all voxels -----*/
1280 for (ixyz = 0; ixyz < nxyz; ixyz++)
1281 {
1282 /*----- Apply mask? -----*/
1283 if (mask != NULL)
1284 {
1285 extract_ts_array (mask, ixyz, mask_val);
1286 if (mask_val[0] == 0.0) continue;
1287 }
1288
1289
1290 /*----- Extract Y-data for this voxel -----*/
1291 if (option_data->input1D_filename != NULL)
1292 {
1293 for (i = 0; i < N; i++)
1294 y.elts[i] = fmri_data[good_list[i]+NFirst];
1295 }
1296 else
1297 {
1298 extract_ts_array (dset, ixyz, ts_array);
1299 if (fabs(ts_array[good_list[0]+NFirst]) < fthr)
1300 continue;
1301 for (i = 0; i < N; i++)
1302 y.elts[i] = ts_array[good_list[i]+NFirst];
1303 }
1304
1305
1306 /*----- Perform the regression analysis for this voxel-----*/
1307 regression_analysis (N, q, num_idealts,
1308 x_base, xtxinvxt_base, x_ideal, xtxinvxt_ideal,
1309 y, x_bot, x_ave, x_top, rarray,
1310 option_data->output_type, FimParams);
1311
1312
1313 /*----- Save results for this voxel -----*/
1314 if (option_data->input1D_filename == NULL)
1315 save_voxel (ixyz, FimParams, fim_params_vol);
1316
1317
1318 /*----- Report results for this voxel -----*/
1319 if ( ((fabs(FimParams[FIM_Correlation]) > option_data->cdisp)
1320 && (option_data->cdisp >= 0.0))
1321 || (option_data->input1D_filename != NULL) )
1322 {
1323 printf ("\n\nResults for Voxel #%d: \n", ixyz);
1324 report_results (option_data->output_type, FimParams, &label);
1325 printf ("%s \n", label);
1326 }
1327
1328 } /*----- Loop over voxels -----*/
1329
1330
1331 /*----- Dispose of matrices and vectors -----*/
1332 vector_destroy (&y);
1333 for (is = 0; is < MAX_FILES; is++)
1334 {
1335 matrix_destroy (&xtxinvxt_ideal[is]);
1336 matrix_destroy (&x_ideal[is]);
1337 }
1338 matrix_destroy (&xtxinvxt_base);
1339 matrix_destroy (&x_base);
1340 matrix_destroy (&xdata);
1341
1342
1343 /*----- Deallocate memory -----*/
1344 free (ts_array); ts_array = NULL;
1345 free (x_bot); x_bot = NULL;
1346 free (x_ave); x_ave = NULL;
1347 free (x_top); x_top = NULL;
1348 free (good_list); good_list = NULL;
1349 for (is = 0; is < num_idealts; is++)
1350 {
1351 free (rarray[is]); rarray[is] = NULL;
1352 }
1353 free (rarray); rarray = NULL;
1354
1355 }
1356
1357 /*---------------------------------------------------------------------------*/
1358 /*
1359 Attach one sub-brick to output bucket data set.
1360 */
1361
attach_sub_brick(THD_3dim_dataset * new_dset,int ibrick,float * volume,int nxyz,int brick_type,char * brick_label,int nsam,int nfit,int nort,short ** bar,int do_scale)1362 void attach_sub_brick
1363 (
1364 THD_3dim_dataset * new_dset, /* output bucket dataset */
1365 int ibrick, /* sub-brick indices */
1366 float * volume, /* volume of floating point data */
1367 int nxyz, /* total number of voxels */
1368 int brick_type, /* indicates statistical type of sub-brick */
1369 char * brick_label, /* character string label for sub-brick */
1370 int nsam,
1371 int nfit,
1372 int nort, /* degrees of freedom */
1373 short ** bar, /* bar[ib] points to data for sub-brick #ib */
1374 int do_scale
1375 )
1376
1377 {
1378 const float EPSILON = 1.0e-10;
1379 float factor; /* factor is new scale factor for sub-brick #ib */
1380
1381
1382 /*----- allocate memory for output sub-brick -----*/
1383 bar[ibrick] = (short *) malloc (sizeof(short) * nxyz);
1384 MTEST (bar[ibrick]);
1385 if (do_scale) {
1386 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
1387 MRI_short, bar[ibrick]);
1388
1389 if (factor < EPSILON) factor = 0.0;
1390 else factor = 1.0 / factor;
1391
1392 EDIT_misfit_report( DSET_FILECODE(new_dset) , ibrick ,
1393 nxyz , factor , bar[ibrick] , volume ) ;
1394 } else { /* ZSS Nov 2011 */
1395 EDIT_coerce_type(nxyz, MRI_float, volume,
1396 MRI_short, bar[ibrick]);
1397 factor = 0.0;
1398 }
1399 /*----- edit the sub-brick -----*/
1400 EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
1401 EDIT_BRICK_FACTOR (new_dset, ibrick, factor);
1402
1403 if (brick_type == FUNC_COR_TYPE)
1404 EDIT_BRICK_TO_FICO (new_dset, ibrick, nsam, nfit, nort);
1405
1406
1407 /*----- attach bar[ib] to be sub-brick #ibrick -----*/
1408 EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]);
1409
1410 }
1411
1412 /*---------------------------------------------------------------------------*/
1413 /*
1414 Routine to write one bucket data set.
1415 */
1416
write_bucket_data(int argc,char ** argv,FIM_options * option_data,float ** fim_params_vol)1417 void write_bucket_data
1418 (
1419 int argc, /* number of input arguments */
1420 char ** argv, /* array of input arguments */
1421 FIM_options * option_data, /* fim program options */
1422 float ** fim_params_vol /* array of volumes of fim output parameters */
1423 )
1424
1425 {
1426 THD_3dim_dataset * old_dset = NULL; /* prototype dataset */
1427 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */
1428 char output_prefix[THD_MAX_NAME]; /* prefix name for bucket dataset */
1429 char output_session[THD_MAX_NAME]; /* directory for bucket dataset */
1430 int nbricks; /* number of sub-bricks in bucket dataset */
1431 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */
1432
1433 int brick_type; /* indicates statistical type of sub-brick */
1434 int brick_coef; /* regression coefficient index for sub-brick */
1435 char brick_label[THD_MAX_NAME]; /* character string label for sub-brick */
1436
1437 int ierror; /* number of errors in editing data */
1438 float * volume; /* volume of floating point data */
1439
1440 int N; /* number of usable data points */
1441 int q; /* number of parameters in the ideal model */
1442 int num_idealts; /* number of ideal time series */
1443 int ip; /* parameter index */
1444 int nxyz; /* total number of voxels */
1445 int ibrick; /* sub-brick index */
1446 int nsam;
1447 int nfit;
1448 int nort; /* degrees of freedom */
1449 char label[THD_MAX_NAME]; /* general label for sub-bricks */
1450
1451
1452 /*----- read prototype dataset -----*/
1453 old_dset = THD_open_dataset (option_data->input_filename);
1454 /* ZSS May 08, changed from THD_open_one_dataset */
1455
1456
1457 /*----- Initialize local variables -----*/
1458 nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;
1459 num_idealts = option_data->num_idealts;
1460 q = option_data->q;
1461 N = option_data->N;
1462
1463
1464 /*----- Calculate number of sub-bricks in the bucket -----*/
1465 nbricks = 0;
1466 for (ip = 0; ip < MAX_OUTPUT_TYPE; ip++)
1467 if (option_data->output_type[ip]) nbricks++;
1468
1469
1470 strcpy (output_prefix, option_data->bucket_filename);
1471 strcpy (output_session, "./");
1472
1473
1474 /*----- allocate memory -----*/
1475 bar = (short **) malloc (sizeof(short *) * nbricks);
1476 MTEST (bar);
1477
1478
1479 /*-- make an empty copy of prototype dataset, for eventual output --*/
1480 new_dset = EDIT_empty_copy (old_dset);
1481
1482
1483 /*----- Record history of dataset -----*/
1484 tross_Copy_History( old_dset , new_dset ) ;
1485 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
1486 sprintf (label, "Output prefix: %s", output_prefix);
1487 tross_Append_History ( new_dset, label);
1488
1489
1490 /*----- delete prototype dataset -----*/
1491 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
1492
1493
1494 /*----- Modify some structural properties. Note that the nbricks
1495 just make empty sub-bricks, without any data attached. -----*/
1496 ierror = EDIT_dset_items (new_dset,
1497 ADN_prefix, output_prefix,
1498 ADN_directory_name, output_session,
1499 ADN_type, HEAD_FUNC_TYPE,
1500 ADN_func_type, FUNC_BUCK_TYPE,
1501 ADN_datum_all, MRI_short ,
1502 ADN_ntt, 0, /* no time */
1503 ADN_nvals, nbricks,
1504 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
1505 ADN_none ) ;
1506
1507 if( ierror > 0 )
1508 {
1509 fprintf(stderr,
1510 "*** %d errors in attempting to create bucket dataset!\n",
1511 ierror);
1512 exit(1);
1513 }
1514
1515 if (!THD_ok_overwrite() && THD_is_file(DSET_HEADNAME(new_dset)))
1516 {
1517 fprintf(stderr,
1518 "*** Output dataset file %s already exists--cannot continue!\n",
1519 DSET_HEADNAME(new_dset));
1520 exit(1);
1521 }
1522
1523
1524 /*----- Attach individual sub-bricks to the bucket dataset -----*/
1525 ibrick = 0;
1526 for (ip = 0; ip < MAX_OUTPUT_TYPE; ip++)
1527 {
1528 if (option_data->output_type[ip] == 0) continue;
1529
1530 strcpy (brick_label, OUTPUT_TYPE_name[ip]);
1531
1532 if (ip == FIM_Correlation)
1533 {
1534 brick_type = FUNC_COR_TYPE;
1535 nsam = N; nort = q;
1536 if (num_idealts > 1) nfit = 2;
1537 else nfit = 1;
1538 }
1539 else if (ip == FIM_SpearmanCC)
1540 {
1541 #if 0
1542 brick_type = FUNC_THR_TYPE;
1543 #else
1544 brick_type = FUNC_COR_TYPE;
1545 #endif
1546 nsam = N; nort = q;
1547 if (num_idealts > 1) nfit = 2;
1548 else nfit = 1;
1549 }
1550 else if (ip == FIM_QuadrantCC)
1551 {
1552 #if 0
1553 brick_type = FUNC_THR_TYPE;
1554 #else
1555 brick_type = FUNC_COR_TYPE;
1556 #endif
1557 nsam = N; nort = q;
1558 if (num_idealts > 1) nfit = 2;
1559 else nfit = 1;
1560 }
1561 else
1562 {
1563 brick_type = FUNC_FIM_TYPE;
1564 nsam = 0; nfit = 0; nort = 0;
1565 }
1566
1567 volume = fim_params_vol[ip];
1568 attach_sub_brick (new_dset, ibrick, volume, nxyz,
1569 brick_type, brick_label, nsam, nfit, nort, bar,
1570 ip == FIM_BestIndex ? 0 : 1);
1571
1572 ibrick++;
1573 }
1574
1575
1576 /*----- write bucket data set -----*/
1577 THD_load_statistics (new_dset);
1578 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
1579 fprintf(stderr,"Wrote bucket dataset ");
1580 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
1581
1582
1583 /*----- deallocate memory -----*/
1584 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
1585
1586 }
1587
1588
1589 /*---------------------------------------------------------------------------*/
1590 /*
1591 Write out the user requested output files.
1592 */
1593
output_results(int argc,char ** argv,FIM_options * option_data,float ** fim_params_vol)1594 void output_results
1595 (
1596 int argc, /* number of input arguments */
1597 char ** argv, /* array of input arguments */
1598 FIM_options * option_data, /* fim algorithm options */
1599 float ** fim_params_vol /* array of volumes of fim output parameters */
1600 )
1601
1602 {
1603 int q; /* number of parameters in baseline model */
1604 int num_idealts; /* number of ideal time series */
1605 int ib; /* sub-brick index */
1606 int is; /* ideal index */
1607 int ts_length; /* length of impulse reponse function */
1608 int N; /* number of usable data points */
1609
1610
1611 /*----- Initialize local variables -----*/
1612 q = option_data->polort + 1;
1613 num_idealts = option_data->num_idealts;
1614 N = option_data->N;
1615
1616
1617 /*----- Write the bucket dataset -----*/
1618 if (option_data->bucket_filename != NULL)
1619 write_bucket_data (argc, argv, option_data, fim_params_vol);
1620
1621 }
1622
1623
1624 /*---------------------------------------------------------------------------*/
1625
terminate_program(FIM_options ** option_data,MRI_IMAGE ** ort_array,int ** ort_list,MRI_IMAGE ** ideal_array,int ** ideal_list,float *** fim_params_vol)1626 void terminate_program
1627 (
1628 FIM_options ** option_data, /* fim program options */
1629 MRI_IMAGE ** ort_array, /* ort time series arrays */
1630 int ** ort_list, /* list of ort time series */
1631 MRI_IMAGE ** ideal_array, /* ideal time series arrays */
1632 int ** ideal_list, /* list of ideal time series */
1633 float *** fim_params_vol /* array of volumes of fim output parameters */
1634 )
1635
1636 {
1637 int num_idealts; /* number of ideal time series */
1638 int ip; /* parameter index */
1639 int is; /* ideal index */
1640
1641
1642 /*----- Initialize local variables -----*/
1643 num_idealts = (*option_data)->num_idealts;
1644
1645
1646 /*----- Deallocate memory for option data -----*/
1647 free (*option_data); *option_data = NULL;
1648
1649
1650 /*----- Deallocate memory for ideal time series -----*/
1651 /*
1652 for (is = 0; is < num_idealts; is++)
1653 { free (ideal[is]); ideal[is] = NULL; }
1654 */
1655
1656
1657 /*----- Deallocate space for volume data -----*/
1658 if (*fim_params_vol != NULL)
1659 {
1660 for (ip = 0; ip < MAX_OUTPUT_TYPE; ip++)
1661 {
1662 if ((*fim_params_vol)[ip] != NULL)
1663 { free ((*fim_params_vol)[ip]); (*fim_params_vol)[ip] = NULL; }
1664 }
1665
1666 free (*fim_params_vol); *fim_params_vol = NULL;
1667 }
1668
1669 }
1670
1671
1672 /*---------------------------------------------------------------------------*/
1673
main(int argc,char ** argv)1674 int main
1675 (
1676 int argc, /* number of input arguments */
1677 char ** argv /* array of input arguments */
1678 )
1679
1680 {
1681 FIM_options * option_data; /* fim algorithm options */
1682 THD_3dim_dataset * dset_time = NULL; /* input 3d+time data set */
1683 THD_3dim_dataset * mask_dset = NULL; /* input mask data set */
1684 float * fmri_data = NULL; /* input fMRI time series data */
1685 int fmri_length; /* length of fMRI time series */
1686 MRI_IMAGE * ort_array[MAX_FILES]; /* ideal time series arrays */
1687 int * ort_list[MAX_FILES]; /* list of ideal time series */
1688 MRI_IMAGE * ideal_array[MAX_FILES]; /* ideal time series arrays */
1689 int * ideal_list[MAX_FILES]; /* list of ideal time series */
1690
1691 float ** fim_params_vol = NULL;
1692 /* array of volumes of fim output parameters */
1693
1694
1695 /*----- Identify software -----*/
1696 #if 0
1697 printf ("\n\n");
1698 printf ("Program: %s \n", PROGRAM_NAME);
1699 printf ("Author: %s \n", PROGRAM_AUTHOR);
1700 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
1701 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
1702 printf ("\n");
1703 #endif
1704
1705 PRINT_VERSION("3dfim+") ; AUTHOR(PROGRAM_AUTHOR) ;
1706 mainENTRY("3dfim+ main") ; machdep() ;
1707
1708 WARNING_message("This program (3dfim+) is very old, may not be useful, and will not be maintained.") ;
1709
1710 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
1711
1712 { int new_argc ; char ** new_argv ;
1713 addto_args( argc , argv , &new_argc , &new_argv ) ;
1714 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
1715 }
1716
1717 #ifdef ALLOW_MCW_MALLOC
1718 enable_mcw_malloc();
1719 #endif
1720
1721 /*----- Program initialization -----*/
1722 initialize_program (argc, argv, &option_data, &dset_time, &mask_dset,
1723 &fmri_data, &fmri_length,
1724 ort_array, ort_list, ideal_array, ideal_list,
1725 &fim_params_vol);
1726
1727
1728 /*----- Perform fim analysis -----*/
1729 calculate_results (option_data, dset_time, mask_dset,
1730 fmri_data, fmri_length,
1731 ort_array, ort_list, ideal_array, ideal_list,
1732 fim_params_vol);
1733
1734 /*----- Deallocate memory for input datasets -----*/
1735 if (dset_time != NULL)
1736 { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; }
1737 if (mask_dset != NULL)
1738 { THD_delete_3dim_dataset (mask_dset, False); mask_dset = NULL; }
1739
1740
1741 /*----- Write requested output files -----*/
1742 if (option_data->input1D_filename == NULL)
1743 output_results (argc, argv, option_data, fim_params_vol);
1744
1745
1746 /*----- Terminate program -----*/
1747 terminate_program (&option_data,
1748 ort_array, ort_list, ideal_array, ideal_list,
1749 &fim_params_vol);
1750
1751 exit(0);
1752 }
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762