00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 2000-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 perform wavelet analysis of an FMRI 3d+time dataset.
00010 
00011   File:    3dWavelets.c
00012   Author:  B. Douglas Ward
00013   Date:    28 March 2000
00014 
00015   Mod:     Added call to AFNI_logger.
00016   Date:    15 August 2001
00017 
00018   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
00019   Date:    02 December 2002
00020 */
00021 
00022 /*---------------------------------------------------------------------------*/
00023 /*
00024   Software identification.
00025 */
00026 
00027 #define PROGRAM_NAME "3dWavelets"                    /* name of this program */
00028 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00029 #define PROGRAM_INITIAL "28 March 2000"   /* date of initial program release */
00030 #define PROGRAM_LATEST "02 December 2002"   /* date of last program revision */
00031 
00032 /*---------------------------------------------------------------------------*/
00033 /*
00034   Include header files and source code.
00035 */
00036 
00037 #include "mrilib.h"
00038 #include "Wavelets.h"
00039 #include "Wavelets.c"
00040 
00041 
00042 /*---------------------------------------------------------------------------*/
00043 /*
00044   Global constants. 
00045 */
00046 
00047 #define MAX_NAME_LENGTH THD_MAX_NAME    /* max. string length for file names */
00048 #define MAX_FILTERS 20               /* max. number of filter specifications */
00049 
00050 
00051 /*---------------------------------------------------------------------------*/
00052 /*
00053   Data structure for operator inputs.
00054 */
00055 
00056 typedef struct WA_options
00057 { 
00058   int NFirst;              /* first image from input 3d+time dataset to use */
00059   int NLast;               /* last image from input 3d+time dataset to use */
00060   int N;                   /* number of usable data points from input data */
00061   int f;                   /* number of parameters removed by the filter */
00062   int p;                   /* number of parameters in the full model */
00063   int q;                   /* number of parameters in the baseline model */
00064 
00065   float fdisp;             /* minimum F-statistic for display */ 
00066   char * input_filename;   /* input 3d+time dataset */
00067   char * mask_filename;    /* input mask dataset */
00068   char * input1D_filename; /* input fMRI measurement time series */
00069 
00070   int wavelet_type;        /* indicates type of wavelet */        
00071 
00072   int num_stop_filters;             /* number of wavelet stop filters */
00073   int stop_band[MAX_FILTERS];       /* wavelet filter stop band */ 
00074   int stop_mintr[MAX_FILTERS];      /* min. time for stop band */
00075   int stop_maxtr[MAX_FILTERS];      /* max. time for stop band */
00076   float * stop_filter;              /* select wavelet coefs. to stop */
00077 
00078   int num_base_filters;             /* number of wavelet base filters */
00079   int base_band[MAX_FILTERS];       /* wavelet filter base band */ 
00080   int base_mintr[MAX_FILTERS];      /* min. time for base band */
00081   int base_maxtr[MAX_FILTERS];      /* max. time for base band */
00082   float * base_filter;              /* select wavelet coefs. for baseline */
00083 
00084   int num_sgnl_filters;             /* number of wavelet signal filters */
00085   int sgnl_band[MAX_FILTERS];       /* wavelet filter signal band */ 
00086   int sgnl_mintr[MAX_FILTERS];      /* min. time for signal band */
00087   int sgnl_maxtr[MAX_FILTERS];      /* max. time for signal band */
00088   float * sgnl_filter;              /* select wavelet coefs. for signal */
00089 
00090   char * bucket_filename;   /* bucket dataset file name */
00091   char * coefts_filename;   /* forward wavelet transform coefficients output */
00092   char * fitts_filename;    /* full model time series 3d+time output */
00093   char * sgnlts_filename;   /* signal model time series 3d+time output */
00094   char * errts_filename;    /* residual error time series 3d+time output */
00095 
00096   int fout;                 /* flag to output F-statistics */
00097   int rout;                 /* flag to output R^2 statistics */
00098   int cout;                 /* flag to output wavelet coefficients */
00099   int vout;                 /* flag to output variance map */
00100   int stat_first;           /* flag to output full model stats first */
00101 
00102 } WA_options;
00103 
00104 
00105 /*---------------------------------------------------------------------------*/
00106 /*
00107   Routine to display 3dWavelets help menu.
00108 */
00109 
00110 void display_help_menu()
00111 {
00112   printf (
00113     "Program to perform wavelet analysis of an FMRI 3d+time dataset.        \n"
00114     "                                                                       \n"
00115     "Usage:                                                                 \n"
00116     "3dWavelets                                                             \n"
00117     "-type wname          wname = name of wavelet to use for the analysis   \n"
00118     "                     At present, there are only two choices for wname: \n"
00119     "                        Haar  -->  Haar wavelets                       \n"
00120     "                        Daub  -->  Daubechies wavelets                 \n"
00121     "-input fname         fname = filename of 3d+time input dataset         \n"
00122     "[-input1D dname]     dname = filename of single (fMRI) .1D time series \n"
00123     "[-mask mname]        mname = filename of 3d mask dataset               \n"
00124     "[-nfirst fnum]       fnum = number of first dataset image to use in    \n"
00125     "                       the wavelet analysis. (default = 0)             \n"
00126     "[-nlast  lnum]       lnum = number of last dataset image to use in     \n"
00127     "                       the wavelet analysis. (default = last)          \n"
00128     "[-fdisp fval]        Write (to screen) results for those voxels        \n"
00129     "                       whose F-statistic is >= fval                    \n"
00130     "                                                                       \n"
00131     "Filter options:                                                        \n"
00132     "[-filt_stop band mintr maxtr] Specify wavelet coefs. to set to zero    \n"
00133     "[-filt_base band mintr maxtr] Specify wavelet coefs. for baseline model\n"
00134     "[-filt_sgnl band mintr maxtr] Specify wavelet coefs. for signal model  \n"
00135     "     where  band  = frequency band                                     \n"
00136     "            mintr = min. value for time window (in TR)                 \n"
00137     "            maxtr = max. value for time window (in TR)                 \n"
00138     "                                                                       \n"
00139     "Output options:                                                        \n"
00140     "[-coefts cprefix]   cprefix = prefix of 3d+time output dataset which   \n"
00141     "                       will contain the forward wavelet transform      \n"
00142     "                       coefficients                                    \n"
00143     "                                                                       \n"
00144     "[-fitts  fprefix]   fprefix = prefix of 3d+time output dataset which   \n"
00145     "                       will contain the full model time series fit     \n"
00146     "                       to the input data                               \n"
00147     "                                                                       \n"
00148     "[-sgnlts sprefix]   sprefix = prefix of 3d+time output dataset which   \n"
00149     "                       will contain the signal model time series fit   \n"
00150     "                       to the input data                               \n"
00151     "                                                                       \n"
00152     "[-errts  eprefix]   eprefix = prefix of 3d+time output dataset which   \n"
00153     "                       will contain the residual error time series     \n"
00154     "                       from the full model fit to the input data       \n"
00155     "                                                                       \n"
00156     "The following options control the contents of the bucket dataset:      \n"
00157     "[-fout]            Flag to output the F-statistics                     \n"
00158     "[-rout]            Flag to output the R^2 statistics                   \n"
00159     "[-cout]            Flag to output the full model wavelet coefficients  \n"
00160     "[-vout]            Flag to output the sample variance (MSE) map        \n"
00161     "                                                                       \n"
00162     "[-stat_first]      Flag to specify that the full model statistics will \n"
00163     "                     appear prior to the wavelet coefficients in the   \n"
00164     "                     bucket dataset output                             \n"
00165     "                                                                       \n"
00166     "[-bucket bprefix]  bprefix = prefix of AFNI 'bucket' dataset containing\n"
00167     "                     parameters of interest, such as the F-statistic   \n"
00168     "                     for significance of the wavelet signal model.     \n"
00169     );
00170   
00171   exit(0);
00172 }
00173 
00174 
00175 /*---------------------------------------------------------------------------*/
00176 /*
00177   Routine to initialize the input options.
00178 */
00179  
00180 void initialize_options 
00181 (
00182   WA_options * option_data    /* wavelet analysis program options */
00183 )
00184  
00185 {
00186   int is;                     /* filter index */
00187 
00188 
00189   /*----- Initialize default values -----*/
00190   option_data->NFirst = 0;
00191   option_data->NLast  = 32767;
00192   option_data->N      = 0;
00193   option_data->f      = 0;
00194   option_data->p      = 0;
00195   option_data->q      = 0;
00196   option_data->fdisp = -1.0;
00197   option_data->wavelet_type = 0;
00198 
00199 
00200   /*----- Initialize filter specifications -----*/
00201   option_data->num_stop_filters = 0;
00202   option_data->num_base_filters = 0;
00203   option_data->num_sgnl_filters = 0;
00204   for (is = 0;  is < MAX_FILTERS;  is++)
00205     {
00206       option_data->stop_band[is] = 0;
00207       option_data->stop_mintr[is] = 0;
00208       option_data->stop_maxtr[is] = 0;
00209       option_data->base_band[is] = 0;
00210       option_data->base_mintr[is] = 0;
00211       option_data->base_maxtr[is] = 0;
00212       option_data->sgnl_band[is] = 0;
00213       option_data->sgnl_mintr[is] = 0;
00214       option_data->sgnl_maxtr[is] = 0;
00215     }
00216   option_data->stop_filter = NULL;
00217   option_data->base_filter = NULL;
00218   option_data->sgnl_filter = NULL;
00219 
00220 
00221   /*----- Initialize output flags -----*/
00222   option_data->fout = 0;
00223   option_data->rout = 0;
00224   option_data->cout = 0;
00225   option_data->vout = 0;
00226   option_data->stat_first = 0;
00227 
00228 
00229   /*----- Initialize file name character strings -----*/
00230   option_data->input_filename = NULL;
00231   option_data->mask_filename = NULL;  
00232   option_data->input1D_filename = NULL;
00233   option_data->bucket_filename = NULL;
00234   option_data->coefts_filename = NULL;
00235   option_data->fitts_filename = NULL;
00236   option_data->sgnlts_filename = NULL;
00237   option_data->errts_filename = NULL;
00238 
00239 }
00240 
00241 
00242 /*---------------------------------------------------------------------------*/
00243 /*
00244   Routine to get user specified input options.
00245 */
00246 
00247 void get_options
00248 (
00249   int argc,                        /* number of input arguments */
00250   char ** argv,                    /* array of input arguments */ 
00251   WA_options * option_data         /* wavelet program options */
00252 )
00253 
00254 {
00255   int nopt = 1;                     /* input option argument counter */
00256   int ival, index;                  /* integer input */
00257   float fval;                       /* float input */
00258   char message[MAX_NAME_LENGTH];    /* error message */
00259   int ifilter;                      /* filter index */
00260 
00261 
00262   /*-- addto the arglist, if user wants to --*/
00263    machdep() ; 
00264   { int new_argc ; char ** new_argv ;
00265     addto_args( argc , argv , &new_argc , &new_argv ) ;
00266     if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00267   }
00268   
00269 
00270   /*----- does user request help menu? -----*/
00271   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00272 
00273 
00274   /*----- Add to program log -----*/
00275   AFNI_logger (PROGRAM_NAME,argc,argv); 
00276 
00277   
00278   /*----- initialize the input options -----*/
00279   initialize_options (option_data); 
00280 
00281   
00282   /*----- main loop over input options -----*/
00283   while (nopt < argc )
00284     {
00285 
00286       /*-----   -type wname   -----*/
00287       if (strncmp(argv[nopt], "-type", 5) == 0)
00288         {
00289           nopt++;
00290           if (nopt >= argc)  WA_error ("need argument after -type ");
00291           for (ival = 0;  ival < MAX_WAVELET_TYPE;  ival++)
00292             if (strncmp(argv[nopt], WAVELET_TYPE_name[ival], 4) == 0)  break;
00293           if (ival >= MAX_WAVELET_TYPE)
00294             {
00295               sprintf(message,"Unrecognized wavelet type: %s\n", argv[nopt]);
00296               WA_error (message);
00297             }
00298           option_data->wavelet_type = ival;
00299           nopt++;
00300           continue;
00301         }
00302       
00303 
00304       /*-----   -input filename   -----*/
00305       if (strncmp(argv[nopt], "-input", 8) == 0)
00306         {
00307           nopt++;
00308           if (nopt >= argc)  WA_error ("need argument after -input ");
00309           option_data->input_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00310           MTEST (option_data->input_filename);
00311           strcpy (option_data->input_filename, argv[nopt]);
00312           nopt++;
00313           continue;
00314         }
00315       
00316 
00317       /*-----   -input1D filename   -----*/
00318       if (strncmp(argv[nopt], "-input1D", 8) == 0)
00319         {
00320           nopt++;
00321           if (nopt >= argc)  WA_error ("need argument after -input1D ");
00322           option_data->input1D_filename = 
00323             malloc (sizeof(char)*MAX_NAME_LENGTH);
00324           MTEST (option_data->input1D_filename);
00325           strcpy (option_data->input1D_filename, argv[nopt]);
00326           nopt++;
00327           continue;
00328         }
00329       
00330 
00331       /*-----   -mask filename   -----*/
00332       if (strncmp(argv[nopt], "-mask", 5) == 0)
00333         {
00334           nopt++;
00335           if (nopt >= argc)  WA_error ("need argument after -mask ");
00336           option_data->mask_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00337           MTEST (option_data->mask_filename);
00338           strcpy (option_data->mask_filename, argv[nopt]);
00339           nopt++;
00340           continue;
00341         }
00342       
00343 
00344       /*-----   -nfirst num  -----*/
00345       if (strncmp(argv[nopt], "-nfirst", 7) == 0)
00346         {
00347           nopt++;
00348           if (nopt >= argc)  WA_error ("need argument after -nfirst ");
00349           sscanf (argv[nopt], "%d", &ival);
00350           if (ival < 0)
00351             WA_error ("illegal argument after -nfirst ");
00352           option_data->NFirst = ival;
00353           nopt++;
00354           continue;
00355         }
00356 
00357 
00358       /*-----   -nlast num  -----*/
00359       if (strncmp(argv[nopt], "-nlast", 6) == 0)
00360         {
00361           nopt++;
00362           if (nopt >= argc)  WA_error ("need argument after -nlast ");
00363           sscanf (argv[nopt], "%d", &ival);
00364           if (ival < 0)
00365             WA_error ("illegal argument after -nlast ");
00366           option_data->NLast = ival;
00367           nopt++;
00368           continue;
00369         }
00370 
00371 
00372       /*-----   -fdisp fval   -----*/
00373       if (strncmp(argv[nopt], "-fdisp", 6) == 0)
00374         {
00375           nopt++;
00376           if (nopt >= argc)  WA_error ("need argument after -fdisp ");
00377           sscanf (argv[nopt], "%f", &fval); 
00378           option_data->fdisp = fval;
00379           nopt++;
00380           continue;
00381         }
00382       
00383 
00384       /*-----   -filt_stop band mintr maxtr   -----*/
00385       if (strncmp(argv[nopt], "-filt_stop", 10) == 0)
00386         {
00387           nopt++;
00388           if (nopt+2 >= argc)  
00389             WA_error ("need 3 arguments after -filt_stop");
00390 
00391           ifilter = option_data->num_stop_filters;
00392 
00393           sscanf (argv[nopt], "%d", &ival);
00394           if (ival < -1)
00395             WA_error ("-filt_stop band mintr maxtr   Require: -1 <= band");
00396           option_data->stop_band[ifilter] = ival;
00397           nopt++;
00398 
00399           sscanf (argv[nopt], "%d", &ival);
00400           if (ival < 0)
00401             WA_error ("-filt_stop band mintr maxtr   Require: 0 <= mintr");
00402           option_data->stop_mintr[ifilter] = ival;
00403           nopt++;
00404 
00405           sscanf (argv[nopt], "%d", &ival);
00406           if (ival < 0)
00407             WA_error ("-filt_stop band mintr maxtr   Require: 0 <= maxtr");
00408           option_data->stop_maxtr[ifilter] = ival;
00409           nopt++;
00410 
00411           option_data->num_stop_filters++;
00412           continue;
00413         }
00414       
00415 
00416       /*-----   -filt_base band mintr maxtr   -----*/
00417       if (strncmp(argv[nopt], "-filt_base", 10) == 0)
00418         {
00419           nopt++;
00420           if (nopt+2 >= argc)  
00421             WA_error ("need 3 arguments after -filt_base");
00422 
00423           ifilter = option_data->num_base_filters;
00424 
00425           sscanf (argv[nopt], "%d", &ival);
00426           if (ival < -1)
00427             WA_error ("-filt_base band mintr maxtr   Require: -1 <= band");
00428           option_data->base_band[ifilter] = ival;
00429           nopt++;
00430 
00431           sscanf (argv[nopt], "%d", &ival);
00432           if (ival < 0)
00433             WA_error ("-filt_base band mintr maxtr   Require: 0 <= mintr");
00434           option_data->base_mintr[ifilter] = ival;
00435           nopt++;
00436 
00437           sscanf (argv[nopt], "%d", &ival);
00438           if (ival < 0)
00439             WA_error ("-filt_base band mintr maxtr   Require: 0 <= maxtr");
00440           option_data->base_maxtr[ifilter] = ival;
00441           nopt++;
00442 
00443           option_data->num_base_filters++;
00444           continue;
00445         }
00446       
00447 
00448       /*-----   -filt_sgnl band mintr maxtr   -----*/
00449       if (strncmp(argv[nopt], "-filt_sgnl", 10) == 0)
00450         {
00451           nopt++;
00452           if (nopt+2 >= argc)  
00453             WA_error ("need 3 arguments after -filt_sgnl");
00454 
00455           ifilter = option_data->num_sgnl_filters;
00456 
00457           sscanf (argv[nopt], "%d", &ival);
00458           if (ival < -1)
00459             WA_error ("-filt_sgnl band mintr maxtr   Require: -1 <= band");
00460           option_data->sgnl_band[ifilter] = ival;
00461           nopt++;
00462 
00463           sscanf (argv[nopt], "%d", &ival);
00464           if (ival < 0)
00465             WA_error ("-filt_sgnl band mintr maxtr   Require: 0 <= mintr");
00466           option_data->sgnl_mintr[ifilter] = ival;
00467           nopt++;
00468 
00469           sscanf (argv[nopt], "%d", &ival);
00470           if (ival < 0)
00471             WA_error ("-filt_sgnl band mintr maxtr   Require: 0 <= maxtr");
00472           option_data->sgnl_maxtr[ifilter] = ival;
00473           nopt++;
00474 
00475           option_data->num_sgnl_filters++;
00476           continue;
00477         }
00478       
00479 
00480       /*-----   -fout   -----*/
00481       if (strncmp(argv[nopt], "-fout", 5) == 0)
00482         {
00483           option_data->fout = 1;
00484           nopt++;
00485           continue;
00486         }
00487       
00488 
00489       /*-----   -rout   -----*/
00490       if (strncmp(argv[nopt], "-rout", 5) == 0)
00491         {
00492           option_data->rout = 1;
00493           nopt++;
00494           continue;
00495         }
00496       
00497 
00498       /*-----   -cout   -----*/
00499       if (strncmp(argv[nopt], "-cout", 5) == 0)
00500         {
00501           option_data->cout = 1;
00502           nopt++;
00503           continue;
00504         }
00505       
00506 
00507       /*-----   -vout   -----*/
00508       if (strncmp(argv[nopt], "-vout", 5) == 0)
00509         {
00510           option_data->vout = 1;
00511           nopt++;
00512           continue;
00513         }
00514       
00515 
00516       /*-----   -stat_first   -----*/
00517       if (strncmp(argv[nopt], "-stat_first", 11) == 0)
00518         {
00519           option_data->stat_first = 1;
00520           nopt++;
00521           continue;
00522         }
00523       
00524 
00525       /*-----   -bucket filename   -----*/
00526       if (strncmp(argv[nopt], "-bucket", 6) == 0)
00527         {
00528           nopt++;
00529           if (nopt >= argc)  WA_error ("need file prefixname after -bucket ");
00530           option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00531           MTEST (option_data->bucket_filename);
00532           strcpy (option_data->bucket_filename, argv[nopt]);
00533           nopt++;
00534           continue;
00535         }
00536       
00537 
00538       /*-----   -coefts filename   -----*/
00539       if (strncmp(argv[nopt], "-coefts", 8) == 0)
00540         {
00541           nopt++;
00542           if (nopt >= argc)  WA_error ("need file prefixname after -coefts ");
00543           option_data->coefts_filename = 
00544             malloc (sizeof(char)*MAX_NAME_LENGTH);
00545           MTEST (option_data->coefts_filename);
00546           strcpy (option_data->coefts_filename, argv[nopt]);
00547           nopt++;
00548           continue;
00549         }
00550       
00551 
00552       /*-----   -fitts filename   -----*/
00553       if (strncmp(argv[nopt], "-fitts", 7) == 0)
00554         {
00555           nopt++;
00556           if (nopt >= argc)  WA_error ("need file prefixname after -fitts ");
00557           option_data->fitts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00558           MTEST (option_data->fitts_filename);
00559           strcpy (option_data->fitts_filename, argv[nopt]);
00560           nopt++;
00561           continue;
00562         }
00563       
00564 
00565       /*-----   -sgnlts filename   -----*/
00566       if (strncmp(argv[nopt], "-sgnlts", 7) == 0)
00567         {
00568           nopt++;
00569           if (nopt >= argc)  WA_error ("need file prefixname after -sgnlts ");
00570           option_data->sgnlts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00571           MTEST (option_data->sgnlts_filename);
00572           strcpy (option_data->sgnlts_filename, argv[nopt]);
00573           nopt++;
00574           continue;
00575         }
00576       
00577 
00578       /*-----   -errts filename   -----*/
00579       if (strncmp(argv[nopt], "-errts", 6) == 0)
00580         {
00581           nopt++;
00582           if (nopt >= argc)  WA_error ("need file prefixname after -errts ");
00583           option_data->errts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00584           MTEST (option_data->errts_filename);
00585           strcpy (option_data->errts_filename, argv[nopt]);
00586           nopt++;
00587           continue;
00588         }
00589       
00590 
00591       /*----- unknown command -----*/
00592       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00593       WA_error (message);
00594       
00595     }
00596 
00597   
00598 }
00599 
00600 
00601 
00602 /*---------------------------------------------------------------------------*/
00603 /*
00604   Read time series from specified file name.  This file name may have
00605   a column selector attached.
00606 */
00607 
00608 float * read_time_series 
00609 (
00610   char * ts_filename,          /* time series file name (plus column index) */
00611   int * ts_length              /* output value for time series length */
00612 )
00613 
00614 {
00615   char message[MAX_NAME_LENGTH];   /* error message */
00616   char * cpt;                      /* pointer to column suffix */
00617   char filename[THD_MAX_NAME];     /* time series file name w/o column index */
00618   char subv[THD_MAX_NAME];         /* string containing column index */
00619   MRI_IMAGE * im, * flim;  /* pointers to image structures
00620                               -- used to read 1D ASCII */
00621   float * far;             /* pointer to MRI_IMAGE floating point data */
00622   int nx;                  /* number of time points in time series */
00623   int ny;                  /* number of columns in time series file */
00624   int iy;                  /* time series file column index */
00625   int ipt;                 /* time point index */
00626   float * ts_data = NULL;  /* input time series data */
00627 
00628 
00629   /*----- Read the time series file -----*/
00630   flim = mri_read_1D(ts_filename) ;
00631   if (flim == NULL)
00632     {
00633       sprintf (message,  "Unable to read time series file: %s",  ts_filename);
00634       WA_error (message);
00635     }
00636 
00637   
00638   /*----- Set pointer to data, and set dimensions -----*/
00639   far = MRI_FLOAT_PTR(flim);
00640   nx = flim->nx;
00641   ny = flim->ny; iy = 0 ;
00642   if( ny > 1 ){
00643     fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny);
00644   }
00645 
00646   
00647   /*----- Save the time series data -----*/
00648   *ts_length = nx;
00649   ts_data = (float *) malloc (sizeof(float) * nx);
00650   MTEST (ts_data);
00651   for (ipt = 0;  ipt < nx;  ipt++)
00652     ts_data[ipt] = far[ipt + iy*nx];   
00653   
00654   
00655   mri_free (flim);  flim = NULL;
00656 
00657   return (ts_data);
00658 }
00659 
00660 
00661 /*---------------------------------------------------------------------------*/
00662 /*
00663   Read the input data files.
00664 */
00665 
00666 void read_input_data
00667 (
00668   WA_options * option_data,         /* wavelet program options */
00669   THD_3dim_dataset ** dset_time,    /* input 3d+time data set */
00670   THD_3dim_dataset ** mask_dset,    /* input mask data set */
00671   float ** fmri_data,               /* input fMRI time series data */
00672   int * fmri_length                 /* length of fMRI time series */
00673 )
00674 
00675 {
00676   char message[MAX_NAME_LENGTH];    /* error message */
00677 
00678 
00679   /*----- Read the input fMRI measurement data -----*/
00680   if (option_data->input1D_filename != NULL)
00681     {
00682       /*----- Read the input fMRI 1D time series -----*/
00683       *fmri_data = read_time_series (option_data->input1D_filename, 
00684                                      fmri_length);
00685       if (*fmri_data == NULL)  
00686         { 
00687           sprintf (message,  "Unable to read time series file: %s", 
00688                    option_data->input1D_filename);
00689           WA_error (message);
00690         }  
00691       *dset_time = NULL;
00692     }
00693 
00694   else if (option_data->input_filename != NULL)
00695     {
00696       /*----- Read the input 3d+time dataset -----*/
00697       *dset_time = THD_open_one_dataset (option_data->input_filename);
00698       if (!ISVALID_3DIM_DATASET(*dset_time))  
00699         { 
00700           sprintf (message,  "Unable to open data file: %s", 
00701                    option_data->input_filename);
00702           WA_error (message);
00703         }  
00704       THD_load_datablock ((*dset_time)->dblk);
00705 
00706       if (option_data->mask_filename != NULL)
00707         {
00708           /*----- Read the input mask dataset -----*/
00709           *mask_dset = THD_open_dataset (option_data->mask_filename);
00710           if (!ISVALID_3DIM_DATASET(*mask_dset))  
00711             { 
00712               sprintf (message,  "Unable to open mask file: %s", 
00713                        option_data->mask_filename);
00714               WA_error (message);
00715             }  
00716           THD_load_datablock ((*mask_dset)->dblk);
00717         }
00718     }
00719   else
00720     WA_error ("Must specify input data");
00721 
00722 }
00723 
00724 
00725 /*---------------------------------------------------------------------------*/
00726 /*
00727   Routine to check whether one output file already exists.
00728 */
00729 
00730 void check_one_output_file 
00731 (
00732   THD_3dim_dataset * dset_time,     /* input 3d+time data set */
00733   char * filename                   /* name of output file */
00734 )
00735 
00736 {
00737   char message[MAX_NAME_LENGTH];      /* error message */
00738   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00739   int ierror;                         /* number of errors in editing data */
00740 
00741   
00742   /*----- make an empty copy of input dataset -----*/
00743   new_dset = EDIT_empty_copy( dset_time ) ;
00744   
00745   
00746   ierror = EDIT_dset_items( new_dset ,
00747                             ADN_prefix , filename ,
00748                             ADN_label1 , filename ,
00749                             ADN_self_name , filename ,
00750                             ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 
00751                                                            GEN_FUNC_TYPE ,
00752                             ADN_none ) ;
00753   
00754   if( ierror > 0 )
00755     {
00756       sprintf (message,
00757                "*** %d errors in attempting to create output dataset!\n", 
00758                ierror);
00759       WA_error (message);
00760     }
00761   
00762   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00763     {
00764       sprintf (message,
00765                "Output dataset file %s already exists "
00766                " -- cannot continue!\a\n",
00767                new_dset->dblk->diskptr->header_name);
00768       WA_error (message);
00769     }
00770   
00771   /*----- deallocate memory -----*/   
00772   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00773   
00774 }
00775 
00776 
00777 /*---------------------------------------------------------------------------*/
00778 /*
00779   Routine to check whether output files already exist.
00780 */
00781 
00782 void check_output_files 
00783 (
00784   WA_options * option_data,       /* wavelet analysis program options */
00785   THD_3dim_dataset * dset_time    /* input 3d+time data set */
00786 )
00787 
00788 {
00789 
00790   if (option_data->bucket_filename != NULL)   
00791     check_one_output_file (dset_time, option_data->bucket_filename);
00792 
00793   if (option_data->coefts_filename != NULL)   
00794     check_one_output_file (dset_time, option_data->coefts_filename);
00795 
00796   if (option_data->fitts_filename != NULL)   
00797     check_one_output_file (dset_time, option_data->fitts_filename);
00798 
00799   if (option_data->sgnlts_filename != NULL)   
00800     check_one_output_file (dset_time, option_data->sgnlts_filename);
00801 
00802   if (option_data->errts_filename != NULL)   
00803     check_one_output_file (dset_time, option_data->errts_filename);
00804   
00805 }
00806 
00807 
00808 /*---------------------------------------------------------------------------*/
00809 /*
00810   Routine to check for valid inputs.
00811 */
00812   
00813 void check_for_valid_inputs 
00814 (
00815   WA_options * option_data,       /* wavelet analysis program options */
00816   THD_3dim_dataset * dset_time,   /* input 3d+time data set */
00817   THD_3dim_dataset * mask_dset    /* input mask data set */
00818 )
00819 
00820 {
00821   char message[MAX_NAME_LENGTH];    /* error message */
00822 
00823 
00824   /*----- If mask is used, check for compatible dimensions -----*/
00825   if (mask_dset != NULL)
00826     {
00827       if ( (DSET_NX(dset_time) != DSET_NX(mask_dset))
00828            || (DSET_NY(dset_time) != DSET_NY(mask_dset))
00829            || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) )
00830         {
00831           sprintf (message, "%s and %s have incompatible dimensions",
00832                    option_data->input_filename, option_data->mask_filename);
00833           WA_error (message);
00834         }
00835 
00836       if (DSET_NVALS(mask_dset) != 1 )
00837         WA_error ("Must specify 1 sub-brick from mask dataset");
00838     }
00839 
00840 
00841   /*----- Check whether any of the output files already exist -----*/
00842   check_output_files (option_data, dset_time);
00843 
00844 }
00845 
00846 
00847 /*---------------------------------------------------------------------------*/
00848 /*
00849   Routine to initialize the filters and control variables.
00850 */
00851   
00852 void initialize_filters 
00853 (
00854   WA_options * option_data,       /* wavelet analysis program options */
00855   THD_3dim_dataset * dset_time,   /* input 3d+time data set */
00856   THD_3dim_dataset * mask_dset,   /* input mask data set */
00857   int fmri_length                 /* length of input fMRI time series */
00858 )
00859 
00860 {
00861   char message[MAX_NAME_LENGTH];  /* error message */
00862   int nt;                  /* number of images in input 3d+time dataset */
00863   int NFirst;              /* first image from input 3d+time dataset to use */
00864   int NLast;               /* last image from input 3d+time dataset to use */
00865   int N;                   /* number of usable time points */
00866   int ifilter;             /* filter number index */
00867   int i;                   /* filter coefficient index */
00868   int f, q, p;             /* numbers of model parameters */
00869 
00870 
00871   /*----- Initialize local variables -----*/
00872   if (option_data->input1D_filename != NULL)
00873     nt = fmri_length;
00874   else
00875     nt = DSET_NUM_TIMES (dset_time);
00876 
00877  
00878   /*----- Determine the number of usable time points -----*/
00879   NFirst = option_data->NFirst;
00880 
00881   NLast = option_data->NLast;   
00882   if (NLast > nt-1)  NLast = nt-1;
00883 
00884   N = NLast - NFirst + 1;
00885   N = powerof2(my_log2(N));
00886   NLast = N + NFirst + 1;
00887 
00888   option_data->N = N;
00889   option_data->NLast = NLast;
00890 
00891 
00892   /*----- Initialize the filters -----*/
00893   option_data->stop_filter = 
00894     FWT_1d_stop_filter (option_data->num_stop_filters, option_data->stop_band,
00895                         option_data->stop_mintr,  option_data->stop_maxtr,
00896                         option_data->NFirst, option_data->N);
00897 
00898   option_data->base_filter = 
00899     FWT_1d_pass_filter (option_data->num_base_filters, option_data->base_band,
00900                         option_data->base_mintr,  option_data->base_maxtr,
00901                         option_data->NFirst, option_data->N);
00902                                           
00903   option_data->sgnl_filter = 
00904     FWT_1d_pass_filter (option_data->num_sgnl_filters, option_data->sgnl_band,
00905                         option_data->sgnl_mintr,  option_data->sgnl_maxtr,
00906                         option_data->NFirst, option_data->N);
00907 
00908                                           
00909   /*----- Count numbers of model parameters -----*/
00910   f = 0;
00911   for (i = 0;  i < N;  i++)
00912     if (option_data->stop_filter[i] == 0.0)
00913       {
00914         f++;
00915         option_data->base_filter[i] = 0.0;
00916         option_data->sgnl_filter[i] = 0.0;
00917       }
00918   option_data->f = f;
00919 
00920   q = 0;
00921   for (i = 0;  i < N;  i++)
00922     if (option_data->base_filter[i] == 1.0)
00923       {
00924         q++;
00925         option_data->sgnl_filter[i] = 1.0;
00926       }
00927   option_data->q = q;
00928 
00929   p = 0;
00930   for (i = 0;  i < N;  i++)
00931     if (option_data->sgnl_filter[i] == 1.0)
00932       {
00933         p++;
00934       }
00935   option_data->p = p;    
00936 
00937 
00938   return;
00939 }
00940 
00941 
00942 /*---------------------------------------------------------------------------*/
00943 /*
00944   Allocate memory for output volumes.
00945 */
00946 
00947 void allocate_memory 
00948 (
00949   WA_options * option_data,         /* wavelet analysis options */
00950   THD_3dim_dataset * dset,          /* input 3d+time data set */
00951 
00952   float *** coef_vol,       /* array of volumes of signal model parameters */
00953   float ** mse_vol,         /* volume of full model mean square error */
00954   float ** ffull_vol,       /* volume of full model F-statistics */
00955   float ** rfull_vol,       /* volume of full model R^2 stats. */
00956 
00957   float *** coefts_vol,     /* volumes for forward wavelet transform coefs. */
00958   float *** fitts_vol,      /* volumes for wavelet filtered time series data */
00959   float *** sgnlts_vol,     /* volumes for signal model fit to input data */
00960   float *** errts_vol       /* volumes for residual errors */
00961 )
00962 
00963 {
00964   int nxyz;                  /* total number of voxels */
00965   int ixyz;                  /* voxel index */
00966   int it;                    /* time point index */
00967   int N;                     /* number of usable data points */
00968   int ip;                    /* parameter index */
00969   int p;                     /* number of parameters in the full model */
00970 
00971 
00972   /*----- Initialize local variables -----*/
00973   nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
00974   N = option_data->N;
00975   p = option_data->p;
00976 
00977 
00978   /*----- Allocate memory space for volume data -----*/
00979   *coef_vol  = (float **) malloc (sizeof(float *) * p);   MTEST(*coef_vol);
00980 
00981   for (ip = 0;  ip < p;  ip++)
00982     {
00983       (*coef_vol)[ip]  = (float *) malloc (sizeof(float) * nxyz);
00984       MTEST((*coef_vol)[ip]);    
00985        for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00986         {
00987           (*coef_vol)[ip][ixyz]  = 0.0;
00988         }
00989     }
00990 
00991 
00992   *mse_vol   = (float *) malloc (sizeof(float) * nxyz);   MTEST(*mse_vol);
00993   *ffull_vol = (float *) malloc (sizeof(float) * nxyz);   MTEST(*ffull_vol);
00994   *rfull_vol = (float *) malloc (sizeof(float) * nxyz);   MTEST(*rfull_vol);
00995   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00996     {
00997       (*mse_vol)[ixyz]   = 0.0;
00998       (*ffull_vol)[ixyz] = 0.0;
00999       (*rfull_vol)[ixyz] = 0.0;
01000     }
01001 
01002 
01003   /*----- Allocate memory for forward wavelet transform coefficients -----*/
01004   if (option_data->coefts_filename != NULL)
01005     {
01006       *coefts_vol = (float **) malloc (sizeof(float **) * N);
01007       MTEST (*coefts_vol);
01008       for (it = 0;  it < N;  it++)
01009         {
01010           (*coefts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01011           MTEST ((*coefts_vol)[it]);
01012           for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01013             (*coefts_vol)[it][ixyz] = 0.0;
01014         }
01015     }
01016 
01017   /*----- Allocate memory for filtered time series -----*/
01018   if (option_data->fitts_filename != NULL)
01019     {
01020       *fitts_vol = (float **) malloc (sizeof(float **) * N);
01021       MTEST (*fitts_vol);
01022       for (it = 0;  it < N;  it++)
01023         {
01024           (*fitts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01025           MTEST ((*fitts_vol)[it]);
01026           for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01027             (*fitts_vol)[it][ixyz] = 0.0;
01028         }
01029     }
01030 
01031   /*----- Allocate memory for signal model time series -----*/
01032   if (option_data->sgnlts_filename != NULL)
01033     {
01034       *sgnlts_vol = (float **) malloc (sizeof(float **) * N);
01035       MTEST (*sgnlts_vol);
01036       for (it = 0;  it < N;  it++)
01037         {
01038           (*sgnlts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01039           MTEST ((*sgnlts_vol)[it]);
01040           for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01041             (*sgnlts_vol)[it][ixyz] = 0.0;
01042         }
01043     }
01044 
01045   /*----- Allocate memory for residual time series -----*/
01046   if (option_data->errts_filename != NULL)
01047     {
01048       *errts_vol = (float **) malloc (sizeof(float **) * N);
01049       MTEST (*errts_vol);
01050       for (it = 0;  it < N;  it++)
01051         {
01052           (*errts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz);
01053           MTEST ((*errts_vol)[it]);
01054           for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01055             (*errts_vol)[it][ixyz] = 0.0;
01056         }
01057     }
01058 }
01059 
01060 
01061 /*---------------------------------------------------------------------------*/
01062 /*
01063   Perform all program initialization.
01064 */
01065 
01066 void initialize_program 
01067 (
01068   int argc,                         /* number of input arguments */
01069   char ** argv,                     /* array of input arguments */ 
01070   WA_options ** option_data,        /* wavelet analysis options */
01071   THD_3dim_dataset ** dset_time,    /* input 3d+time data set */
01072   THD_3dim_dataset ** mask_dset,    /* input mask data set */
01073   float ** fmri_data,               /* input fMRI time series data */
01074   int * fmri_length,                /* length of fMRI time series */
01075 
01076   float *** coef_vol,        /* array of volumes of signal model parameters */
01077   float ** mse_vol,          /* volume of full model mean square error */
01078   float ** ffull_vol,        /* volume of full model F-statistics */
01079   float ** rfull_vol,        /* volume of full model R^2 stats. */
01080 
01081   float *** coefts_vol,      /* volumes for forward wavelet transform coefs. */
01082   float *** fitts_vol,       /* volumes for wavelet filtered time series  */
01083   float *** sgnlts_vol,      /* volumes for signal model fit to input data */
01084   float *** errts_vol        /* volumes for residual errors */
01085 )
01086      
01087 {
01088 
01089   /*----- Allocate memory -----*/
01090   *option_data = (WA_options *) malloc (sizeof(WA_options));
01091 
01092    
01093   /*----- Get command line inputs -----*/
01094   get_options (argc, argv, *option_data);
01095 
01096 
01097   /*----- Read input data -----*/
01098   read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length);
01099 
01100 
01101   /*----- Check for valid inputs -----*/
01102   check_for_valid_inputs (*option_data, *dset_time, *mask_dset);
01103   
01104 
01105   /*----- Initialize the filters and control variables  -----*/
01106   initialize_filters (*option_data, *dset_time, *mask_dset, *fmri_length);
01107   
01108 
01109   /*----- Allocate memory for output volumes -----*/
01110   if ((*option_data)->input1D_filename == NULL)
01111     allocate_memory (*option_data, *dset_time, 
01112                      coef_vol, mse_vol, ffull_vol, rfull_vol,
01113                      coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
01114 
01115 }
01116 
01117 
01118 /*---------------------------------------------------------------------------*/
01119 /*
01120   Get the time series for one voxel from the AFNI 3d+time data set.
01121 */
01122 
01123 void extract_ts_array 
01124 (
01125   THD_3dim_dataset * dset_time,      /* input 3d+time dataset */
01126   int iv,                            /* get time series for this voxel */
01127   float * ts_array                   /* time series data for voxel #iv */
01128 )
01129 
01130 {
01131   MRI_IMAGE * im = NULL;   /* intermediate float data */
01132   float * ar = NULL;       /* pointer to float data */
01133   int ts_length;           /* length of input 3d+time data set */
01134   int it;                  /* time index */
01135 
01136 
01137   /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/
01138   im = THD_extract_series (iv, dset_time, 0);
01139 
01140 
01141   /*----- Verify extraction -----*/
01142   if (im == NULL)  WA_error ("Unable to extract data from 3d+time dataset");
01143 
01144 
01145   /*----- Now extract time series from MRI_IMAGE -----*/
01146   ts_length = DSET_NUM_TIMES (dset_time);
01147   ar = MRI_FLOAT_PTR (im);
01148   for (it = 0;  it < ts_length;  it++)
01149     {
01150       ts_array[it] = ar[it];
01151     }
01152 
01153 
01154   /*----- Release memory -----*/
01155   mri_free (im);   im = NULL;
01156 
01157 }
01158 
01159 
01160 /*---------------------------------------------------------------------------*/
01161 /*
01162   Save results for this voxel.
01163 */
01164 
01165 void save_voxel 
01166 (
01167   WA_options * option_data,    /* wavelet analysis options */
01168   int iv,                      /* current voxel index */      
01169   float * coef,                /* regression parameters */
01170   float mse,                   /* full model mean square error */
01171   float ffull,                 /* full model F-statistic */
01172   float rfull,                 /* full model R^2 stat. */
01173 
01174   float * coefts,             /* forward wavelet transform coefficients */
01175   float * fitts,              /* wavelet filtered time series */
01176   float * sgnlts,              /* signal model fitted time series */
01177   float * errts,               /* signal model residual error time series */
01178 
01179   float ** coef_vol,        /* array of volumes of signal model parameters */
01180   float * mse_vol,          /* volume of full model mean square error */
01181   float * ffull_vol,        /* volume of full model F-statistics */
01182   float * rfull_vol,        /* volume of full model R^2 stats. */
01183 
01184   float ** coefts_vol,     /* volumes for forward wavelet transform coefs. */
01185   float ** fitts_vol,      /* volumes for wavelet filtered time series data */
01186   float ** sgnlts_vol,      /* volumes for signal model fit to input data */
01187   float ** errts_vol        /* volumes for residual errors */
01188 
01189 )
01190 
01191 {
01192   int ip;                   /* parameter index */ 
01193   int p;                    /* total number of parameters */
01194   int it;                   /* time point index */
01195   int N;                    /* number of usable data points */
01196 
01197 
01198   /*----- Initialize local variables -----*/
01199   p = option_data->p;
01200   N = option_data->N;
01201 
01202 
01203   /*----- Saved wavelet coefficients -----*/
01204   for (ip = 0;  ip < p;  ip++)
01205     {
01206       coef_vol[ip][iv]  = coef[ip];
01207     }
01208 
01209 
01210   /*----- Save full model mean square error -----*/
01211   mse_vol[iv] = mse;
01212 
01213 
01214   /*----- Save regression F-statistic -----*/
01215   ffull_vol[iv] = ffull;
01216 
01217 
01218   /*----- Save R^2 values -----*/
01219   rfull_vol[iv] = rfull;
01220 
01221 
01222   /*----- Save the forward wavelet transform coefficients -----*/
01223   if (coefts_vol != NULL)
01224     for (it = 0;  it < N;  it++)
01225       coefts_vol[it][iv] = coefts[it];
01226 
01227 
01228   /*----- Save the wavelet filtered time series -----*/
01229   if (fitts_vol != NULL)
01230     for (it = 0;  it < N;  it++)
01231       fitts_vol[it][iv] = fitts[it];
01232 
01233 
01234   /*----- Save the fitted signal model time series -----*/
01235   if (sgnlts_vol != NULL)
01236     for (it = 0;  it < N;  it++)
01237       sgnlts_vol[it][iv] = sgnlts[it];
01238 
01239 
01240   /*----- Save the residual errors time series -----*/
01241   if (errts_vol != NULL)
01242     for (it = 0;  it < N;  it++)
01243       errts_vol[it][iv] = errts[it];
01244 
01245 
01246 }
01247 
01248 
01249 /*---------------------------------------------------------------------------*/
01250 /*
01251   Calculate the wavelet signal model and associated statistics.
01252 */
01253 
01254 void calculate_results 
01255 (
01256   WA_options * option_data,         /* wavelet analysis options */
01257   THD_3dim_dataset * dset,          /* input 3d+time data set */
01258   THD_3dim_dataset * mask,          /* input mask data set */
01259   float * fmri_data,                /* input fMRI time series data */
01260   int fmri_length,                  /* length of fMRI time series */
01261 
01262   float ** coef_vol,        /* array of volumes of signal model parameters */
01263   float * mse_vol,          /* volume of full model mean square error */
01264   float * ffull_vol,        /* volume of F-statistic for the full model */
01265   float * rfull_vol,        /* volume of R^2 for the full model */
01266 
01267   float ** coefts_vol,     /* volumes for forward wavelet transform coefs. */
01268   float ** fitts_vol,      /* volumes for wavelet filtered time series */
01269   float ** sgnlts_vol,      /* volumes for full model fit to input data */
01270   float ** errts_vol        /* volumes for residual errors */
01271 )
01272   
01273 {
01274   float * ts_array = NULL;    /* array of measured data for one voxel */
01275   float mask_val[1];          /* value of mask at current voxel */
01276 
01277   int f;
01278   int p;                      /* number of parameters in the full model */
01279   int q;                      /* number of parameters in the baseline model */
01280   int m;                      /* parameter index */
01281   int n;                      /* data point index */
01282 
01283   float * coef = NULL;        /* regression parameters */
01284   float sse_base;             /* baseline model error sum of squares */ 
01285   float sse_full;             /* full model error sum of squares */
01286   float ffull;                /* full model F-statistic */
01287   float rfull;                /* full model R^2 stat. */
01288 
01289   int ixyz;                   /* voxel index */
01290   int nxyz;                   /* number of voxels per dataset */
01291 
01292   int nt;                  /* number of images in input 3d+time dataset */
01293   int NFirst;              /* first image from input 3d+time dataset to use */
01294   int NLast;               /* last image from input 3d+time dataset to use */
01295   int N;                   /* number of usable data points */
01296 
01297   int i;                   /* data point index */
01298 
01299   float * coefts = NULL;   /* filtered FWT coefficients */
01300   float * fitts  = NULL;   /* filterd time series */
01301   float * sgnlts = NULL;   /* signal model fitted time series */
01302   float * errts  = NULL;   /* residual error time series */
01303 
01304   char * label = NULL;     /* string containing stat. summary of results */
01305 
01306 
01307   /*----- Initialize local variables -----*/
01308   if (option_data->input1D_filename != NULL)
01309     {
01310       nxyz = 1;
01311       nt = fmri_length;
01312     }
01313   else
01314     {
01315       nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;       
01316       nt = DSET_NUM_TIMES (dset);
01317     }
01318 
01319   NFirst = option_data->NFirst;
01320   NLast = option_data->NLast;   
01321   N = option_data->N;
01322   f = option_data->f;
01323   q = option_data->q;
01324   p = option_data->p;
01325 
01326 
01327   /*----- Allocate memory for arrays -----*/
01328   ts_array = (float *) malloc (sizeof(float) * nt);   MTEST (ts_array);
01329   coef     = (float *) malloc (sizeof(float) * p);    MTEST (coef);
01330   coefts  = (float *) malloc (sizeof(float) * N);    MTEST (coefts);
01331   fitts   = (float *) malloc (sizeof(float) * N);    MTEST (fitts);
01332   sgnlts   = (float *) malloc (sizeof(float) * N);    MTEST (sgnlts);
01333   errts    = (float *) malloc (sizeof(float) * N);    MTEST (errts);
01334 
01335   
01336   /*----- Loop over all voxels -----*/
01337   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01338     {
01339       /*----- Apply mask? -----*/
01340       if (mask != NULL)
01341         {
01342           extract_ts_array (mask, ixyz, mask_val);
01343           if (mask_val[0] == 0.0)  continue; 
01344         }
01345       
01346       
01347       /*----- Extract Y-data for this voxel -----*/
01348       if (option_data->input1D_filename != NULL)
01349         {
01350           for (i = 0;  i < N;  i++)
01351             ts_array[i] = fmri_data[i+NFirst];
01352         }
01353       else
01354         extract_ts_array (dset, ixyz, ts_array);
01355 
01356           
01357       /*----- Perform the wavelet analysis for this voxel-----*/
01358       wavelet_analysis (option_data->wavelet_type, 
01359                         f, option_data->stop_filter,
01360                         q, option_data->base_filter,
01361                         p, option_data->sgnl_filter,
01362                         N, ts_array+NFirst, coef,
01363                         &sse_base, &sse_full, &ffull, &rfull,
01364                         coefts, fitts, sgnlts, errts);
01365 
01366       
01367       /*----- Save results for this voxel -----*/
01368       if (option_data->input1D_filename == NULL)
01369         save_voxel (option_data, ixyz, coef, sse_full/(N-f-p), ffull, rfull, 
01370                     coefts, fitts, sgnlts, errts, 
01371                     coef_vol, mse_vol, ffull_vol, rfull_vol, 
01372                     coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
01373       else
01374         {  
01375           /*----- If only one voxel, save results as .1D files -----*/
01376           ts_fprint ("WA.coefts.1D", N, coefts);
01377           ts_fprint ("WA.fitts.1D",  N, fitts);
01378           ts_fprint ("WA.sgnlts.1D", N, sgnlts);
01379           ts_fprint ("WA.errts.1D",  N, errts);
01380         }
01381       
01382           
01383       /*----- Report results for this voxel -----*/
01384       if ( ((ffull > option_data->fdisp) && (option_data->fdisp >= 0.0))
01385            || (option_data->input1D_filename != NULL) )
01386         {
01387           printf ("\n\nResults for Voxel #%d: \n", ixyz);
01388           report_results (N, NFirst, f, p, q, 
01389                           option_data->base_filter, option_data->sgnl_filter,
01390                           coef, sse_base, sse_full, ffull, rfull, 
01391                           &label);
01392           printf ("%s \n", label);
01393         }
01394 
01395           
01396     }  /*----- Loop over voxels -----*/
01397 
01398 
01399 
01400   /*----- Release memory -----*/
01401   free (ts_array);  ts_array = NULL;
01402   free (coef);      coef     = NULL;
01403   free (coefts);    coefts   = NULL;
01404   free (fitts);     fitts    = NULL;
01405   free (sgnlts);    sgnlts   = NULL;
01406   free (errts);     errts    = NULL;
01407 
01408 }
01409 
01410 
01411 /*---------------------------------------------------------------------------*/
01412 /*
01413   Convert one volume to another type, autoscaling:
01414      nxy   = # voxels
01415      itype = input datum type
01416      ivol  = pointer to input volume
01417      otype = output datum type
01418      ovol  = pointer to output volume (again, must be pre-malloc-ed)
01419   Return value is the scaling factor used (0.0 --> no scaling).
01420 */
01421 
01422 float EDIT_coerce_autoscale_new( int nxyz ,
01423                                  int itype,void *ivol , int otype,void *ovol )
01424 {
01425   float fac=0.0 , top ;
01426   
01427   if( MRI_IS_INT_TYPE(otype) ){
01428     top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01429     if (top == 0.0)  fac = 0.0;
01430     else  fac = MRI_TYPE_maxval[otype]/top;
01431   }
01432   
01433   EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01434   return ( fac );
01435 }
01436 
01437 
01438 /*---------------------------------------------------------------------------*/
01439 /*
01440   Routine to write one AFNI 3d+time data set. 
01441 */
01442 
01443 
01444 void write_ts_array 
01445 (
01446   int argc,                              /* number of input arguments */
01447   char ** argv,                          /* array of input arguments */ 
01448   WA_options * option_data,              /* wavelet analysis options */
01449   int ts_length,                         /* length of time series data */  
01450   float ** vol_array,                    /* output time series volume data */
01451   char * output_filename                 /* output afni data set file name */
01452 )
01453 
01454 {
01455   const float EPSILON = 1.0e-10;
01456 
01457   THD_3dim_dataset * dset = NULL;        /* input afni data set pointer */
01458   THD_3dim_dataset * new_dset = NULL;    /* output afni data set pointer */
01459   int ib;                                /* sub-brick index */ 
01460   int ierror;                            /* number of errors in editing data */
01461   int nxyz;                              /* total number of voxels */ 
01462   float factor;             /* factor is new scale factor for sub-brick #ib */
01463   char * input_filename;    /* input afni data set file name */
01464   short ** bar = NULL;      /* bar[ib] points to data for sub-brick #ib */
01465   float * fbuf;             /* float buffer */
01466   float * volume;           /* pointer to volume of data */
01467   char label[80];           /* label for output file */ 
01468   
01469 
01470   /*----- Initialize local variables -----*/
01471   input_filename = option_data->input_filename;
01472   dset = THD_open_one_dataset (input_filename);
01473   nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
01474 
01475  
01476   /*----- allocate memory -----*/
01477   bar  = (short **) malloc (sizeof(short *) * ts_length);   MTEST (bar);
01478   fbuf = (float *)  malloc (sizeof(float)   * ts_length);   MTEST (fbuf);
01479   
01480   
01481   /*-- make an empty copy of the prototype dataset, for eventual output --*/
01482   new_dset = EDIT_empty_copy (dset);
01483 
01484 
01485   /*----- Record history of dataset -----*/
01486   tross_Copy_History( dset , new_dset ) ;
01487 
01488   { char * commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
01489     sprintf (label, "Output prefix: %s", output_filename);
01490     if( commandline != NULL )
01491        tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
01492     else
01493        tross_Append_History ( new_dset, label);
01494     free(commandline) ;
01495   }
01496 
01497   /*----- Delete prototype dataset -----*/
01498   THD_delete_3dim_dataset (dset, False);  dset = NULL ;
01499   
01500 
01501   ierror = EDIT_dset_items (new_dset,
01502                             ADN_prefix,      output_filename,
01503                             ADN_label1,      output_filename,
01504                             ADN_self_name,   output_filename,
01505                             ADN_malloc_type, DATABLOCK_MEM_MALLOC,  
01506                             ADN_datum_all,   MRI_short,   
01507                             ADN_nvals,       ts_length,
01508                             ADN_ntt,         ts_length,
01509                             ADN_none);
01510   
01511   if( ierror > 0 ){
01512     fprintf(stderr,
01513           "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01514     exit(1) ;
01515   }
01516   
01517   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01518     fprintf(stderr,
01519             "*** Output dataset file %s already exists--cannot continue!\a\n",
01520             new_dset->dblk->diskptr->header_name ) ;
01521     exit(1) ;
01522   }
01523 
01524   
01525   /*----- attach bricks to new data set -----*/
01526   for (ib = 0;  ib < ts_length;  ib++)
01527     {
01528 
01529       /*----- Set pointer to appropriate volume -----*/
01530       volume = vol_array[ib];
01531       
01532       /*----- Allocate memory for output sub-brick -----*/
01533       bar[ib]  = (short *) malloc (sizeof(short) * nxyz);
01534       MTEST (bar[ib]);
01535 
01536       /*----- Convert data type to short for this sub-brick -----*/
01537       factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01538                                           MRI_short, bar[ib]);
01539       if (factor < EPSILON)  factor = 0.0;
01540       else factor = 1.0 / factor;
01541       fbuf[ib] = factor;
01542 
01543       /*----- attach bar[ib] to be sub-brick #ib -----*/
01544       mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib)); 
01545     }
01546 
01547 
01548   /*----- write afni data set -----*/
01549 
01550   (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01551 
01552   THD_load_statistics (new_dset);
01553   THD_write_3dim_dataset (NULL, NULL, new_dset, True);
01554   fprintf(stderr,"--- Wrote 3d+time dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
01555   
01556 
01557   /*----- deallocate memory -----*/   
01558   THD_delete_3dim_dataset (new_dset, False);   new_dset = NULL ;
01559   free (fbuf);   fbuf = NULL;
01560 
01561 }
01562 
01563 
01564 /*---------------------------------------------------------------------------*/
01565 /*
01566   Attach one sub-brick to output bucket data set.
01567 */
01568 
01569 void attach_sub_brick
01570 (
01571   THD_3dim_dataset * new_dset,      /* output bucket dataset */
01572   int ibrick,               /* sub-brick indices */
01573   float * volume,           /* volume of floating point data */
01574   int nxyz,                 /* total number of voxels */
01575   int brick_type,           /* indicates statistical type of sub-brick */
01576   char * brick_label,       /* character string label for sub-brick */
01577   int dof, 
01578   int ndof, 
01579   int ddof,                 /* degrees of freedom */
01580   short ** bar              /* bar[ib] points to data for sub-brick #ib */  
01581 )
01582 
01583 {
01584   const float EPSILON = 1.0e-10;
01585   float factor;             /* factor is new scale factor for sub-brick #ib */
01586 
01587 
01588   /*----- allocate memory for output sub-brick -----*/
01589   bar[ibrick]  = (short *) malloc (sizeof(short) * nxyz);
01590   MTEST (bar[ibrick]);
01591   factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01592                                       MRI_short, bar[ibrick]);
01593   
01594   if (factor < EPSILON)  factor = 0.0;
01595   else factor = 1.0 / factor;
01596   
01597 
01598   /*----- edit the sub-brick -----*/
01599   EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
01600   EDIT_BRICK_FACTOR (new_dset, ibrick, factor);
01601 
01602   if (brick_type == FUNC_TT_TYPE)
01603     EDIT_BRICK_TO_FITT (new_dset, ibrick, dof);
01604   else if (brick_type == FUNC_FT_TYPE)
01605     EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof);
01606   
01607   
01608   /*----- attach bar[ib] to be sub-brick #ibrick -----*/
01609   EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]);
01610 
01611 }
01612 
01613 /*---------------------------------------------------------------------------*/
01614 /*
01615   Routine to write one bucket data set.
01616 */
01617 
01618 void write_bucket_data 
01619 (
01620   int argc,                         /* number of input arguments */
01621   char ** argv,                     /* array of input arguments */ 
01622   WA_options * option_data,         /* wavelet analysis options */
01623 
01624   float ** coef_vol,       /* array of volumes of signal model parameters */
01625   float * mse_vol,         /* volume of full model mean square error */
01626   float * ffull_vol,       /* volume of full model F-statistics */
01627   float * rfull_vol        /* volume of full model R^2 statistics */
01628 )
01629 
01630 {
01631   THD_3dim_dataset * old_dset = NULL;      /* prototype dataset */
01632   THD_3dim_dataset * new_dset = NULL;      /* output bucket dataset */
01633   char output_prefix[MAX_NAME_LENGTH];     /* prefix name for bucket dataset */
01634   char output_session[MAX_NAME_LENGTH];    /* directory for bucket dataset */
01635   int nbricks;              /* number of sub-bricks in bucket dataset */
01636   short ** bar = NULL;      /* bar[ib] points to data for sub-brick #ib */
01637 
01638   int brick_type;           /* indicates statistical type of sub-brick */
01639   int brick_coef;           /* regression coefficient index for sub-brick */
01640   char brick_label[MAX_NAME_LENGTH]; /* character string label for sub-brick */
01641 
01642   int ierror;               /* number of errors in editing data */
01643   float * volume;           /* volume of floating point data */
01644 
01645   int N;                    /* number of usable data points */
01646   int NFirst;
01647   int f;                    /* number of parameters removed by stop filter */
01648   int p;                    /* number of parameters in the signal model */
01649   int q;                    /* number of parameters in the rdcd model */
01650   int nxyz;                 /* total number of voxels */
01651   int nt;                   /* number of images in input 3d+time dataset */
01652   int it;
01653   int ilag;                 /* lag index */
01654   int icoef;                /* coefficient index */
01655   int ibrick;               /* sub-brick index */
01656   int dof, ndof, ddof;      /* degrees of freedom */
01657   char label[MAX_NAME_LENGTH];   /* general label for sub-bricks */
01658 
01659 
01660   /*----- read prototype dataset -----*/
01661   old_dset = THD_open_one_dataset (option_data->input_filename);
01662 
01663     
01664   /*----- Initialize local variables -----*/
01665   nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz;  
01666   nt = DSET_NUM_TIMES (old_dset);
01667 
01668 
01669   f = option_data->f;
01670   q = option_data->q;
01671   p = option_data->p;
01672   N = option_data->N;
01673   NFirst = option_data->NFirst;
01674   printf ("f = %d   q = %d   p = %d   N = %d \n", f, q, p, N);
01675   
01676 
01677   /*----- Calculate number of sub-bricks in the bucket -----*/
01678   nbricks = p * option_data->cout
01679     + option_data->rout + option_data->fout + option_data->vout;
01680   printf ("nbricks = %d \n", nbricks);
01681 
01682 
01683   strcpy (output_prefix, option_data->bucket_filename);
01684   strcpy (output_session, "./");
01685   
01686  
01687   /*----- allocate memory -----*/
01688   bar  = (short **) malloc (sizeof(short *) * nbricks);
01689   MTEST (bar);
01690   
01691 
01692   /*-- make an empty copy of prototype dataset, for eventual output --*/
01693   new_dset = EDIT_empty_copy (old_dset);
01694 
01695 
01696   /*----- Record history of dataset -----*/
01697   tross_Copy_History( old_dset , new_dset ) ;
01698   tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
01699   sprintf (label, "Output prefix: %s", output_prefix);
01700   tross_Append_History ( new_dset, label);
01701 
01702   
01703   /*----- delete prototype dataset -----*/ 
01704   THD_delete_3dim_dataset( old_dset , False );  old_dset = NULL ;
01705 
01706 
01707   /*----- Modify some structural properties.  Note that the nbricks
01708           just make empty sub-bricks, without any data attached. -----*/
01709   ierror = EDIT_dset_items (new_dset,
01710                             ADN_prefix,          output_prefix,
01711                             ADN_directory_name,  output_session,
01712                             ADN_type,            HEAD_FUNC_TYPE,
01713                             ADN_func_type,       FUNC_BUCK_TYPE,
01714                             ADN_datum_all,       MRI_short ,   
01715                             ADN_ntt,             0,               /* no time */
01716                             ADN_nvals,           nbricks,
01717                             ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,  
01718                             ADN_none ) ;
01719   
01720   if( ierror > 0 )
01721     {
01722       fprintf(stderr, 
01723               "*** %d errors in attempting to create bucket dataset!\n", 
01724               ierror);
01725       exit(1);
01726     }
01727   
01728   if (THD_is_file(DSET_HEADNAME(new_dset))) 
01729     {
01730       fprintf(stderr,
01731               "*** Output dataset file %s already exists--cannot continue!\n",
01732               DSET_HEADNAME(new_dset));
01733       exit(1);
01734     }
01735   
01736 
01737   /*----- Attach individual sub-bricks to the bucket dataset -----*/
01738   ibrick = 0;
01739 
01740 
01741   /*----- Full model wavelet coefficients -----*/
01742   if (option_data->cout)
01743     {
01744       /*----- User can choose to place full model stats first -----*/
01745       if (option_data->stat_first)  
01746         ibrick += option_data->vout + option_data->rout + option_data->fout;
01747 
01748       icoef = 0;
01749       for (it = 0;  it < N;  it++)        
01750         {                      
01751           if (option_data->sgnl_filter[it])
01752             {                             
01753               /*----- Wavelet coefficient -----*/
01754               brick_type = FUNC_FIM_TYPE;
01755               
01756               {
01757                 int band, mintr, maxtr;
01758                 
01759                 if (it == 0)
01760                   {
01761                     band = -1;
01762                     mintr = 0;
01763                     maxtr = N-1;
01764                   }
01765                 else
01766                   {
01767                     band = my_log2(it);
01768                     mintr = (it - powerof2(band)) * powerof2(my_log2(N)-band);
01769                     maxtr = mintr + powerof2(my_log2(N)-band) - 1;
01770                   }
01771             
01772                 mintr += NFirst;
01773                 maxtr += NFirst;
01774             
01775                 if (option_data->base_filter[it])
01776                   sprintf (brick_label, "B(%2d)[%3d,%3d]", band, mintr, maxtr);
01777                 else
01778                   sprintf (brick_label, "S(%2d)[%3d,%3d]", band, mintr, maxtr);
01779               }
01780 
01781               volume = coef_vol[icoef];           
01782               attach_sub_brick (new_dset, ibrick, volume, nxyz, 
01783                                 brick_type, brick_label, 0, 0, 0, bar);
01784               
01785               ibrick++;
01786               icoef++;
01787             }
01788           
01789         }  /* End loop over Wavelet coefficients */
01790     }  /* if (option_data->cout) */
01791 
01792 
01793       
01794   /*----- Statistics for full model -----*/
01795   if (option_data->stat_first)  ibrick = 0;
01796 
01797   /*----- Full model MSE -----*/
01798   if (option_data->vout)
01799     {
01800       brick_type = FUNC_FIM_TYPE;
01801       sprintf (brick_label, "Full MSE");
01802       volume = mse_vol;
01803       attach_sub_brick (new_dset, ibrick, volume, nxyz, 
01804                         brick_type, brick_label, 0, 0, 0, bar);
01805       ibrick++;
01806     }
01807 
01808   /*----- Full model R^2 -----*/
01809   if (option_data->rout)
01810     {
01811       brick_type = FUNC_THR_TYPE;
01812       sprintf (brick_label, "Full R^2");
01813       volume = rfull_vol;
01814       attach_sub_brick (new_dset, ibrick, volume, nxyz, 
01815                         brick_type, brick_label, 0, 0, 0, bar);
01816       ibrick++;
01817     }
01818 
01819   /*----- Full model F-stat -----*/
01820   if (option_data->fout)
01821     {
01822       brick_type = FUNC_FT_TYPE;
01823       ndof = p - q;
01824       ddof = N - f - p;
01825       sprintf (brick_label, "Full F-stat");
01826       volume = ffull_vol;
01827       attach_sub_brick (new_dset, ibrick, volume, nxyz, 
01828                         brick_type, brick_label, 0, ndof, ddof, bar);
01829       ibrick++;
01830     }
01831 
01832 
01833   /*----- write bucket data set -----*/
01834   THD_load_statistics (new_dset);
01835   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01836   fprintf(stderr,"Writing `bucket' dataset ");
01837   fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01838 
01839   
01840   /*----- deallocate memory -----*/   
01841   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01842 
01843 }
01844 
01845 
01846 /*---------------------------------------------------------------------------*/
01847 /*
01848   Write out the user requested output files.
01849 */
01850 
01851 void output_results
01852 (
01853   int argc,                         /* number of input arguments */
01854   char ** argv,                     /* array of input arguments */ 
01855   WA_options * option_data,         /* wavelet analysis options */
01856 
01857   float ** coef_vol,        /* array of volumes of signal model parameters */
01858   float * mse_vol,          /* volume of full model mean square error */
01859   float * ffull_vol,        /* volume of full model F-statistics */
01860   float * rfull_vol,        /* volume of full model R^2 statistics */
01861   float ** coefts_vol,     /* volumes for forward wavelet transform coefs. */
01862   float ** fitts_vol,      /* volumes for wavelet filtered time series */
01863   float ** sgnlts_vol,      /* volumes for signal model time series */
01864   float ** errts_vol        /* volumes for residual errors */
01865 )
01866 
01867 {
01868   int N;                    /* number of usable data points */
01869 
01870 
01871   /*----- Initialize local variables -----*/
01872   N = option_data->N;
01873 
01874 
01875   /*----- Write the bucket dataset -----*/
01876   if (option_data->bucket_filename != NULL)
01877     write_bucket_data (argc, argv, option_data,  coef_vol, 
01878                        mse_vol, ffull_vol, rfull_vol);
01879 
01880 
01881 
01882   /*----- Write the forward wavelet transform coefficients -----*/
01883   if (option_data->coefts_filename != NULL)
01884     write_ts_array (argc, argv, option_data, N, coefts_vol, 
01885                     option_data->coefts_filename);
01886 
01887 
01888   /*----- Write the wavelet filtered 3d+time dataset -----*/
01889   if (option_data->fitts_filename != NULL)
01890     write_ts_array (argc, argv, option_data, N, fitts_vol, 
01891                     option_data->fitts_filename);
01892 
01893 
01894   /*----- Write the signal model 3d+time dataset -----*/
01895   if (option_data->sgnlts_filename != NULL)
01896     write_ts_array (argc, argv, option_data, N, sgnlts_vol, 
01897                     option_data->sgnlts_filename);
01898 
01899 
01900   /*----- Write the residual errors 3d+time dataset -----*/
01901   if (option_data->errts_filename != NULL)
01902     write_ts_array (argc, argv, option_data, N, errts_vol, 
01903                     option_data->errts_filename);
01904 
01905 
01906 
01907 }
01908 
01909 
01910 /*---------------------------------------------------------------------------*/
01911 
01912 void terminate_program
01913 (
01914   WA_options ** option_data,         /* wavelet analysis options */
01915   float ** fmri_data,                /* input fMRI time series data */
01916 
01917   float *** coef_vol,       /* array of volumes of signal model parameters */
01918   float ** mse_vol,         /* volume of full model mean square error */
01919   float ** ffull_vol,       /* volume of full model F-statistics */
01920   float ** rfull_vol,       /* volume of full model R^2 statistics */
01921 
01922   float *** coefts_vol,     /* volumes for forward wavelet transform coefs. */
01923   float *** fitts_vol,      /* volumes for wavelet filtered time series data */
01924   float *** sgnlts_vol,     /* volumes for signal model fit to input data */
01925   float *** errts_vol       /* volumes for residual errors */
01926 )
01927 
01928 {
01929   int p;                    /* number of parameters in full model */
01930   int ip;                   /* parameter index */
01931   int is;                   /* stimulus index */
01932   int it;                   /* time index */
01933   int N;                    /* number of usable data points */
01934 
01935 
01936   /*----- Initialize local variables -----*/
01937   p = (*option_data)->p;
01938   N = (*option_data)->N;
01939 
01940 
01941   /*----- Deallocate memory for option data -----*/   
01942   if ((*option_data)->stop_filter != NULL)  free ((*option_data)->stop_filter);
01943   if ((*option_data)->base_filter != NULL)  free ((*option_data)->base_filter);
01944   if ((*option_data)->sgnl_filter != NULL)  free ((*option_data)->sgnl_filter);
01945   if ((*option_data)->input_filename != NULL) 
01946     free ((*option_data)->input_filename);
01947   if ((*option_data)->mask_filename != NULL)  
01948     free ((*option_data)->mask_filename);  
01949   if ((*option_data)->input1D_filename != NULL) 
01950     free ((*option_data)->input1D_filename);
01951   if ((*option_data)->bucket_filename != NULL) 
01952     free ((*option_data)->bucket_filename);
01953   if ((*option_data)->coefts_filename != NULL) 
01954     free ((*option_data)->coefts_filename);
01955   if ((*option_data)->fitts_filename != NULL) 
01956     free ((*option_data)->fitts_filename);
01957   if ((*option_data)->sgnlts_filename != NULL) 
01958     free ((*option_data)->sgnlts_filename);
01959   if ((*option_data)->errts_filename != NULL) 
01960     free ((*option_data)->errts_filename);
01961   free (*option_data);  *option_data = NULL;
01962 
01963 
01964   /*----- Deallocate memory for fMRI time series data -----*/
01965   if (*fmri_data != NULL)    { free (*fmri_data);  *fmri_data = NULL; }
01966 
01967 
01968   /*----- Deallocate space for volume data -----*/
01969   if (*coef_vol  != NULL) 
01970     {
01971       for (ip = 0;  ip < p;  ip++)
01972         {
01973           if ((*coef_vol)[ip] != NULL)
01974             { free ((*coef_vol)[ip]);   (*coef_vol)[ip] = NULL; }
01975         }
01976       free (*coef_vol);   *coef_vol  = NULL;
01977     }
01978   
01979   if (*mse_vol   != NULL)    { free (*mse_vol);    *mse_vol   = NULL; }
01980   if (*ffull_vol != NULL)    { free (*ffull_vol);  *ffull_vol = NULL; }
01981   if (*rfull_vol != NULL)    { free (*rfull_vol);  *rfull_vol = NULL; } 
01982 
01983 
01984   /*----- Deallocate space for forward wavelet transform coefficients -----*/
01985   if (*coefts_vol != NULL)
01986     {
01987       for (it = 0;  it < N;  it++)
01988         { free ((*coefts_vol)[it]);   (*coefts_vol)[it] = NULL; }
01989       free (*coefts_vol);   *coefts_vol = NULL;
01990     }
01991 
01992 
01993   /*----- Deallocate space for wavelet filtered time series -----*/
01994   if (*fitts_vol != NULL)
01995     {
01996       for (it = 0;  it < N;  it++)
01997         { free ((*fitts_vol)[it]);   (*fitts_vol)[it] = NULL; }
01998       free (*fitts_vol);   *fitts_vol = NULL;
01999     }
02000 
02001 
02002   /*----- Deallocate space for signal model time series -----*/
02003   if (*sgnlts_vol != NULL)
02004     {
02005       for (it = 0;  it < N;  it++)
02006         { free ((*sgnlts_vol)[it]);   (*sgnlts_vol)[it] = NULL; }
02007       free (*sgnlts_vol);   *sgnlts_vol = NULL;
02008     }
02009 
02010 
02011   /*----- Deallocate space for residual errors time series -----*/
02012   if (*errts_vol != NULL)
02013     {
02014       for (it = 0;  it < N;  it++)
02015         { free ((*errts_vol)[it]);   (*errts_vol)[it] = NULL; }
02016       free (*errts_vol);   *errts_vol = NULL;
02017     }
02018 
02019 
02020 }
02021 
02022 
02023 /*---------------------------------------------------------------------------*/
02024 
02025 int main 
02026 (
02027   int argc,                /* number of input arguments */
02028   char ** argv             /* array of input arguments */ 
02029 )
02030 
02031 {
02032   WA_options * option_data;               /* wavelet analysis options */
02033   THD_3dim_dataset * dset_time = NULL;    /* input 3d+time data set */
02034   THD_3dim_dataset * mask_dset = NULL;    /* input mask data set */
02035   float * fmri_data = NULL;               /* input fMRI time series data */
02036   int fmri_length;                        /* length of fMRI time series */
02037 
02038   float ** coef_vol = NULL;   /* array of volumes for model parameters */
02039   float * mse_vol   = NULL;   /* volume of full model mean square error */
02040   float * ffull_vol = NULL;   /* volume of full model F-statistics */
02041   float * rfull_vol = NULL;   /* volume of full model R^2 stats. */
02042 
02043   float ** coefts_vol = NULL; /* volumes for wavelet transform coefs.*/
02044   float ** fitts_vol = NULL;  /* volumes for wavelet filtered time series */
02045   float ** sgnlts_vol = NULL;  /* volumes for signal model fit to input data */
02046   float ** errts_vol = NULL;   /* volumes for residual errors */
02047 
02048   
02049   /*----- Identify software -----*/
02050   printf ("\n\n");
02051   printf ("Program: %s \n", PROGRAM_NAME);
02052   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
02053   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
02054   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
02055   printf ("\n");
02056 
02057   
02058   /*----- Program initialization -----*/
02059   initialize_program (argc, argv, &option_data, 
02060                       &dset_time, &mask_dset, &fmri_data, &fmri_length,
02061                       &coef_vol, &mse_vol, &ffull_vol, &rfull_vol, 
02062                       &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol);
02063 
02064 
02065   /*----- Perform wavelet analysis -----*/
02066   calculate_results (option_data, dset_time, mask_dset, fmri_data, fmri_length,
02067                      coef_vol, mse_vol, ffull_vol, rfull_vol, 
02068                      coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
02069   
02070 
02071   /*----- Deallocate memory for input datasets -----*/   
02072   if (dset_time != NULL)  
02073     { THD_delete_3dim_dataset (dset_time, False);  dset_time = NULL; }
02074   if (mask_dset != NULL)  
02075     { THD_delete_3dim_dataset (mask_dset, False);  mask_dset = NULL; }
02076 
02077 
02078   /*----- Write requested output files -----*/
02079   if (option_data->input1D_filename == NULL)
02080     output_results (argc, argv, option_data, coef_vol, mse_vol, 
02081                     ffull_vol, rfull_vol, 
02082                     coefts_vol, fitts_vol, sgnlts_vol, errts_vol);
02083 
02084 
02085   /*----- Terminate program -----*/
02086   terminate_program (&option_data, &fmri_data, 
02087                      &coef_vol, &mse_vol, &ffull_vol, &rfull_vol,
02088                      &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol);
02089 
02090   exit(0);
02091 }
02092 
02093 
02094 
02095 
02096