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