00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   Program to calculate the cross-correlation of an ideal reference waveform  
00010   with the measured FMRI time series for each voxel.                          
00011   This is the stand-alone (batch command) version of the afni fim+ function.
00012 
00013   File:    3dfim+.c
00014   Author:  B. Douglas Ward
00015   Date:    28 April 2000
00016 
00017 
00018   Mod:     Use output_type array in regression_analysis routine to avoid
00019            some unnecessary calculations.
00020   Date:    18 May 2000
00021 
00022   Mod:     Added call to AFNI_logger.
00023   Date:    15 August 2001
00024 
00025   Mod:     Add FICO-ness of sub-bricks for Spearman and Quadrant correlation.
00026   Date     29 Oct 2004 - RWCox
00027 
00028 */
00029 
00030 /*---------------------------------------------------------------------------*/
00031 
00032 #define PROGRAM_NAME "3dfim+"                        /* name of this program */
00033 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00034 #define PROGRAM_INITIAL "28 April 2000"   /* date of initial program release */
00035 #define PROGRAM_LATEST  "29 October 2004"  /* date of latest program revision */
00036 
00037 /*---------------------------------------------------------------------------*/
00038 
00039 #define MAX_FILES 20                        /* maximum number of ideal files */
00040                                             /* = maximum number of ort files */
00041 #define RA_error FIM_error
00042 
00043 /*---------------------------------------------------------------------------*/
00044 
00045 #include "mrilib.h"
00046 #include "matrix.h"
00047 
00048 #include "fim+.c"
00049 
00050 /*---------------------------------------------------------------------------*/
00051 
00052 typedef struct FIM_options
00053 { 
00054   int NFirst;              /* first image from input 3d+time dataset to use */
00055   int NLast;               /* last image from input 3d+time dataset to use */
00056   int N;                   /* number of usable data points from input data */
00057   int polort;              /* degree of polynomial for baseline model */
00058   int num_ortts;           /* number of ort time series */
00059   int num_idealts;         /* number of ideal time series */
00060   int q;                   /* number of parameters in the baseline model */
00061   int p;                   /* number of parameters in the baseline model 
00062                               plus number of ideals */
00063 
00064   float fim_thr;           /* threshold for internal fim mask */
00065   float cdisp;             /* minimum correlation coefficient for display */ 
00066 
00067   char * input_filename;   /* input 3d+time dataset filename */
00068   char * mask_filename;    /* input mask dataset filename */
00069   char * input1D_filename; /* input fMRI measurement time series */
00070 
00071   int num_ort_files;                  /* number of ort files */
00072   char * ort_filename[MAX_FILES];     /* input ort time series file names */
00073   int num_ideal_files;                /* number of ideal files */
00074   char * ideal_filename[MAX_FILES];   /* input ideal time series file names */
00075   char * bucket_filename;             /* output bucket dataset file name */
00076 
00077   int output_type[MAX_OUTPUT_TYPE];   /* output type options */
00078 
00079 } FIM_options;
00080 
00081 
00082 /*---------------------------------------------------------------------------*/
00083 /*
00084   Get the time series for one voxel from the AFNI 3d+time data set.
00085 */
00086 
00087 void extract_ts_array 
00088 (
00089   THD_3dim_dataset * dset_time,      /* input 3d+time dataset */
00090   int iv,                            /* get time series for this voxel */
00091   float * ts_array                   /* time series data for voxel #iv */
00092 );
00093 
00094 /*---------------------------------------------------------------------------*/
00095 /*
00096    Print error message and stop.
00097 */
00098 
00099 void FIM_error (char * message)
00100 {
00101   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00102   exit(1);
00103 }
00104 
00105 
00106 /*---------------------------------------------------------------------------*/
00107 /*
00108   Routine to display 3dfim+ help menu.
00109 */
00110 
00111 void display_help_menu()
00112 {
00113   printf (
00114 "Program to calculate the cross-correlation of an ideal reference waveform  \n"
00115 "with the measured FMRI time series for each voxel.                         \n"
00116     "                                                                       \n"
00117     "Usage:                                                                 \n"
00118     "3dfim+                                                                 \n"
00119     "-input fname       fname = filename of input 3d+time dataset           \n"
00120     "[-input1D dname]   dname = filename of single (fMRI) .1D time series   \n"
00121     "[-mask mname]      mname = filename of 3d mask dataset                 \n"
00122     "[-nfirst fnum]     fnum = number of first dataset image to use in      \n"
00123     "                     the cross-correlation procedure. (default = 0)    \n"
00124     "[-nlast  lnum]     lnum = number of last dataset image to use in       \n"
00125     "                     the cross-correlation procedure. (default = last) \n"
00126     "[-polort pnum]     pnum = degree of polynomial corresponding to the    \n"
00127     "                     baseline model  (pnum = 0, 1, etc.)               \n"
00128     "                     (default: pnum = 1)                               \n"
00129     "[-fim_thr p]       p = fim internal mask threshold value (0 <= p <= 1) \n"
00130     "                     (default: p = 0.0999)                             \n"
00131     "[-cdisp cval]      Write (to screen) results for those voxels          \n"
00132     "                     whose correlation stat. > cval  (0 <= cval <= 1)  \n"
00133     "                     (default: disabled)                               \n"
00134     "[-ort_file sname]  sname = input ort time series file name             \n"
00135     "-ideal_file rname  rname = input ideal time series file name           \n"
00136     "                                                                       \n"
00137     "            Note:  The -ort_file and -ideal_file commands may be used  \n"
00138     "                   more than once.                                     \n"
00139     "            Note:  If files sname or rname contain multiple columns,   \n"
00140     "                   then ALL columns will be used as ort or ideal       \n"
00141     "                   time series.  However, individual columns or        \n"
00142     "                   a subset of columns may be selected using a file    \n"
00143     "                   name specification like 'fred.1D[0,3,5]', which     \n"
00144     "                   indicates that only columns #0, #3, and #5 will     \n"
00145     "                   be used for input.                                  \n"
00146     );
00147 
00148   printf("\n"
00149     "[-out param]       Flag to output the specified parameter, where       \n"
00150     "                   the string 'param' may be any one of the following: \n"
00151     "                                                                       \n"
00152     "%12s       L.S. fit coefficient for Best Ideal                \n"
00153     "%12s       Index number for Best Ideal                        \n"
00154     "%12s       P-P amplitude of signal response / Baseline        \n"
00155     "%12s       Average of baseline model response                 \n"
00156     "%12s       Best Ideal product-moment correlation coefficient  \n"
00157     "%12s       P-P amplitude of signal response / Average         \n"
00158     "%12s       Baseline + average of signal response              \n"
00159     "%12s       P-P amplitude of signal response / Topline         \n"
00160     "%12s       Baseline + P-P amplitude of signal response        \n"
00161     "%12s       Std. Dev. of residuals from best fit               \n"
00162     "%9sAll       This specifies all of the above parameters       \n"
00163     "%12s       Spearman correlation coefficient                   \n"
00164     "%12s       Quadrant correlation coefficient                   \n"
00165     "                                                                       \n"
00166     "            Note:  Multiple '-out' commands may be used.               \n"
00167     "            Note:  If a parameter name contains imbedded spaces, the   \n"
00168     "                   entire parameter name must be enclosed by quotes,   \n"
00169     "                   e.g.,  -out '%8s'                                   \n"
00170     "                                                                       \n"
00171     "[-bucket bprefix]  Create one AFNI 'bucket' dataset containing the     \n"
00172     "                   parameters of interest, as specified by the above   \n"
00173     "                   '-out' commands.                                    \n"
00174     "                   The output 'bucket' dataset is written to a file    \n"
00175     "                   with the prefix name bprefix.                       \n"
00176     ,
00177     OUTPUT_TYPE_name[FIM_FitCoef],
00178     OUTPUT_TYPE_name[FIM_BestIndex],
00179     OUTPUT_TYPE_name[FIM_PrcntChange],
00180     OUTPUT_TYPE_name[FIM_Baseline],
00181     OUTPUT_TYPE_name[FIM_Correlation],
00182     OUTPUT_TYPE_name[FIM_PrcntFromAve],
00183     OUTPUT_TYPE_name[FIM_Average],
00184     OUTPUT_TYPE_name[FIM_PrcntFromTop],
00185     OUTPUT_TYPE_name[FIM_Topline],
00186     OUTPUT_TYPE_name[FIM_SigmaResid],
00187     "",
00188     OUTPUT_TYPE_name[FIM_SpearmanCC],
00189     OUTPUT_TYPE_name[FIM_QuadrantCC],
00190     OUTPUT_TYPE_name[FIM_FitCoef]
00191 );
00192   
00193   exit(0);
00194 }
00195 
00196 
00197 /*---------------------------------------------------------------------------*/
00198 /*
00199   Routine to initialize the input options.
00200 */
00201  
00202 void initialize_options 
00203 (
00204   FIM_options * option_data    /* fim program options */
00205 )
00206  
00207 {
00208   int is;                     /* index */
00209 
00210 
00211   /*----- Initialize default values -----*/
00212   option_data->NFirst = 0;
00213   option_data->NLast  = 32767;
00214   option_data->N      = 0;
00215   option_data->polort = 1;
00216   option_data->num_ortts = 0;
00217   option_data->num_idealts = 0;
00218   option_data->q = 0;
00219   option_data->p = 0;
00220 
00221   option_data->fim_thr = 0.0999;
00222   option_data->cdisp = -1.0;
00223 
00224   option_data->num_ort_files = 0;
00225   option_data->num_ideal_files = 0;
00226 
00227 
00228   /*----- Initialize output flags -----*/
00229   for (is = 0;  is < MAX_OUTPUT_TYPE;  is++)
00230     option_data->output_type[is] = 0;
00231 
00232 
00233   /*----- Initialize file names -----*/
00234   option_data->input_filename = NULL;
00235   option_data->mask_filename = NULL;  
00236   option_data->input1D_filename = NULL;
00237   option_data->bucket_filename = NULL;
00238 
00239   for (is = 0;  is < MAX_FILES;  is++)
00240     {  
00241       option_data->ort_filename[is] = NULL;
00242       option_data->ideal_filename[is] = NULL;
00243     }
00244 
00245 }
00246 
00247 
00248 /*---------------------------------------------------------------------------*/
00249 /*
00250   Routine to get user specified input options.
00251 */
00252 
00253 void get_options
00254 (
00255   int argc,                        /* number of input arguments */
00256   char ** argv,                    /* array of input arguments */ 
00257   FIM_options * option_data        /* fim program options */
00258 )
00259 
00260 {
00261   int nopt = 1;                     /* input option argument counter */
00262   int ival, index;                  /* integer input */
00263   float fval;                       /* float input */
00264   char message[THD_MAX_NAME];       /* error message */
00265   int k;                            /* ideal time series index */
00266 
00267 
00268   /*----- does user request help menu? -----*/
00269   if (argc < 2 || strcmp(argv[1], "-help") == 0)  display_help_menu();  
00270 
00271   
00272   /*----- add to program log -----*/
00273   AFNI_logger (PROGRAM_NAME,argc,argv); 
00274 
00275   
00276   /*----- initialize the input options -----*/
00277   initialize_options (option_data); 
00278 
00279   
00280   /*----- main loop over input options -----*/
00281   while (nopt < argc )
00282     {
00283 
00284       /*-----   -input filename   -----*/
00285       if (strcmp(argv[nopt], "-input") == 0)
00286         {
00287           nopt++;
00288           if (nopt >= argc)  FIM_error ("need argument after -input ");
00289           option_data->input_filename = malloc (sizeof(char)*THD_MAX_NAME);
00290           MTEST (option_data->input_filename);
00291           strcpy (option_data->input_filename, argv[nopt]);
00292           nopt++;
00293           continue;
00294         }
00295       
00296 
00297       /*-----   -mask filename   -----*/
00298       if (strcmp(argv[nopt], "-mask") == 0)
00299         {
00300           nopt++;
00301           if (nopt >= argc)  FIM_error ("need argument after -mask ");
00302           option_data->mask_filename = malloc (sizeof(char)*THD_MAX_NAME);
00303           MTEST (option_data->mask_filename);
00304           strcpy (option_data->mask_filename, argv[nopt]);
00305           nopt++;
00306           continue;
00307         }
00308       
00309 
00310       /*-----   -input1D filename   -----*/
00311       if (strcmp(argv[nopt], "-input1D") == 0)
00312         {
00313           nopt++;
00314           if (nopt >= argc)  FIM_error ("need argument after -input1D ");
00315           option_data->input1D_filename = 
00316             malloc (sizeof(char)*THD_MAX_NAME);
00317           MTEST (option_data->input1D_filename);
00318           strcpy (option_data->input1D_filename, argv[nopt]);
00319           nopt++;
00320           continue;
00321         }
00322       
00323 
00324       /*-----   -nfirst num  -----*/
00325       if (strcmp(argv[nopt], "-nfirst") == 0)
00326         {
00327           nopt++;
00328           if (nopt >= argc)  FIM_error ("need argument after -nfirst ");
00329           sscanf (argv[nopt], "%d", &ival);
00330           if (ival < 0)
00331             FIM_error ("illegal argument after -nfirst ");
00332           option_data->NFirst = ival;
00333           nopt++;
00334           continue;
00335         }
00336 
00337 
00338       /*-----   -nlast num  -----*/
00339       if (strcmp(argv[nopt], "-nlast") == 0)
00340         {
00341           nopt++;
00342           if (nopt >= argc)  FIM_error ("need argument after -nlast ");
00343           sscanf (argv[nopt], "%d", &ival);
00344           if (ival < 0)
00345             FIM_error ("illegal argument after -nlast ");
00346           option_data->NLast = ival;
00347           nopt++;
00348           continue;
00349         }
00350 
00351 
00352       /*-----   -polort num  -----*/
00353       if (strcmp(argv[nopt], "-polort") == 0)
00354         {
00355           nopt++;
00356           if (nopt >= argc)  FIM_error ("need argument after -polort ");
00357           sscanf (argv[nopt], "%d", &ival);
00358 
00359 #undef PMAX
00360 #ifdef USE_LEGENDRE
00361 # define PMAX 19
00362 #else
00363 # define PMAX  2
00364 #endif
00365           if ((ival < 0) || (ival > PMAX))
00366             FIM_error ("illegal argument after -polort ");
00367 
00368 #ifdef USE_LEGENDRE
00369      if( ival > 2 )
00370        fprintf(stderr,
00371             "** WARNING: -polort > 2 is a new untested option: 29 Mar 2005\n") ;
00372 #endif
00373 
00374           option_data->polort = ival;
00375           nopt++;
00376           continue;
00377         }
00378 
00379       
00380       /*-----   -fim_thr r  -----*/
00381       if (strcmp(argv[nopt], "-fim_thr") == 0)
00382         {
00383           nopt++;
00384           if (nopt >= argc)  FIM_error ("need argument after -fim_thr ");
00385           sscanf (argv[nopt], "%f", &fval); 
00386           if ((fval < 0.0) || (fval > 1.0))
00387             FIM_error ("illegal argument after -fim_thr ");
00388           option_data->fim_thr = fval;
00389           nopt++;
00390           continue;
00391         }
00392       
00393 
00394       /*-----   -cdisp cval   -----*/
00395       if (strcmp(argv[nopt], "-cdisp") == 0)
00396         {
00397           nopt++;
00398           if (nopt >= argc)  FIM_error ("need argument after -cdisp ");
00399           sscanf (argv[nopt], "%f", &fval); 
00400           if ((fval < 0.0) || (fval > 1.0))
00401             FIM_error ("illegal argument after -cdisp ");
00402           option_data->cdisp = fval;
00403           nopt++;
00404           continue;
00405         }
00406       
00407 
00408       /*-----   -ort_file sname   -----*/
00409       if (strcmp(argv[nopt], "-ort_file") == 0)
00410         {
00411           nopt++;
00412           if (nopt >= argc)  FIM_error ("need argument after -ort_file");
00413 
00414           k = option_data->num_ort_files;
00415           if (k+1 > MAX_FILES)
00416             {
00417               sprintf (message, "Too many ( > %d ) ort time series files. ", 
00418                        MAX_FILES);
00419               FIM_error (message);
00420             }
00421 
00422           option_data->ort_filename[k] 
00423             = malloc (sizeof(char)*THD_MAX_NAME);
00424           MTEST (option_data->ort_filename[k]);
00425           strcpy (option_data->ort_filename[k], argv[nopt]);
00426           option_data->num_ort_files++;
00427           nopt++;
00428           continue;
00429         }
00430       
00431 
00432       /*-----   -ideal_file rname   -----*/
00433       if (strcmp(argv[nopt], "-ideal_file") == 0)
00434         {
00435           nopt++;
00436           if (nopt >= argc)  FIM_error ("need argument after -ideal_file");
00437 
00438           k = option_data->num_ideal_files;
00439           if (k+1 > MAX_FILES)
00440             {
00441               sprintf (message, "Too many ( > %d ) ideal time series files. ", 
00442                        MAX_FILES);
00443               FIM_error (message);
00444             }
00445 
00446           option_data->ideal_filename[k] 
00447             = malloc (sizeof(char)*THD_MAX_NAME);
00448           MTEST (option_data->ideal_filename[k]);
00449           strcpy (option_data->ideal_filename[k], argv[nopt]);
00450           option_data->num_ideal_files++;
00451           nopt++;
00452           continue;
00453         }
00454       
00455 
00456       /*-----   -out option   -----*/
00457       if (strcmp(argv[nopt], "-out") == 0)
00458         {
00459           nopt++;
00460           if (nopt >= argc)  FIM_error ("need argument after -out ");
00461           for (ival = 0;  ival < MAX_OUTPUT_TYPE;  ival++)
00462             if (strcmp(argv[nopt], OUTPUT_TYPE_name[ival]) == 0)  break;
00463           if (ival < MAX_OUTPUT_TYPE)
00464             {
00465               option_data->output_type[ival] = 1;
00466               nopt++;
00467               continue;
00468             }
00469           else if (strcmp(argv[nopt], "All") == 0)
00470             {
00471               for (ival = 0;  ival < MAX_OUTPUT_TYPE-2;  ival++)
00472                 option_data->output_type[ival] = 1;
00473               nopt++;
00474               continue;
00475             }
00476           else
00477             {
00478               sprintf(message,"Unrecognized output type: %s\n", argv[nopt]);
00479               FIM_error (message);
00480             }
00481         }
00482       
00483 
00484       /*-----   -bucket filename   -----*/
00485       if (strcmp(argv[nopt], "-bucket") == 0)
00486         {
00487           nopt++;
00488           if (nopt >= argc)  FIM_error ("need file prefixname after -bucket ");
00489           option_data->bucket_filename = malloc (sizeof(char)*THD_MAX_NAME);
00490           MTEST (option_data->bucket_filename);
00491           strcpy (option_data->bucket_filename, argv[nopt]);
00492           nopt++;
00493           continue;
00494         }
00495       
00496 
00497       /*----- unknown command -----*/
00498       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00499       FIM_error (message);
00500       
00501     }
00502 
00503 }
00504 
00505 
00506 /*---------------------------------------------------------------------------*/
00507 /*
00508   Read one time series from specified file name.  This file name may have
00509   a column selector attached.
00510 */
00511 
00512 float * read_one_time_series 
00513 (
00514   char * ts_filename,          /* time series file name (plus column index) */
00515   int * ts_length              /* output value for time series length */
00516 )
00517 
00518 {
00519   char message[THD_MAX_NAME];    /* error message */
00520   char * cpt;                    /* pointer to column suffix */
00521   char subv[THD_MAX_NAME];       /* string containing column index */
00522   MRI_IMAGE * im, * flim;  /* pointers to image structures 
00523                               -- used to read 1D ASCII */
00524   float * far;             /* pointer to MRI_IMAGE floating point data */
00525   int nx;                  /* number of time points in time series */
00526   int ny;                  /* number of columns in time series file */
00527   int iy;                  /* time series file column index */
00528   int ipt;                 /* time point index */
00529   float * ts_data;         /* input time series data */
00530 
00531 
00532   /*----- Read the time series file -----*/
00533   flim = mri_read_1D(ts_filename) ;
00534   if (flim == NULL)
00535     {                      /* filename -> ts_filename   11 Sep 2003 [rickr] */
00536       sprintf (message,  "Unable to read time series file: %s",  ts_filename);
00537       FIM_error (message);
00538     }
00539 
00540   
00541   /*----- Set pointer to data, and set dimensions -----*/
00542   far = MRI_FLOAT_PTR(flim);
00543   nx = flim->nx;
00544   ny = flim->ny; iy = 0 ;
00545   if( ny > 1 ){
00546     fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny);
00547   }
00548   
00549 
00550   /*----- Save the time series data -----*/
00551   *ts_length = nx;
00552   ts_data = (float *) malloc (sizeof(float) * nx);
00553   MTEST (ts_data);
00554   for (ipt = 0;  ipt < nx;  ipt++)
00555     ts_data[ipt] = far[ipt + iy*nx];   
00556   
00557   
00558   mri_free (flim);  flim = NULL;
00559 
00560   return (ts_data);
00561 }
00562 
00563 
00564 /*---------------------------------------------------------------------------*/
00565 /*
00566   Read multiple time series from specified file name.  This file name may have
00567   a column selector attached.
00568 */
00569 
00570 MRI_IMAGE * read_time_series 
00571 (
00572   char * ts_filename,      /* time series file name (plus column selectors) */
00573   int ** column_list       /* list of selected columns */
00574 )
00575 
00576 {
00577   char message[THD_MAX_NAME];    /* error message */
00578   char * cpt;                    /* pointer to column suffix */
00579   char filename[THD_MAX_NAME];   /* time series file name w/o column index */
00580   char subv[THD_MAX_NAME];       /* string containing column index */
00581   MRI_IMAGE * im, * flim;  /* pointers to image structures 
00582                               -- used to read 1D ASCII */
00583   float * far;             /* pointer to MRI_IMAGE floating point data */
00584   int nx;                  /* number of time points in time series */
00585   int ny;                  /* number of columns in time series file */
00586 
00587 
00588   /*----- Read the time series file -----*/
00589   flim = mri_read_1D(ts_filename) ;
00590   if (flim == NULL)
00591     {
00592       sprintf (message,  "Unable to read time series file: %s",  ts_filename);
00593       FIM_error (message);
00594     }
00595 
00596   
00597   /*----- Set pointer to data, and set dimensions -----*/
00598   far = MRI_FLOAT_PTR(flim);
00599   nx = flim->nx;
00600   ny = flim->ny;
00601   *column_list = NULL;   /* mri_read_1D does column selection */
00602 
00603   return (flim);
00604 }
00605 
00606 
00607 /*---------------------------------------------------------------------------*/
00608 /*
00609   Read the input data files.
00610 */
00611 
00612 void read_input_data
00613 (
00614   FIM_options * option_data,        /* fim program options */
00615   THD_3dim_dataset ** dset_time,    /* input 3d+time data set */
00616   THD_3dim_dataset ** mask_dset,    /* input mask data set */
00617   float ** fmri_data,               /* input fMRI time series data */
00618   int * fmri_length,                /* length of fMRI time series */
00619   MRI_IMAGE ** ort_array,           /* ort time series arrays */
00620   int ** ort_list,                  /* list of ort time series */
00621   MRI_IMAGE ** ideal_array,         /* ideal time series arrays */
00622   int ** ideal_list                 /* list of ideal time series */
00623 )
00624 
00625 {
00626   char message[THD_MAX_NAME];    /* error message */
00627 
00628   int num_ort_files;       /* number of ort time series files */
00629   int num_ideal_files;     /* number of ideal time series files */
00630   int is;                  /* time series file index */
00631   int polort;              /* degree of polynomial for baseline model */
00632   int num_ortts;           /* number of ort time series */
00633   int num_idealts;         /* number of ideal time series */
00634   int q;                   /* number of parameters in the baseline model */
00635   int p;                   /* number of parameters in the baseline model 
00636                               plus number of ideals */
00637 
00638 
00639   /*----- Initialize local variables -----*/
00640   polort          = option_data->polort;
00641   num_ort_files   = option_data->num_ort_files;
00642   num_ideal_files = option_data->num_ideal_files;
00643 
00644 
00645   /*----- Read the input fMRI measurement data -----*/
00646   if (option_data->input1D_filename != NULL)
00647     {
00648       /*----- Read the input fMRI 1D time series -----*/
00649       *fmri_data = read_one_time_series (option_data->input1D_filename, 
00650                                          fmri_length);
00651       if (*fmri_data == NULL)  
00652         { 
00653           sprintf (message,  "Unable to read time series file: %s", 
00654                    option_data->input1D_filename);
00655           FIM_error (message);
00656         }  
00657       *dset_time = NULL;
00658     }
00659 
00660   else if (option_data->input_filename != NULL)
00661     {
00662       /*----- Read the input 3d+time dataset -----*/
00663       *dset_time = THD_open_one_dataset (option_data->input_filename);
00664       if (!ISVALID_3DIM_DATASET(*dset_time))  
00665         { 
00666           sprintf (message,  "Unable to open data file: %s", 
00667                    option_data->input_filename);
00668           FIM_error (message);
00669         }  
00670       THD_load_datablock ((*dset_time)->dblk);
00671 
00672       if (option_data->mask_filename != NULL)
00673         {
00674           /*----- Read the input mask dataset -----*/
00675           *mask_dset = THD_open_dataset (option_data->mask_filename);
00676           if (!ISVALID_3DIM_DATASET(*mask_dset))  
00677             { 
00678               sprintf (message,  "Unable to open mask file: %s", 
00679                        option_data->mask_filename);
00680               FIM_error (message);
00681             }  
00682           THD_load_datablock ((*mask_dset)->dblk);
00683         }
00684     }
00685   else
00686     FIM_error ("Must specify input measurement data");
00687 
00688 
00689   /*----- Read the input ort time series files -----*/
00690   for (is = 0;  is < num_ort_files;  is++)
00691     {
00692       ort_array[is] = read_time_series (option_data->ort_filename[is], 
00693                                         &(ort_list[is]));
00694 
00695       if (ort_array[is] == NULL)
00696         {
00697           sprintf (message,  "Unable to read ort time series file: %s", 
00698                    option_data->ort_filename[is]);
00699           FIM_error (message);
00700         }
00701     }
00702 
00703   
00704   /*----- Read the input ideal time series files -----*/
00705   for (is = 0;  is < num_ideal_files;  is++)
00706     {
00707       ideal_array[is] = read_time_series (option_data->ideal_filename[is], 
00708                                           &(ideal_list[is]));
00709 
00710       if (ideal_array[is] == NULL)
00711         {
00712           sprintf (message,  "Unable to read ideal time series file: %s", 
00713                    option_data->ideal_filename[is]);
00714           FIM_error (message);
00715         }
00716     }
00717 
00718   
00719   /*----- Count number of ort and number of ideal time series -----*/
00720   num_ortts = 0;
00721   for (is = 0;  is < num_ort_files;  is++)
00722     {
00723       if (ort_list[is] == NULL)
00724         num_ortts += ort_array[is]->ny;
00725       else
00726         num_ortts += ort_list[is][0];
00727     }
00728   q = polort + 1 + num_ortts;
00729 
00730   num_idealts = 0;
00731   for (is = 0;  is < num_ideal_files;  is++)
00732     {
00733       if (ideal_list[is] == NULL)
00734         num_idealts += ideal_array[is]->ny;
00735       else
00736         num_idealts += ideal_list[is][0];
00737     }
00738   p = q + num_idealts;
00739 
00740   option_data->num_ortts = num_ortts;
00741   option_data->num_idealts = num_idealts;
00742   option_data->q = q;
00743   option_data->p = p;
00744 
00745 }
00746 
00747 
00748 /*---------------------------------------------------------------------------*/
00749 /*
00750   Routine to check whether one output file already exists.
00751 */
00752 
00753 void check_one_output_file 
00754 (
00755   THD_3dim_dataset * dset_time,     /* input 3d+time data set */
00756   char * filename                   /* name of output file */
00757 )
00758 
00759 {
00760   char message[THD_MAX_NAME];      /* error message */
00761   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00762   int ierror;                         /* number of errors in editing data */
00763 
00764   
00765   /*----- make an empty copy of input dataset -----*/
00766   new_dset = EDIT_empty_copy( dset_time ) ;
00767   
00768   
00769   ierror = EDIT_dset_items( new_dset ,
00770                             ADN_prefix , filename ,
00771                             ADN_label1 , filename ,
00772                             ADN_self_name , filename ,
00773                             ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 
00774                                                            GEN_FUNC_TYPE ,
00775                             ADN_none ) ;
00776   
00777   if( ierror > 0 )
00778     {
00779       sprintf (message,
00780                "*** %d errors in attempting to create output dataset!\n", 
00781                ierror);
00782       FIM_error (message);
00783     }
00784   
00785   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00786     {
00787       sprintf (message,
00788                "Output dataset file %s already exists "
00789                " -- cannot continue!\a\n",
00790                new_dset->dblk->diskptr->header_name);
00791       FIM_error (message);
00792     }
00793   
00794   /*----- deallocate memory -----*/   
00795   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00796   
00797 }
00798 
00799 
00800 /*---------------------------------------------------------------------------*/
00801 /*
00802   Routine to check whether output files already exist.
00803 */
00804 
00805 void check_output_files 
00806 (
00807   FIM_options * option_data,       /* fim program options */
00808   THD_3dim_dataset * dset_time     /* input 3d+time data set */
00809 )
00810 
00811 {
00812 
00813   if ((option_data->bucket_filename != NULL)
00814       && (option_data->input1D_filename == NULL))
00815     check_one_output_file (dset_time, option_data->bucket_filename);
00816   
00817 }
00818 
00819 
00820 /*---------------------------------------------------------------------------*/
00821 /*
00822   Routine to check for valid inputs.
00823 */
00824   
00825 void check_for_valid_inputs 
00826 (
00827   FIM_options * option_data,      /* fim program options */
00828   THD_3dim_dataset * dset_time,   /* input 3d+time data set */
00829   THD_3dim_dataset * mask_dset,   /* input mask data set */
00830   int fmri_length,                /* length of input fMRI time series */
00831   MRI_IMAGE ** ort_array,         /* ort time series arrays */
00832   MRI_IMAGE ** ideal_array        /* ideal time series arrays */
00833 )
00834 
00835 {
00836   char message[THD_MAX_NAME];  /* error message */
00837   int is;                  /* ideal index */
00838   int num_ort_files;       /* number of ort time series files */
00839   int num_ideal_files;     /* number of ideal time series files */
00840   int num_idealts;         /* number of ideal time series */
00841   int nt;                  /* number of images in input 3d+time dataset */
00842   int NFirst;              /* first image from input 3d+time dataset to use */
00843   int NLast;               /* last image from input 3d+time dataset to use */
00844   int N;                   /* number of usable time points */
00845 
00846 
00847   /*----- Initialize local variables -----*/
00848   if (option_data->input1D_filename != NULL)
00849     nt = fmri_length;
00850   else
00851     nt = DSET_NUM_TIMES (dset_time);
00852 
00853   num_ort_files   = option_data->num_ort_files;
00854   num_ideal_files = option_data->num_ideal_files;
00855   num_idealts     = option_data->num_idealts;
00856 
00857 
00858   NFirst = option_data->NFirst;
00859 
00860   NLast = option_data->NLast;   
00861   if (NLast > nt-1)  NLast = nt-1;
00862   option_data->NLast = NLast;
00863 
00864   N = NLast - NFirst + 1;
00865   option_data->N = N;
00866 
00867 
00868   /*----- Check number of ideal time series -----*/
00869   if (num_idealts < 1)  FIM_error ("No ideal time series?");
00870   if (num_idealts < 2)  option_data->output_type[FIM_BestIndex] = 0;
00871 
00872 
00873   /*----- If mask is used, check for compatible dimensions -----*/
00874   if (mask_dset != NULL)
00875     {
00876       if ( (DSET_NX(dset_time) != DSET_NX(mask_dset))
00877            || (DSET_NY(dset_time) != DSET_NY(mask_dset))
00878            || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) )
00879         {
00880           sprintf (message, "%s and %s have incompatible dimensions",
00881                    option_data->input_filename, option_data->mask_filename);
00882           FIM_error (message);
00883         }
00884 
00885       if (DSET_NVALS(mask_dset) != 1 )
00886         FIM_error ("Must specify 1 sub-brick from mask dataset");
00887     }
00888 
00889 
00890   /*----- Check lengths of ort time series -----*/
00891   for (is = 0;  is < num_ort_files;  is++)
00892     {
00893       if (ort_array[is]->nx < NLast+1)
00894         {
00895           sprintf (message, "Input ort time series file %s is too short",
00896                    option_data->ort_filename[is]);
00897           FIM_error (message);
00898         }
00899     }
00900 
00901  
00902   /*----- Check lengths of ideal time series -----*/
00903   for (is = 0;  is < num_ideal_files;  is++)
00904     {
00905       if (ideal_array[is]->nx < NLast+1)
00906         {
00907           sprintf (message, "Input ideal time series file %s is too short",
00908                    option_data->ideal_filename[is]);
00909           FIM_error (message);
00910         }
00911     }
00912 
00913  
00914   /*----- Check whether any of the output files already exist -----*/
00915   check_output_files (option_data, dset_time);
00916 
00917 }
00918 
00919 
00920 /*---------------------------------------------------------------------------*/
00921 /*
00922   Allocate memory for output volumes.
00923 */
00924 
00925 void allocate_memory 
00926 (
00927   FIM_options * option_data,  /* fim program options */
00928   THD_3dim_dataset * dset,    /* input 3d+time data set */
00929   float *** fim_params_vol    /* array of volumes containing fim parameters */
00930 )
00931 
00932 {
00933   int ip;                    /* parameter index */
00934   int nxyz;                  /* total number of voxels */
00935   int ixyz;                  /* voxel index */
00936 
00937 
00938   /*----- Initialize local variables -----*/
00939   nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
00940 
00941 
00942   /*----- Allocate memory space for fim parameters -----*/
00943   *fim_params_vol  = (float **) malloc (sizeof(float *) * MAX_OUTPUT_TYPE);   
00944   MTEST(*fim_params_vol);
00945 
00946   for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)
00947     {
00948       if (option_data->output_type[ip])
00949         {
00950           (*fim_params_vol)[ip]  = (float *) malloc (sizeof(float) * nxyz);
00951           MTEST((*fim_params_vol)[ip]);    
00952           for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00953             {
00954               (*fim_params_vol)[ip][ixyz]  = 0.0;
00955             }
00956         }
00957       else
00958         (*fim_params_vol)[ip] = NULL;
00959     }
00960 
00961 }
00962 
00963 
00964 
00965 
00966 /*---------------------------------------------------------------------------*/
00967 /*
00968   Perform all program initialization.
00969 */
00970 
00971 void initialize_program 
00972 (
00973   int argc,                         /* number of input arguments */
00974   char ** argv,                     /* array of input arguments */ 
00975   FIM_options ** option_data,       /* fim algorithm options */
00976   THD_3dim_dataset ** dset_time,    /* input 3d+time data set */
00977   THD_3dim_dataset ** mask_dset,    /* input mask data set */
00978   float ** fmri_data,               /* input fMRI time series data */
00979   int * fmri_length,                /* length of fMRI time series */
00980   MRI_IMAGE ** ort_array,           /* ort time series arrays */
00981   int ** ort_list,                  /* list of ort time series */
00982   MRI_IMAGE ** ideal_array,         /* ideal time series arrays */
00983   int ** ideal_list,                /* list of ideal time series */
00984   float *** fim_params_vol    /* array of volumes containing fim parameters */
00985 )
00986      
00987 {
00988   /*----- Allocate memory -----*/
00989   *option_data = (FIM_options *) malloc (sizeof(FIM_options));
00990 
00991    
00992   /*----- Get command line inputs -----*/
00993   get_options (argc, argv, *option_data);
00994 
00995 
00996   /*----- Read input data -----*/
00997   read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length,
00998                    ort_array, ort_list, ideal_array, ideal_list);
00999 
01000 
01001   /*----- Check for valid inputs -----*/
01002   check_for_valid_inputs (*option_data, *dset_time, *mask_dset, 
01003                           *fmri_length, ort_array, ideal_array);
01004   
01005 
01006   /*----- Allocate memory for output volumes -----*/
01007   if ((*option_data)->input1D_filename == NULL)
01008     allocate_memory (*option_data, *dset_time, fim_params_vol);
01009 
01010 }
01011 
01012 
01013 /*---------------------------------------------------------------------------*/
01014 /*
01015   Get the time series for one voxel from the AFNI 3d+time data set.
01016 */
01017 
01018 void extract_ts_array 
01019 (
01020   THD_3dim_dataset * dset_time,      /* input 3d+time dataset */
01021   int iv,                            /* get time series for this voxel */
01022   float * ts_array                   /* time series data for voxel #iv */
01023 )
01024 
01025 {
01026   MRI_IMAGE * im;          /* intermediate float data */
01027   float * ar;              /* pointer to float data */
01028   int ts_length;           /* length of input 3d+time data set */
01029   int it;                  /* time index */
01030 
01031 
01032   /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/
01033   im = THD_extract_series (iv, dset_time, 0);
01034 
01035 
01036   /*----- Verify extraction -----*/
01037   if (im == NULL)  FIM_error ("Unable to extract data from 3d+time dataset");
01038 
01039 
01040   /*----- Now extract time series from MRI_IMAGE -----*/
01041   ts_length = DSET_NUM_TIMES (dset_time);
01042   ar = MRI_FLOAT_PTR (im);
01043   for (it = 0;  it < ts_length;  it++)
01044     {
01045       ts_array[it] = ar[it];
01046     }
01047 
01048 
01049   /*----- Release memory -----*/
01050   mri_free (im);   im = NULL;
01051 
01052 }
01053 
01054 
01055 /*---------------------------------------------------------------------------*/
01056 /*
01057   Save results for this voxel.
01058 */
01059 
01060 void save_voxel 
01061 (
01062   int iv,                      /* current voxel index */      
01063   float * fim_params,          /* array of fim parameters */
01064   float ** fim_params_vol      /* array of volumes of fim output parameters */
01065 )
01066 
01067 {
01068   int ip;                   /* parameter index */ 
01069 
01070 
01071   /*----- Saved user requested fim parameters -----*/
01072   for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)
01073     {
01074       if (fim_params_vol[ip] != NULL)
01075         fim_params_vol[ip][iv]  = fim_params[ip];
01076     }
01077 
01078 }
01079 
01080 
01081 /*---------------------------------------------------------------------------*/
01082 /*
01083   Set fim threshold level.
01084 */
01085 
01086 float set_fim_thr_level 
01087 (
01088   int NFirst,                /* first usable data point */
01089   float fim_thr,             /* fim threshold (as proportion of mean) */
01090   THD_3dim_dataset * dset    /* input 3d+time data set */
01091 )
01092 
01093 {
01094   int nt;                    /* number of time points */
01095   int nxyz;                  /* total number of voxels */
01096   int ixyz;                  /* voxel index */
01097   double sum;                /* sum of voxel intensities */
01098   float fthr;                /* fim threshold (as intensity level) */
01099   float * ts_array = NULL;   /* time series data for individual voxel */
01100 
01101 
01102   /*----- Initialize local variables -----*/
01103   nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01104   nt = DSET_NUM_TIMES (dset);
01105 
01106   ts_array = (float *) malloc (sizeof(float) * nt);
01107   MTEST (ts_array);
01108 
01109   sum = 0.0;  /* Ides March 2004 [rickr] */
01110   /*----- Sum values of all voxels at initial time point -----*/
01111   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01112     {
01113       extract_ts_array (dset, ixyz, ts_array);
01114       sum += fabs(ts_array[NFirst]);
01115     }
01116 
01117 
01118   /*----- Set fim intensity level threshold -----*/
01119   fthr = fim_thr * sum / nxyz;
01120 
01121 
01122   free (ts_array);   ts_array = NULL;
01123 
01124   return (fthr);
01125 }
01126 
01127 
01128 /*---------------------------------------------------------------------------*/
01129 /*
01130   Calculate the best cross correlation for each voxel.
01131 */
01132 
01133 void calculate_results 
01134 (
01135   FIM_options * option_data,  /* fim program options */
01136   THD_3dim_dataset * dset,    /* input 3d+time data set */
01137   THD_3dim_dataset * mask,    /* input mask data set */
01138   float * fmri_data,          /* input fMRI time series data */
01139   int fmri_length,            /* length of fMRI time series */
01140   MRI_IMAGE ** ort_array,     /* ort time series arrays */
01141   int ** ort_list,            /* list of ort time series */
01142   MRI_IMAGE ** ideal_array,   /* ideal time series arrays */
01143   int ** ideal_list,          /* list of ideal time series */
01144   float ** fim_params_vol     /* array of volumes of fim output parameters */
01145 )
01146   
01147 {
01148   float * ts_array = NULL;    /* array of measured data for one voxel */
01149   float mask_val[1];          /* value of mask at current voxel */
01150   float fthr;                 /* internal mask threshold level */
01151 
01152   int q;                      /* number of parameters in the baseline model */
01153   int p;                      /* number of parameters in the baseline model 
01154                                  plus number of ideals */
01155   int m;                      /* parameter index */
01156   int n;                      /* data point index */
01157 
01158 
01159   matrix xdata;               /* independent variable matrix */
01160   matrix x_base;              /* extracted X matrix    for baseline model */
01161   matrix xtxinvxt_base;       /* matrix:  (1/(X'X))X'  for baseline model */
01162   matrix x_ideal[MAX_FILES];  /* extracted X matrices  for ideal models */
01163   matrix xtxinvxt_ideal[MAX_FILES];     
01164                               /* matrix:  (1/(X'X))X'  for ideal models */
01165   vector y;                   /* vector of measured data */       
01166 
01167   int ixyz;                   /* voxel index */
01168   int nxyz;                   /* number of voxels per dataset */
01169 
01170   int nt;                  /* number of images in input 3d+time dataset */
01171   int NFirst;              /* first image from input 3d+time dataset to use */
01172   int NLast;               /* last image from input 3d+time dataset to use */
01173   int N;                   /* number of usable data points */
01174 
01175   int num_ort_files;       /* number of ort time series files */
01176   int num_ideal_files;     /* number of ideal time series files */
01177   int polort;              /* degree of polynomial ort */
01178   int num_ortts;           /* number of ort time series */
01179   int num_idealts;         /* number of ideal time series */
01180   
01181   int i;                   /* data point index */
01182   int is;                  /* ideal index */
01183   int ilag;                /* time lag index */
01184   float stddev;            /* normalized parameter standard deviation */
01185   char * label;            /* string containing stat. summary of results */
01186 
01187   float * x_bot = NULL;    /* minimum of stimulus time series */
01188   float * x_ave = NULL;    /* average of stimulus time series */
01189   float * x_top = NULL;    /* maximum of stimulus time series */
01190   int * good_list = NULL;  /* list of good (usable) time points */ 
01191   float ** rarray = NULL;  /* ranked arrays of ideal time series */
01192   float FimParams[MAX_OUTPUT_TYPE];  /* output fim parameters */
01193 
01194 
01195   /*----- Initialize matrices and vectors -----*/
01196   matrix_initialize (&xdata);
01197   matrix_initialize (&x_base);
01198   matrix_initialize (&xtxinvxt_base);
01199   for (is =0;  is < MAX_FILES;  is++)
01200     {
01201       matrix_initialize (&x_ideal[is]);
01202       matrix_initialize (&xtxinvxt_ideal[is]);
01203     }
01204   vector_initialize (&y);
01205 
01206 
01207   /*----- Initialize local variables -----*/
01208   if (option_data->input1D_filename != NULL)
01209     {
01210       nxyz = 1;
01211       nt = fmri_length;
01212     }
01213   else
01214     {
01215       nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;       
01216       nt = DSET_NUM_TIMES (dset);
01217     }
01218 
01219   NFirst = option_data->NFirst;
01220   NLast = option_data->NLast;   
01221   N = option_data->N;
01222   p = option_data->p;
01223   q = option_data->q;
01224 
01225   polort          = option_data->polort;
01226   num_idealts     = option_data->num_idealts;
01227   num_ort_files   = option_data->num_ort_files;
01228   num_ideal_files = option_data->num_ideal_files;
01229 
01230 
01231   /*----- Allocate memory -----*/
01232   ts_array = (float *) malloc (sizeof(float) * nt);      MTEST (ts_array);
01233   x_bot = (float *)    malloc (sizeof(float) * p);       MTEST (x_bot);
01234   x_ave = (float *)    malloc (sizeof(float) * p);       MTEST (x_ave);
01235   x_top = (float *)    malloc (sizeof(float) * p);       MTEST (x_top);
01236   good_list = (int *)  malloc (sizeof(int) * N);         MTEST (good_list);
01237   rarray = (float **)  malloc (sizeof(float *) * num_idealts);  MTEST (rarray);
01238 
01239 
01240   /*----- Initialize the independent variable matrix -----*/
01241   N = init_indep_var_matrix (q, p, NFirst, N, polort, 
01242                              num_ort_files, num_ideal_files, 
01243                              ort_array, ort_list, ideal_array, ideal_list, 
01244                              x_bot, x_ave, x_top, good_list, &xdata);
01245   option_data->N = N;
01246 
01247 
01248   /*----- Initialize fim threshold level -----*/
01249   if (option_data->input1D_filename == NULL)
01250     fthr = set_fim_thr_level (good_list[0]+NFirst, option_data->fim_thr, dset);
01251 
01252 
01253   /*----- Initialization for the regression analysis -----*/
01254   init_regression_analysis (q, p, N, num_idealts, xdata, &x_base, 
01255                             &xtxinvxt_base, x_ideal, xtxinvxt_ideal, rarray);
01256 
01257 
01258   vector_create (N, &y);
01259 
01260   
01261   /*----- Loop over all voxels -----*/
01262   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01263     {
01264       /*----- Apply mask? -----*/
01265       if (mask != NULL)
01266         {
01267           extract_ts_array (mask, ixyz, mask_val);
01268           if (mask_val[0] == 0.0)  continue; 
01269         }
01270       
01271       
01272       /*----- Extract Y-data for this voxel -----*/
01273       if (option_data->input1D_filename != NULL)
01274         {
01275           for (i = 0;  i < N;  i++)
01276             y.elts[i] = fmri_data[good_list[i]+NFirst];
01277         }
01278       else
01279         {
01280           extract_ts_array (dset, ixyz, ts_array);
01281           if (fabs(ts_array[good_list[0]+NFirst]) < fthr)  
01282             continue;
01283           for (i = 0;  i < N;  i++)
01284             y.elts[i] = ts_array[good_list[i]+NFirst];
01285         }
01286       
01287 
01288       /*----- Perform the regression analysis for this voxel-----*/
01289       regression_analysis (N, q, num_idealts,
01290                            x_base, xtxinvxt_base, x_ideal, xtxinvxt_ideal, 
01291                            y, x_bot, x_ave, x_top, rarray, 
01292                            option_data->output_type, FimParams);
01293 
01294 
01295       /*----- Save results for this voxel -----*/
01296       if (option_data->input1D_filename == NULL)
01297         save_voxel (ixyz, FimParams, fim_params_vol);
01298       
01299       
01300       /*----- Report results for this voxel -----*/
01301       if ( ((fabs(FimParams[FIM_Correlation]) > option_data->cdisp) 
01302             && (option_data->cdisp >= 0.0))
01303            || (option_data->input1D_filename != NULL) )
01304         {
01305           printf ("\n\nResults for Voxel #%d: \n", ixyz);
01306           report_results (option_data->output_type, FimParams, &label);
01307           printf ("%s \n", label);
01308         }
01309       
01310     }  /*----- Loop over voxels -----*/
01311   
01312 
01313   /*----- Dispose of matrices and vectors -----*/
01314   vector_destroy (&y);
01315   for (is = 0;  is < MAX_FILES;  is++)
01316     {
01317       matrix_destroy (&xtxinvxt_ideal[is]);
01318       matrix_destroy (&x_ideal[is]);
01319     } 
01320   matrix_destroy (&xtxinvxt_base);
01321   matrix_destroy (&x_base);
01322   matrix_destroy (&xdata);
01323 
01324 
01325   /*----- Deallocate memory -----*/
01326   free (ts_array);     ts_array = NULL;
01327   free (x_bot);        x_bot = NULL;
01328   free (x_ave);        x_ave = NULL;
01329   free (x_top);        x_top = NULL;
01330   free (good_list);    good_list = NULL;
01331   for (is = 0;  is < num_idealts;  is++)
01332     {
01333       free (rarray[is]);   rarray[is] = NULL;
01334     }
01335   free (rarray);       rarray = NULL;
01336   
01337 }
01338 
01339 
01340 /*---------------------------------------------------------------------------*/
01341 /*
01342   Convert one volume to another type, autoscaling:
01343      nxy   = # voxels
01344      itype = input datum type
01345      ivol  = pointer to input volume
01346      otype = output datum type
01347      ovol  = pointer to output volume (again, must be pre-malloc-ed)
01348   Return value is the scaling factor used (0.0 --> no scaling).
01349 */
01350 
01351 float EDIT_coerce_autoscale_new( int nxyz ,
01352                                  int itype,void *ivol , int otype,void *ovol )
01353 {
01354   float fac=0.0 , top ;
01355   
01356   if( MRI_IS_INT_TYPE(otype) ){
01357     top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01358     if (top == 0.0)  fac = 0.0;
01359     else  fac = MRI_TYPE_maxval[otype]/top;
01360   }
01361   
01362   EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01363   return ( fac );
01364 }
01365 
01366 
01367 /*---------------------------------------------------------------------------*/
01368 /*
01369   Attach one sub-brick to output bucket data set.
01370 */
01371 
01372 void attach_sub_brick
01373 (
01374   THD_3dim_dataset * new_dset,      /* output bucket dataset */
01375   int ibrick,               /* sub-brick indices */
01376   float * volume,           /* volume of floating point data */
01377   int nxyz,                 /* total number of voxels */
01378   int brick_type,           /* indicates statistical type of sub-brick */
01379   char * brick_label,       /* character string label for sub-brick */
01380   int nsam, 
01381   int nfit, 
01382   int nort,                 /* degrees of freedom */
01383   short ** bar              /* bar[ib] points to data for sub-brick #ib */  
01384 )
01385 
01386 {
01387   const float EPSILON = 1.0e-10;
01388   float factor;             /* factor is new scale factor for sub-brick #ib */
01389 
01390 
01391   /*----- allocate memory for output sub-brick -----*/
01392   bar[ibrick]  = (short *) malloc (sizeof(short) * nxyz);
01393   MTEST (bar[ibrick]);
01394   factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01395                                       MRI_short, bar[ibrick]);
01396   
01397   if (factor < EPSILON)  factor = 0.0;
01398   else factor = 1.0 / factor;
01399   
01400 
01401   /*----- edit the sub-brick -----*/
01402   EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
01403   EDIT_BRICK_FACTOR (new_dset, ibrick, factor);
01404 
01405   if (brick_type == FUNC_COR_TYPE)
01406     EDIT_BRICK_TO_FICO (new_dset, ibrick, nsam, nfit, nort);
01407   
01408   
01409   /*----- attach bar[ib] to be sub-brick #ibrick -----*/
01410   EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]);
01411 
01412 }
01413 
01414 /*---------------------------------------------------------------------------*/
01415 /*
01416   Routine to write one bucket data set.
01417 */
01418 
01419 void write_bucket_data 
01420 (
01421   int argc,                         /* number of input arguments */
01422   char ** argv,                     /* array of input arguments */ 
01423   FIM_options * option_data,        /* fim program options */
01424   float ** fim_params_vol      /* array of volumes of fim output parameters */
01425 )
01426 
01427 {
01428   THD_3dim_dataset * old_dset = NULL;      /* prototype dataset */
01429   THD_3dim_dataset * new_dset = NULL;      /* output bucket dataset */
01430   char output_prefix[THD_MAX_NAME];     /* prefix name for bucket dataset */
01431   char output_session[THD_MAX_NAME];    /* directory for bucket dataset */
01432   int nbricks;              /* number of sub-bricks in bucket dataset */
01433   short ** bar = NULL;      /* bar[ib] points to data for sub-brick #ib */
01434 
01435   int brick_type;           /* indicates statistical type of sub-brick */
01436   int brick_coef;           /* regression coefficient index for sub-brick */
01437   char brick_label[THD_MAX_NAME]; /* character string label for sub-brick */
01438 
01439   int ierror;               /* number of errors in editing data */
01440   float * volume;           /* volume of floating point data */
01441 
01442   int N;                    /* number of usable data points */
01443   int q;                    /* number of parameters in the ideal model */
01444   int num_idealts;           /* number of ideal time series */
01445   int ip;                   /* parameter index */
01446   int nxyz;                 /* total number of voxels */
01447   int ibrick;               /* sub-brick index */
01448   int nsam; 
01449   int nfit; 
01450   int nort;                 /* degrees of freedom */
01451   char label[THD_MAX_NAME];   /* general label for sub-bricks */
01452 
01453 
01454   /*----- read prototype dataset -----*/
01455   old_dset = THD_open_one_dataset (option_data->input_filename);
01456 
01457     
01458   /*----- Initialize local variables -----*/
01459   nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;  
01460   num_idealts = option_data->num_idealts;
01461   q = option_data->q;
01462   N = option_data->N;
01463 
01464 
01465   /*----- Calculate number of sub-bricks in the bucket -----*/
01466   nbricks = 0;
01467   for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)
01468     if (option_data->output_type[ip])  nbricks++;
01469 
01470 
01471   strcpy (output_prefix, option_data->bucket_filename);
01472   strcpy (output_session, "./");
01473   
01474  
01475   /*----- allocate memory -----*/
01476   bar  = (short **) malloc (sizeof(short *) * nbricks);
01477   MTEST (bar);
01478   
01479 
01480   /*-- make an empty copy of prototype dataset, for eventual output --*/
01481   new_dset = EDIT_empty_copy (old_dset);
01482 
01483 
01484   /*----- Record history of dataset -----*/
01485   tross_Copy_History( old_dset , new_dset ) ;
01486   tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
01487   sprintf (label, "Output prefix: %s", output_prefix);
01488   tross_Append_History ( new_dset, label);
01489 
01490   
01491   /*----- delete prototype dataset -----*/ 
01492   THD_delete_3dim_dataset( old_dset , False );  old_dset = NULL ;
01493 
01494 
01495   /*----- Modify some structural properties.  Note that the nbricks
01496           just make empty sub-bricks, without any data attached. -----*/
01497   ierror = EDIT_dset_items (new_dset,
01498                             ADN_prefix,          output_prefix,
01499                             ADN_directory_name,  output_session,
01500                             ADN_type,            HEAD_FUNC_TYPE,
01501                             ADN_func_type,       FUNC_BUCK_TYPE,
01502                             ADN_datum_all,       MRI_short ,   
01503                             ADN_ntt,             0,               /* no time */
01504                             ADN_nvals,           nbricks,
01505                             ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,  
01506                             ADN_none ) ;
01507   
01508   if( ierror > 0 )
01509     {
01510       fprintf(stderr, 
01511               "*** %d errors in attempting to create bucket dataset!\n", 
01512               ierror);
01513       exit(1);
01514     }
01515   
01516   if (THD_is_file(DSET_HEADNAME(new_dset))) 
01517     {
01518       fprintf(stderr,
01519               "*** Output dataset file %s already exists--cannot continue!\n",
01520               DSET_HEADNAME(new_dset));
01521       exit(1);
01522     }
01523   
01524 
01525   /*----- Attach individual sub-bricks to the bucket dataset -----*/
01526   ibrick = 0;
01527   for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)        
01528     {                                 
01529       if (option_data->output_type[ip] == 0)  continue;
01530 
01531       strcpy (brick_label, OUTPUT_TYPE_name[ip]);
01532 
01533       if (ip == FIM_Correlation)
01534         {
01535           brick_type = FUNC_COR_TYPE;
01536           nsam = N;  nort = q;
01537           if (num_idealts > 1)  nfit = 2;
01538           else                  nfit = 1;
01539         }
01540       else if (ip == FIM_SpearmanCC)
01541         {
01542 #if 0
01543           brick_type = FUNC_THR_TYPE;
01544 #else
01545           brick_type = FUNC_COR_TYPE;
01546 #endif
01547           nsam = N;  nort = q;
01548           if (num_idealts > 1)  nfit = 2;
01549           else                  nfit = 1;
01550         } 
01551       else if (ip == FIM_QuadrantCC)
01552         {
01553 #if 0
01554           brick_type = FUNC_THR_TYPE;
01555 #else
01556           brick_type = FUNC_COR_TYPE;
01557 #endif
01558           nsam = N;  nort = q;
01559           if (num_idealts > 1)  nfit = 2;
01560           else                  nfit = 1;
01561         }
01562       else
01563         {
01564           brick_type = FUNC_FIM_TYPE;
01565           nsam = 0;  nfit = 0;  nort = 0;
01566         }
01567 
01568       volume = fim_params_vol[ip];                
01569       attach_sub_brick (new_dset, ibrick, volume, nxyz, 
01570                         brick_type, brick_label, nsam, nfit, nort, bar);  
01571 
01572       ibrick++;
01573     }
01574 
01575 
01576   /*----- write bucket data set -----*/
01577   THD_load_statistics (new_dset);
01578   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01579   fprintf(stderr,"Wrote bucket dataset ");
01580   fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01581 
01582   
01583   /*----- deallocate memory -----*/   
01584   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01585 
01586 }
01587 
01588 
01589 /*---------------------------------------------------------------------------*/
01590 /*
01591   Write out the user requested output files.
01592 */
01593 
01594 void output_results
01595 (
01596   int argc,                         /* number of input arguments */
01597   char ** argv,                     /* array of input arguments */ 
01598   FIM_options * option_data,        /* fim algorithm options */
01599   float ** fim_params_vol      /* array of volumes of fim output parameters */
01600 )
01601 
01602 {
01603   int q;                    /* number of parameters in baseline model */
01604   int num_idealts;           /* number of ideal time series */
01605   int ib;                   /* sub-brick index */
01606   int is;                   /* ideal index */
01607   int ts_length;            /* length of impulse reponse function */
01608   int N;                    /* number of usable data points */
01609 
01610 
01611   /*----- Initialize local variables -----*/
01612   q = option_data->polort + 1;
01613   num_idealts = option_data->num_idealts;
01614   N = option_data->N;
01615 
01616 
01617   /*----- Write the bucket dataset -----*/
01618   if (option_data->bucket_filename != NULL)
01619     write_bucket_data (argc, argv, option_data, fim_params_vol);
01620 
01621 }
01622 
01623 
01624 /*---------------------------------------------------------------------------*/
01625 
01626 void terminate_program
01627 (
01628   FIM_options ** option_data,   /* fim program options */
01629   MRI_IMAGE ** ort_array,           /* ort time series arrays */
01630   int ** ort_list,                  /* list of ort time series */
01631   MRI_IMAGE ** ideal_array,         /* ideal time series arrays */
01632   int ** ideal_list,                /* list of ideal time series */
01633   float *** fim_params_vol      /* array of volumes of fim output parameters */
01634 )
01635 
01636 {
01637   int num_idealts;           /* number of ideal time series */
01638   int ip;                   /* parameter index */
01639   int is;                   /* ideal index */
01640 
01641 
01642   /*----- Initialize local variables -----*/
01643   num_idealts = (*option_data)->num_idealts;
01644 
01645 
01646   /*----- Deallocate memory for option data -----*/   
01647   free (*option_data);  *option_data = NULL;
01648 
01649 
01650   /*----- Deallocate memory for ideal time series -----*/
01651   /*
01652   for (is = 0;  is < num_idealts;  is++)
01653     { free (ideal[is]);  ideal[is] = NULL; } 
01654   */
01655 
01656 
01657   /*----- Deallocate space for volume data -----*/
01658   if (*fim_params_vol != NULL)
01659     {
01660       for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)
01661         {
01662           if ((*fim_params_vol)[ip] != NULL)
01663             { free ((*fim_params_vol)[ip]);   (*fim_params_vol)[ip] = NULL; }
01664         }
01665 
01666       free (*fim_params_vol);   *fim_params_vol  = NULL; 
01667     }
01668 
01669 }
01670 
01671 
01672 /*---------------------------------------------------------------------------*/
01673 
01674 int main 
01675 (
01676   int argc,                /* number of input arguments */
01677   char ** argv             /* array of input arguments */ 
01678 )
01679 
01680 {
01681   FIM_options * option_data;              /* fim algorithm options */
01682   THD_3dim_dataset * dset_time = NULL;    /* input 3d+time data set */
01683   THD_3dim_dataset * mask_dset = NULL;    /* input mask data set */
01684   float * fmri_data = NULL;               /* input fMRI time series data */
01685   int fmri_length;                        /* length of fMRI time series */
01686   MRI_IMAGE * ort_array[MAX_FILES];       /* ideal time series arrays */
01687   int * ort_list[MAX_FILES];              /* list of ideal time series */
01688   MRI_IMAGE * ideal_array[MAX_FILES];     /* ideal time series arrays */
01689   int * ideal_list[MAX_FILES];            /* list of ideal time series */
01690 
01691   float ** fim_params_vol = NULL;
01692                                 /* array of volumes of fim output parameters */
01693 
01694   
01695   /*----- Identify software -----*/
01696   printf ("\n\n");
01697   printf ("Program: %s \n", PROGRAM_NAME);
01698   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
01699   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01700   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01701   printf ("\n");
01702 
01703    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
01704 
01705    { int new_argc ; char ** new_argv ;
01706      addto_args( argc , argv , &new_argc , &new_argv ) ;
01707      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
01708    }
01709 
01710   
01711   /*----- Program initialization -----*/
01712   initialize_program (argc, argv, &option_data, &dset_time, &mask_dset, 
01713                       &fmri_data, &fmri_length, 
01714                       ort_array, ort_list, ideal_array, ideal_list, 
01715                       &fim_params_vol);
01716 
01717 
01718   /*----- Perform fim analysis -----*/
01719   calculate_results (option_data, dset_time, mask_dset, 
01720                      fmri_data, fmri_length,
01721                      ort_array, ort_list, ideal_array, ideal_list, 
01722                      fim_params_vol);
01723   
01724 
01725   /*----- Deallocate memory for input datasets -----*/   
01726   if (dset_time != NULL)  
01727     { THD_delete_3dim_dataset (dset_time, False);  dset_time = NULL; }
01728   if (mask_dset != NULL)  
01729     { THD_delete_3dim_dataset (mask_dset, False);  mask_dset = NULL; }
01730 
01731 
01732   /*----- Write requested output files -----*/
01733   if (option_data->input1D_filename == NULL)
01734     output_results (argc, argv, option_data, fim_params_vol);
01735 
01736 
01737   /*----- Terminate program -----*/
01738   terminate_program (&option_data, 
01739                      ort_array, ort_list, ideal_array, ideal_list, 
01740                      &fim_params_vol); 
01741 
01742   exit(0);
01743 }
01744 
01745 
01746 
01747 
01748 
01749 
01750 
01751 
01752