00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008    This program generates an AFNI 3d+time data set.  The time series for
00009    each voxel is generated according to a user specified signal+noise model.
00010 
00011    File:     3dTSgen.c
00012    Author:   B. Douglas Ward
00013    Date:     17 June 1997
00014    
00015    Mod:      Function srand48 used for random number initialization.
00016    Date:     29 August 1997
00017 
00018    Mod:      Extensive changes required to implement the 'bucket' dataset.
00019    Date:     09 January 1998
00020 
00021    Mod:      Added the -inTR option.
00022              Removed duplicate readin of dset_time, and removed the
00023                input of dset_time's bricks (never used).
00024    Date:     22 July 1998 -- RWCox
00025              
00026    Mod:      Correction to initialization of dataset parameters.
00027    Date:     12 November 1998
00028 
00029    Mod:      Added changes for incorporating History notes.
00030    Date:     09 September 1999
00031 
00032 */
00033 
00034 /*---------------------------------------------------------------------------*/
00035 
00036 #define PROGRAM_NAME "3dTSgen"                       /* name of this program */
00037 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00038 #define PROGRAM_DATE "09 September 1999"         /* date of last program mod */
00039 
00040 
00041 #include <stdio.h>
00042 #include <math.h>
00043 #include <stdlib.h>
00044 
00045 #include "mrilib.h"
00046 
00047 #include "matrix.h"
00048 #include "simplex.h"
00049 #include "NLfit.h"
00050 
00051 #include "matrix.c"
00052 #include "simplex.c"
00053 #include "NLfit.c"
00054 
00055 
00056 typedef struct NL_options
00057 { 
00058   char * bucket_filename;      /* file name for bucket dataset */
00059   int numbricks;               /* number of sub-bricks in bucket dataset */
00060   int * brick_type;            /* indicates type of sub-brick */
00061   int * brick_coef;            /* regression coefficient number for sub-brick*/
00062   char ** brick_label;         /* character string label for sub-brick */
00063 
00064 } NL_options;
00065 
00066 /***** 22 July 1998 -- RWCox:
00067        Modified to allow DELT to be set from the TR of the input file *****/
00068 
00069 static float DELT = 1.0;   /* default */
00070 static int   inTR = 0 ;    /* set to 1 if -inTR option is used */
00071 static float dsTR = 0.0 ;  /* TR of the input file */
00072 
00073 
00074 static char * commandline = NULL ;         /* command line for history notes */
00075 
00076 
00077 /*---------------------------------------------------------------------------*/
00078 /*
00079   Routine to display 3dTSgen help menu.
00080 */
00081 
00082 void display_help_menu()
00083 {
00084   printf 
00085     (
00086      "This program generates an AFNI 3d+time data set.  The time series for \n"
00087      "each voxel is generated according to a user specified signal + noise  \n"
00088      "model.                                                              \n\n"
00089      "Usage:                                                                \n"
00090      "3dTSgen                                                               \n"
00091      "-input fname       fname = filename of prototype 3d + time data file  \n"
00092      "[-inTR]            set the TR of the created timeseries to be the TR  \n"
00093      "                     of the prototype dataset                         \n"
00094      "                     [The default is to compute with TR = 1.]         \n"
00095      "                     [The model functions are called for a  ]         \n"
00096      "                     [time grid of 0, TR, 2*TR, 3*TR, ....  ]         \n"
00097      "-signal slabel     slabel = name of (non-linear) signal model         \n"
00098      "-noise  nlabel     nlabel = name of (linear) noise model              \n"
00099      "-sconstr k c d     constraints for kth signal parameter:              \n"
00100      "                      c <= gs[k] <= d                                 \n"
00101      "-nconstr k c d     constraints for kth noise parameter:               \n"
00102      "                      c+b[k] <= gn[k] <= d+b[k]                       \n"
00103      "-sigma  s          s = std. dev. of additive Gaussian noise           \n"
00104      "[-voxel num]       screen output for voxel #num                       \n"
00105      "-output fname      fname = filename of output 3d + time data file     \n"
00106      "                                                                      \n"
00107      "                                                                      \n"
00108      "The following commands generate individual AFNI 1 sub-brick datasets: \n"
00109      "                                                                      \n"
00110      "[-scoef k fname]   write kth signal parameter gs[k];                  \n"
00111      "                     output 'fim' is written to prefix filename fname \n"
00112      "[-ncoef k fname]   write kth noise parameter gn[k];                   \n"
00113      "                     output 'fim' is written to prefix filename fname \n"
00114      "                                                                      \n"
00115      "                                                                      \n"
00116      "The following commands generate one AFNI 'bucket' type dataset:       \n"
00117      "                                                                      \n"
00118      "[-bucket n prefixname]   create one AFNI 'bucket' dataset containing  \n"
00119      "                           n sub-bricks; n=0 creates default output;  \n"
00120      "                           output 'bucket' is written to prefixname   \n"
00121      "The mth sub-brick will contain:                                       \n"
00122      "[-brick m scoef k label]   kth signal parameter regression coefficient\n"
00123      "[-brick m ncoef k label]   kth noise parameter regression coefficient \n"
00124     );
00125   
00126   exit(0);
00127 }
00128 
00129 
00130 /*---------------------------------------------------------------------------*/
00131      
00132 /** macro to test a malloc-ed pointer for validity **/
00133      
00134 #define MTEST(ptr) \
00135      if((ptr)==NULL) \
00136      ( fprintf(stderr,"*** Cannot allocate memory for statistics!\n"         \
00137                "*** Try using the -workmem option to reduce memory needs,\n" \
00138                "*** or create more swap space in the operating system.\n"    \
00139                ), exit(0) )
00140      
00141 
00142 /*---------------------------------------------------------------------------*/
00143 /*
00144   Routine to initialize the input options.
00145 */
00146  
00147 void initialize_options 
00148 (
00149   vfp * nmodel,            /* pointer to noise model */
00150   vfp * smodel,            /* pointer to signal model */  
00151   char *** npname,         /* noise parameter names */
00152   char *** spname,         /* signal parameter names */
00153   float ** min_nconstr,    /* minimum parameter constraints for noise model */
00154   float ** max_nconstr,    /* maximum parameter constraints for noise model */
00155   float ** min_sconstr,    /* minimum parameter constraints for signal model */
00156   float ** max_sconstr,    /* maximum parameter constraints for signal model */
00157   float * sigma,           /* std. dev. of additive Gaussian noise */
00158   int * nvoxel,            /* screen output for voxel #nvoxel */
00159   char ** input_filename,     /* file name of prototype 3d+time dataset */
00160   char ** output_filename,    /* file name for output 3d+time dataset */
00161   char *** ncoef_filename,    /* file name for noise model parameters */
00162   char *** scoef_filename,    /* file name for signal model parameters */
00163   NL_options * option_data    /* bucket dataset options */
00164 )
00165  
00166 {
00167   int ip;                     /* parameter index */
00168 
00169 
00170   /*----- initialize default values -----*/
00171   *sigma = 0.0;
00172   *nvoxel = -1;
00173   *smodel = NULL;
00174   *nmodel = NULL;
00175 
00176 
00177   /*----- allocate memory for parameter names -----*/
00178   *npname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00179   if (*npname == NULL)  
00180     NLfit_error ("Unable to allocate memory for noise parameter names");
00181   for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
00182     {
00183       (*npname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00184       if ((*npname)[ip] == NULL)  
00185         NLfit_error ("Unable to allocate memory for noise parameter names");
00186     }
00187 
00188   *spname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00189   if (*spname == NULL)  
00190     NLfit_error ("Unable to allocate memory for signal parameter names");
00191   for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
00192     {
00193       (*spname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00194       if ((*spname)[ip] == NULL)  
00195         NLfit_error ("Unable to allocate memory for signal parameter names");
00196     }
00197   
00198 
00199   /*----- allocate memory for parameter constraints -----*/
00200   *min_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00201   if (*min_nconstr == NULL)  
00202     NLfit_error ("Unable to allocate memory for min_nconstr");
00203   *max_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00204   if (*max_nconstr == NULL)
00205     NLfit_error ("Unable to allocate memory for max_nconstr");
00206   *min_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00207   if (*min_sconstr == NULL)  
00208     NLfit_error ("Unable to allocate memory for min_sconstr");
00209   *max_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS);
00210   if (*max_sconstr == NULL)
00211     NLfit_error ("Unable to allocate memory for max_sconstr");
00212 
00213 
00214   /*----- allocate memory space and initialize pointers for filenames -----*/
00215   *input_filename = NULL;
00216   *output_filename = NULL;  
00217   *ncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00218   if (*ncoef_filename == NULL)
00219     NLfit_error ("Unable to allocate memory for ncoef_filename");
00220   *scoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS);
00221   if (*scoef_filename == NULL)
00222     NLfit_error ("Unable to allocate memory for scoef_filename");
00223   for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
00224     {
00225       (*ncoef_filename)[ip] = NULL;
00226       (*scoef_filename)[ip] = NULL;
00227     }
00228 
00229 
00230   /*----- initialize bucket dataset options -----*/
00231   option_data->bucket_filename = NULL;
00232   option_data->numbricks = -1;
00233   option_data->brick_type = NULL;
00234   option_data->brick_coef = NULL;
00235   option_data->brick_label = NULL;
00236 
00237 }
00238 
00239 
00240 /*---------------------------------------------------------------------------*/
00241 /*
00242   Routine to get user specified input options.
00243 */
00244 
00245 void get_options
00246 (
00247   int argc,                /* number of input arguments */
00248   char ** argv,            /* array of input arguments */ 
00249   char ** nname,           /* name of noise model */
00250   char ** sname,           /* name of signal model */
00251   vfp * nmodel,            /* pointer to noise model */
00252   vfp * smodel,            /* pointer to signal model */  
00253   int * r,                 /* number of parameters in the noise model */
00254   int * p,                 /* number of parameters in the signal model */
00255   char *** npname,         /* noise parameter names */
00256   char *** spname,         /* signal parameter names */
00257   float ** min_nconstr,    /* minimum parameter constraints for noise model */
00258   float ** max_nconstr,    /* maximum parameter constraints for noise model */
00259   float ** min_sconstr,    /* minimum parameter constraints for signal model */
00260   float ** max_sconstr,    /* maximum parameter constraints for signal model */
00261   float * sigma,           /* std. dev. of additive Gaussian noise */
00262   int * nvoxel,            /* screen output for voxel #nvoxel */
00263   char ** input_filename,     /* file name of prototype 3d+time dataset */
00264   char ** output_filename,    /* file name for output 3d+time dataset */
00265   char *** ncoef_filename,    /* file name for noise model parameters */
00266   char *** scoef_filename,    /* file name for signal model parameters */
00267 
00268   int * nxyz,                       /* number of voxels in image */
00269   int * ts_length,                  /* length of time series data */  
00270   NL_options * option_data          /* bucket dataset options */
00271 )
00272 
00273 {
00274   const MAX_BRICKS = 100;           /* max. number of bricks in the bucket */
00275   int nopt = 1;                     /* input option argument counter */
00276   int ival, index;                  /* integer input */
00277   float fval;                       /* float input */
00278   char message[MAX_NAME_LENGTH];    /* error message */
00279   int ok;                           /* boolean for specified model exists */
00280   THD_3dim_dataset * dset_time;     /* prototype 3d+time data set */
00281 
00282   NLFIT_MODEL_array * model_array = NULL;   /* array of SO models */
00283   int im;                                   /* model index */
00284   int ibrick;                       /* sub-brick index */
00285   int nbricks;                      /* number of bricks in the bucket */
00286 
00287   
00288   /*----- does user request help menu? -----*/
00289   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00290   
00291 
00292   /*----- initialize the model array -----*/
00293   model_array = NLFIT_get_many_MODELs ();
00294   if ((model_array == NULL) || (model_array->num == 0))
00295     NLfit_error ("Unable to locate any models");
00296 
00297 
00298   /*----- initialize the input options -----*/
00299   initialize_options (nmodel, smodel, npname, spname,
00300                       min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00301                       sigma, nvoxel, input_filename, 
00302                       output_filename, ncoef_filename, scoef_filename,
00303                       option_data); 
00304   
00305   /*----- main loop over input options -----*/
00306   while (nopt < argc )
00307     {
00308 
00309       /*-----   -input filename   -----*/
00310       if (strncmp(argv[nopt], "-input", 6) == 0)
00311         {
00312           nopt++;
00313           if (nopt >= argc)  NLfit_error ("need argument after -input ");
00314           *input_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00315           if (*input_filename == NULL)
00316             NLfit_error ("Unable to allocate memory for input_filename");
00317           strcpy (*input_filename, argv[nopt]);
00318 
00319           /*----- initialize data set parameters -----*/
00320           dset_time = THD_open_one_dataset (*input_filename);
00321           if (dset_time == NULL)  
00322             { 
00323               sprintf (message, 
00324                        "Unable to open data file: %s", *input_filename);
00325               NLfit_error (message);
00326             }
00327           *nxyz =  dset_time->dblk->diskptr->dimsizes[0]
00328             * dset_time->dblk->diskptr->dimsizes[1]
00329             * dset_time->dblk->diskptr->dimsizes[2] ;
00330           *ts_length = DSET_NUM_TIMES(dset_time);
00331 
00332           dsTR = DSET_TIMESTEP(dset_time) ;
00333           if( DSET_TIMEUNITS(dset_time) == UNITS_MSEC_TYPE ) dsTR *= 0.001 ;
00334 
00335           THD_delete_3dim_dataset(dset_time, False);  dset_time = NULL ;
00336 
00337           nopt++;
00338           continue;
00339         }
00340 
00341       /*----- 22 July 1998: the -inTR option -----*/
00342 
00343       if( strncmp(argv[nopt],"-inTR",5) == 0 ){
00344          inTR = 1 ;
00345          nopt++ ; continue ;
00346       }
00347       
00348       /*-----   -signal slabel  -----*/
00349       if (strncmp(argv[nopt], "-signal", 7) == 0)
00350         {
00351           nopt++;
00352           if (nopt >= argc)  NLfit_error ("need argument after -signal ");
00353           *sname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00354           if (*sname == NULL)
00355             NLfit_error ("Unable to allocate memory for signal model name");
00356           strcpy (*sname, argv[nopt]);
00357           initialize_signal_model (model_array, *sname, 
00358                                    smodel, p, *spname,
00359                                    *min_sconstr, *max_sconstr);
00360           nopt++;
00361           continue;
00362         }
00363       
00364       
00365       /*-----   -noise nlabel  -----*/
00366       if (strncmp(argv[nopt], "-noise", 6) == 0)
00367         {
00368           nopt++;
00369           if (nopt >= argc)  NLfit_error ("need argument after -noise ");
00370           *nname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00371           if (*nname == NULL)
00372             NLfit_error ("Unable to allocate memory for noise model name");
00373           strcpy (*nname, argv[nopt]);
00374           initialize_noise_model (model_array, *nname, 
00375                                   nmodel, r, *npname,
00376                                   *min_nconstr, *max_nconstr);
00377           nopt++;
00378           continue;
00379         }
00380       
00381       
00382       /*-----   -sconstr k min max  -----*/
00383       if (strncmp(argv[nopt], "-sconstr", 8) == 0 || strncmp(argv[nopt],"-scnstr",8) == 0 )
00384         {
00385           nopt++;
00386           if (nopt+2 >= argc)  NLfit_error("need 3 arguments after -sconstr ");
00387 
00388           sscanf (argv[nopt], "%d", &ival);
00389           if ((ival < 0) || (ival >= *p))
00390             NLfit_error ("illegal argument after -sconstr ");
00391           index = ival;
00392           nopt++;
00393 
00394           sscanf (argv[nopt], "%f", &fval); 
00395           (*min_sconstr)[index] = fval;
00396           nopt++;
00397 
00398           sscanf (argv[nopt], "%f", &fval); 
00399           (*max_sconstr)[index] = fval;
00400           nopt++;
00401           continue;
00402         }
00403       
00404       
00405       /*-----   -nconstr k min max  -----*/
00406       if (strncmp(argv[nopt], "-nconstr", 8) == 0 || strncmp(argv[nopt],"-ncnstr",8) == 0 )
00407         {
00408           nopt++;
00409           if (nopt+2 >= argc)  NLfit_error("need 3 arguments after -nconstr ");
00410 
00411           sscanf (argv[nopt], "%d", &ival);
00412           if ((ival < 0) || (ival >= *r))
00413             NLfit_error ("illegal argument after -nconstr ");
00414           index = ival;
00415           nopt++;
00416 
00417           sscanf (argv[nopt], "%f", &fval); 
00418           (*min_nconstr)[index] = fval;
00419           nopt++;
00420 
00421           sscanf (argv[nopt], "%f", &fval); 
00422           (*max_nconstr)[index] = fval;
00423           nopt++;
00424           continue;
00425         }
00426       
00427       
00428        /*-----   -sigma s  -----*/
00429       if (strncmp(argv[nopt], "-sigma", 6) == 0)
00430         {
00431           nopt++;
00432           if (nopt >= argc)  NLfit_error ("need argument after -sigma ");
00433           sscanf (argv[nopt], "%f", &fval); 
00434           if (fval < 0.0)
00435             NLfit_error ("illegal argument after -sigma ");
00436           *sigma = fval;
00437           nopt++;
00438           continue;
00439         }
00440       
00441 
00442       /*-----   -voxel num  -----*/
00443       if (strncmp(argv[nopt], "-voxel", 6) == 0)
00444         {
00445           nopt++;
00446           if (nopt >= argc)  NLfit_error ("need argument after -voxel ");
00447           sscanf (argv[nopt], "%d", &ival);
00448           if (ival <= 0)
00449             NLfit_error ("illegal argument after -voxel ");
00450           *nvoxel = ival;
00451           nopt++;
00452           continue;
00453         }
00454       
00455       
00456        /*-----   -output filename   -----*/
00457       if (strncmp(argv[nopt], "-output", 7) == 0)
00458         {
00459           nopt++;
00460           if (nopt >= argc)  NLfit_error ("need argument after -output ");
00461           *output_filename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00462           if (*output_filename == NULL)
00463             NLfit_error ("Unable to allocate memory for output_filename");
00464           strcpy (*output_filename, argv[nopt]);
00465           nopt++;
00466           continue;
00467         }
00468       
00469 
00470      /*-----   -scoef k filename   -----*/
00471       if (strncmp(argv[nopt], "-scoef", 6) == 0)
00472         {
00473           nopt++;
00474           if (nopt+1 >= argc)  NLfit_error ("need 2 arguments after -scoef ");
00475           sscanf (argv[nopt], "%d", &ival);
00476           if ((ival < 0) || (ival >= *p))
00477             NLfit_error ("illegal argument after -scoef ");
00478           index = ival;
00479           nopt++;
00480 
00481           (*scoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00482           if ((*scoef_filename)[index] == NULL)
00483             NLfit_error ("Unable to allocate memory for scoef_filename");
00484           strcpy ((*scoef_filename)[index], argv[nopt]);
00485 
00486           nopt++;
00487           continue;
00488         }
00489       
00490 
00491      /*-----   -ncoef k filename   -----*/
00492       if (strncmp(argv[nopt], "-ncoef", 7) == 0)
00493         {
00494           nopt++;
00495           if (nopt+1 >= argc)  NLfit_error ("need 2 arguments after -ncoef ");
00496           sscanf (argv[nopt], "%d", &ival);
00497           if ((ival < 0) || (ival >= *r))
00498             NLfit_error ("illegal argument after -ncoef ");
00499           index = ival;
00500           nopt++;
00501 
00502           (*ncoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH);
00503           if ((*ncoef_filename)[index] == NULL)
00504             NLfit_error ("Unable to allocate memory for ncoef_filename");
00505           strcpy ((*ncoef_filename)[index], argv[nopt]);
00506 
00507           nopt++;
00508           continue;
00509         }
00510       
00511 
00512       /*----- -bucket n prefixname -----*/
00513       if (strncmp(argv[nopt], "-bucket", 7) == 0)
00514         {
00515           nopt++;
00516           if (nopt+1 >= argc)  NLfit_error ("need 2 arguments after -bucket ");
00517           sscanf (argv[nopt], "%d", &ival);
00518           if ((ival < 0) || (ival > MAX_BRICKS))
00519             NLfit_error ("illegal argument after -bucket ");
00520           nopt++;
00521 
00522           option_data->bucket_filename = 
00523             malloc (sizeof(char) * MAX_NAME_LENGTH);
00524           if (option_data->bucket_filename == NULL)
00525             NLfit_error ("Unable to allocate memory for bucket_filename");
00526           strcpy (option_data->bucket_filename, argv[nopt]);
00527           
00528           /*----- set number of sub-bricks in the bucket -----*/
00529           if (ival == 0)
00530             nbricks = (*p) + (*r);
00531           else
00532             nbricks = ival;
00533           option_data->numbricks = nbricks;
00534           
00535           /*----- allocate memory and initialize bucket dataset options -----*/
00536           option_data->brick_type = malloc (sizeof(int) * nbricks);
00537           option_data->brick_coef = malloc (sizeof(int) * nbricks);
00538           option_data->brick_label = malloc (sizeof(char *) * nbricks);
00539           for (ibrick = 0;  ibrick < nbricks;  ibrick++)
00540             {
00541               option_data->brick_type[ibrick] = -1;
00542               option_data->brick_coef[ibrick] = -1;
00543               option_data->brick_label[ibrick] = 
00544                 malloc (sizeof(char) * MAX_NAME_LENGTH);
00545             }
00546           
00547 
00548           if (ival == 0)   
00549             /*----- throw  (almost) everything into the bucket -----*/
00550             {
00551               for (ibrick = 0;  ibrick < (*r);  ibrick++)
00552                 {
00553                   option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00554                   option_data->brick_coef[ibrick] = ibrick;
00555                   strcpy (option_data->brick_label[ibrick], (*npname)[ibrick]);
00556                 }
00557               
00558               for (ibrick = (*r);  ibrick < (*p) + (*r);  ibrick++)
00559                 {
00560                   option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00561                   option_data->brick_coef[ibrick] = ibrick;
00562                   strcpy (option_data->brick_label[ibrick],
00563                           (*spname)[ibrick-(*r)]);
00564                 }             
00565             }
00566 
00567           nopt++;
00568           continue;
00569         }
00570 
00571 
00572       /*----- -brick m type k label -----*/
00573       if (strncmp(argv[nopt], "-brick", 6) == 0)
00574         {
00575           nopt++;
00576           if (nopt+2 >= argc)  
00577             NLfit_error ("need more arguments after -brick ");
00578           sscanf (argv[nopt], "%d", &ibrick);
00579           if ((ibrick < 0) || (ibrick >= option_data->numbricks))
00580             NLfit_error ("illegal argument after -brick ");
00581           nopt++;
00582 
00583           if (strncmp(argv[nopt], "scoef", 4) == 0)
00584             {
00585               option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00586 
00587               nopt++;
00588               sscanf (argv[nopt], "%d", &ival);
00589               if ((ival < 0) || (ival > (*p)))
00590                 NLfit_error ("illegal argument after scoef ");
00591               option_data->brick_coef[ibrick] = ival + (*r);
00592               
00593               nopt++;
00594               if (nopt >= argc)  
00595                 NLfit_error ("need more arguments after -brick ");
00596               strcpy (option_data->brick_label[ibrick], argv[nopt]); 
00597             }
00598 
00599           else if (strncmp(argv[nopt], "ncoef", 4) == 0)
00600             {
00601               option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
00602 
00603               nopt++;
00604               sscanf (argv[nopt], "%d", &ival);
00605               if ((ival < 0) || (ival > (*r)))
00606                 NLfit_error ("illegal argument after ncoef ");
00607               option_data->brick_coef[ibrick] = ival;
00608               
00609               nopt++;
00610               if (nopt >= argc)  
00611                 NLfit_error ("need more arguments after -brick ");
00612               strcpy (option_data->brick_label[ibrick], argv[nopt]); 
00613             }
00614 
00615           else  NLfit_error ("unable to interpret options after -brick ");
00616                   
00617           printf ("ibrick = %d \n", ibrick);
00618           printf ("brick_type  = %d \n", option_data->brick_type[ibrick]);
00619           printf ("brick_coef  = %d \n", option_data->brick_coef[ibrick]);
00620           printf ("brick_label = %s \n", option_data->brick_label[ibrick]);
00621           
00622           nopt++;
00623           continue;
00624         }
00625      
00626       
00627       /*----- unknown command -----*/
00628 
00629       { char buf[256] ;
00630         sprintf(buf,"unrecognized command line option: %s",argv[nopt]) ;
00631         NLfit_error (buf);
00632       }
00633     }
00634 
00635   
00636   /*----- discard the model array -----*/
00637   DESTROY_MODEL_ARRAY (model_array) ;
00638   
00639 }
00640 
00641 
00642 /*---------------------------------------------------------------------------*/
00643 /*
00644   Routine to check whether one output file already exists.
00645 */
00646 
00647 void check_one_output_file 
00648 (
00649   THD_3dim_dataset * dset,          /* prototype 3d+time data set */
00650   char * filename                   /* name of output file */
00651 )
00652 
00653 {
00654   THD_3dim_dataset * new_dset = NULL;   /* output afni data set pointer */
00655   int ierror;                           /* number of errors in editing data */
00656   
00657   
00658   /*----- make an empty copy of input dataset -----*/
00659   new_dset = EDIT_empty_copy (dset);  
00660   
00661   ierror = EDIT_dset_items( new_dset,
00662                             ADN_prefix, filename ,
00663                             ADN_label1, filename ,
00664                             ADN_self_name, filename ,
00665                             ADN_type, ISHEAD(dset) ? HEAD_FUNC_TYPE : 
00666                                                      GEN_FUNC_TYPE ,
00667                             ADN_none ) ;
00668   
00669   if( ierror > 0 )
00670     {
00671       fprintf(stderr,
00672               "*** %d errors in attempting to create output dataset!\n", 
00673               ierror ) ;
00674       exit(1) ;
00675     }
00676   
00677   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00678     {
00679       fprintf(stderr,
00680               "*** Output dataset file %s already exists"
00681               "--cannot continue!\a\n",
00682               new_dset->dblk->diskptr->header_name ) ;
00683       exit(1) ;
00684     }
00685   
00686   /*----- deallocate memory -----*/   
00687   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00688   
00689 }
00690 
00691 
00692 /*---------------------------------------------------------------------------*/
00693 /*
00694   Routine to check whether output files already exist.
00695 */
00696 
00697 void check_output_files 
00698 (
00699   char * input_filename,         /* file name for input 3d+time data set */
00700   char * output_filename,        /* file name for output 3d+time data set */
00701   char ** ncoef_filename,        /* file name for noise model parameters */
00702   char ** scoef_filename,        /* file name for signal model parameters */
00703   char * bucket_filename         /* file name for bucket dataset */
00704 )
00705 
00706 {
00707   THD_3dim_dataset * dset_time;  /* prototype 3d+time data set */
00708   int ip;                        /* parameter index */
00709   
00710 
00711   dset_time = THD_open_one_dataset (input_filename);
00712 
00713 
00714   if (output_filename != NULL)   
00715     check_one_output_file (dset_time, output_filename);
00716   
00717   for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
00718     {
00719       if (ncoef_filename[ip] != NULL)
00720         check_one_output_file (dset_time, ncoef_filename[ip]);
00721       if (scoef_filename[ip] != NULL)
00722         check_one_output_file (dset_time, scoef_filename[ip]);
00723     }
00724 
00725 
00726   if (bucket_filename != NULL)   
00727     check_one_output_file (dset_time, bucket_filename);
00728 
00729   THD_delete_3dim_dataset (dset_time, False);  dset_time = NULL ;
00730 }
00731 
00732 
00733 /*---------------------------------------------------------------------------*/
00734 /*
00735   Routine to check for valid inputs.
00736 */
00737   
00738 void check_for_valid_inputs 
00739 (
00740   int r,                  /* number of parameters in the noise model */
00741   int p,                  /* number of parameters in the signal model */
00742   float * min_nconstr,    /* minimum parameter constraints for noise model */
00743   float * max_nconstr,    /* maximum parameter constraints for noise model */
00744   float * min_sconstr,    /* minimum parameter constraints for signal model */
00745   float * max_sconstr,    /* maximum parameter constraints for signal model */
00746 
00747   char * input_filename,          /* file name of prototype 3d+time dataset */
00748   char * output_filename,         /* file name for output 3d+time data set */
00749   char ** ncoef_filename,         /* file name for noise model parameters */
00750   char ** scoef_filename,         /* file name for signal model parameters */
00751   char * bucket_filename          /* file name for bucket dataset */
00752 )
00753 
00754 {
00755   int ip;                         /* parameter index */
00756 
00757 
00758   /*----- check for valid constraints -----*/
00759   for (ip = 0;  ip < r;  ip++)
00760     if (min_nconstr[ip] > max_nconstr[ip])
00761       NLfit_error ("Must have minimum constraints <= maximum constraints");
00762   for (ip = 0;  ip < p;  ip++)
00763     if (min_sconstr[ip] > max_sconstr[ip])
00764       NLfit_error("Must have mininum constraints <= maximum constraints");
00765       
00766 
00767   /*----- check whether any of the output files already exist -----*/
00768   check_output_files (input_filename, output_filename, 
00769                       ncoef_filename, scoef_filename, bucket_filename);
00770 
00771 }
00772 
00773 
00774 /*---------------------------------------------------------------------------*/
00775 /*
00776   Perform all program initialization.
00777 */
00778 
00779 void initialize_program 
00780 (
00781   int argc,                /* number of input arguments */
00782   char ** argv,            /* array of input arguments */ 
00783   char ** nname,           /* name of noise model */
00784   char ** sname,           /* name of signal model */
00785   vfp * nmodel,            /* pointer to noise model */
00786   vfp * smodel,            /* pointer to signal model */  
00787   int *r,                  /* number of parameters in the noise model */
00788   int *p,                  /* number of parameters in the signal model */
00789   char *** npname,         /* noise parameter names */
00790   char *** spname,         /* signal parameter names */
00791   float ** min_nconstr,    /* minimum parameter constraints for noise model */
00792   float ** max_nconstr,    /* maximum parameter constraints for noise model */
00793   float ** min_sconstr,    /* minimum parameter constraints for signal model */
00794   float ** max_sconstr,    /* maximum parameter constraints for signal model */
00795   float * sigma,           /* std. dev. of additive Gaussian noise */
00796   int * nvoxel,            /* screen output for voxel #nvoxel */
00797 
00798   char ** input_filename,     /* file name of prototype 3d+time dataset */
00799   char ** output_filename,    /* file name for output 3d+time dataset */
00800   char *** ncoef_filename,    /* file name for noise model parameters */
00801   char *** scoef_filename,    /* file name for signal model parameters */
00802 
00803   int * nxyz,                       /* number of voxels in image */
00804   int * ts_length,                  /* length of time series data */  
00805 
00806   float *** x_array,       /* independent variable matrix */
00807   float ** ts_array,       /* input time series array */
00808   short *** d_array,       /* output times series volume */
00809 
00810   float ** par_full,       /* estimated parameters for the full model */
00811 
00812   float *** ncoef_vol,     /* volume of noise model parameters */
00813   float *** scoef_vol,     /* volume of signal model parameters */
00814 
00815   NL_options * option_data          /* bucket dataset options */
00816 
00817 )
00818      
00819 {
00820   int dimension;           /* dimension of full model */
00821   int ip;                  /* parameter index */
00822   int it;                  /* time index */
00823 
00824 
00825   /*----- save command line for history notes -----*/
00826   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00827 
00828   
00829   /*----- get command line inputs -----*/
00830   get_options (argc, argv, nname, sname, nmodel, smodel, 
00831                r, p, npname, spname,
00832                min_nconstr, max_nconstr, min_sconstr, max_sconstr, 
00833                sigma, nvoxel,
00834                input_filename, output_filename, ncoef_filename, scoef_filename,
00835                nxyz, ts_length, option_data);
00836 
00837   /*----- check for valid inputs -----*/
00838   check_for_valid_inputs (*r, *p, *min_nconstr, *max_nconstr, 
00839                           *min_sconstr, *max_sconstr, *input_filename, 
00840                           *output_filename, *ncoef_filename, *scoef_filename,
00841                           option_data->bucket_filename);
00842 
00843   /*----- allocate space for input time series -----*/
00844   *ts_array = (float *) malloc (sizeof(float) * (*ts_length));
00845   if (*ts_array == NULL)
00846     NLfit_error ("Unable to allocate memory for ts_array");
00847 
00848   
00849   /*----- allocate space for independent variable matrix -----*/
00850   *x_array = (float **) malloc (sizeof(float *) * (*ts_length));
00851   if (*x_array == NULL)
00852     NLfit_error ("Unable to allocate memory for x_array");
00853   for (it = 0;  it < (*ts_length);  it++)
00854     {
00855       (*x_array)[it] = (float *) malloc (sizeof(float) * 3);
00856       if ((*x_array)[it] == NULL)
00857         NLfit_error ("Unable to allocate memory for x_array[it]");
00858     }
00859     
00860   /*----- initialize independent variable matrix -----*/
00861 
00862   if( inTR && dsTR > 0.0 ){   /* 22 July 1998 */
00863      DELT = dsTR ;
00864      fprintf(stderr,"--- computing with TR = %g\n",DELT) ;
00865   }
00866 
00867   for (it = 0;  it < (*ts_length);  it++)  
00868     {
00869       (*x_array)[it][0] = 1.0;
00870       (*x_array)[it][1] = it * DELT;
00871       (*x_array)[it][2] = (it * DELT) * (it * DELT);
00872     }
00873   
00874   /*----- allocate memory space for parameters -----*/
00875   dimension = (*r) + (*p);
00876   *par_full = (float *) malloc (sizeof(float) * dimension);
00877   if (*par_full == NULL)
00878     NLfit_error ("Unable to allocate memory for par_full");
00879 
00880   *ncoef_vol = (float **) malloc (sizeof(float *) * (*r));
00881   if (*ncoef_vol == NULL)
00882     NLfit_error ("Unable to allocate memory for ncoef_vol");
00883   for (ip = 0;  ip < (*r);  ip++)
00884     {
00885       (*ncoef_vol)[ip] = (float *) malloc (sizeof(float) * (*nxyz));
00886       if ((*ncoef_vol)[ip] == NULL)
00887         NLfit_error ("Unable to allocate memory for ncoef_vol[ip]");
00888     }
00889   
00890   *scoef_vol = (float **) malloc (sizeof(float *) * (*p));
00891   if (*scoef_vol == NULL)
00892     NLfit_error ("Unable to allocate memory for scoef_vol");
00893   for (ip = 0;  ip < (*p);  ip++)
00894     {
00895       (*scoef_vol)[ip] = (float *) malloc (sizeof(float) * (*nxyz));
00896       if ((*scoef_vol)[ip] == NULL)
00897         NLfit_error ("Unable to allocate memory for scoef_vol[ip]");
00898     }
00899 
00900   
00901   /*----- allocate space for output time series volume -----*/
00902   *d_array = (short **) malloc (sizeof(short *) * (*ts_length));
00903   if (*d_array == NULL)
00904     NLfit_error ("Unable to allocate memory for d_array");
00905   for (it = 0;  it < *ts_length;  it++)
00906     {
00907       (*d_array)[it] = (short *) malloc (sizeof(short) * (*nxyz));
00908       if ((*d_array)[it] == NULL)
00909         NLfit_error ("Unable to allocate memory for d_array[it]");
00910     }
00911 
00912 
00913   /*----- initialize random number generator -----*/
00914   srand48 (1234567);
00915  
00916 }
00917 
00918 
00919 /*---------------------------------------------------------------------------*/
00920 /*
00921   Generate artificial time series array.
00922 */
00923 
00924 void generate_ts_array 
00925 (
00926   vfp nmodel,             /* pointer to noise model */
00927   vfp smodel,             /* pointer to signal model */
00928   int r,                  /* number of parameters in the noise model */
00929   int p,                  /* number of parameters in the signal model */
00930   float * min_nconstr,    /* minimum parameter constraints for noise model */
00931   float * max_nconstr,    /* maximum parameter constraints for noise model */
00932   float * min_sconstr,    /* minimum parameter constraints for signal model */
00933   float * max_sconstr,    /* maximum parameter constraints for signal model */
00934   int ts_length,          /* length of time series array */
00935   float ** x_array,       /* independent variable matrix */
00936  
00937   float sigma,            /* std. dev. of additive Gaussian noise */
00938 
00939   float * par_true,       /* true parameters for this time series */
00940   float * ts_array        /* artificially generated time series data */
00941 )
00942 
00943 {
00944   int ip;                 /* parameter index */
00945   int it;                 /* time index */
00946   float n1, n2;           /* normal random variates */
00947 
00948 
00949   /*----- select 'true' noise parameters -----*/
00950   for (ip = 0;  ip < r;  ip++)
00951     par_true[ip] = get_random_value (min_nconstr[ip], max_nconstr[ip]);
00952 
00953   /*----- select 'true' signal parameters -----*/
00954   for (ip = 0;  ip < p;  ip++)
00955     par_true[ip+r] = get_random_value (min_sconstr[ip], max_sconstr[ip]);
00956   
00957   /*----- calculate the 'true' time series using the 'true' parameters -----*/
00958   full_model (nmodel, smodel, par_true, par_true + r, 
00959               ts_length, x_array, ts_array);
00960 
00961   /*----- add Gaussian noise -----*/
00962   for (it = 0;  it < ts_length;  it++)
00963     {
00964       normal (&n1, &n2);
00965       ts_array[it] += n1*sigma;
00966     }
00967 
00968 }
00969 
00970 
00971 /*---------------------------------------------------------------------------*/
00972 /*
00973   Save time series into AFNI 3d+time data set.
00974 */
00975 
00976 void save_ts_array 
00977 (
00978   int iv,                  /* save time series for this voxel */
00979   int ts_length,           /* length of time series array */
00980   float * ts_array,        /* time series data for voxel #iv */
00981   short ** d_array         /* output time series volume */
00982 )
00983 
00984 {
00985   int it;
00986 
00987 
00988   for (it = 0;  it < ts_length;  it++)
00989     {
00990       d_array[it][iv] = (short) ts_array[it];
00991     }
00992         
00993 }
00994 
00995 
00996 /*---------------------------------------------------------------------------*/
00997 /*
00998   Print out the given time series.
00999 */
01000 
01001 void write_ts_array 
01002 (
01003   int iv,                  /* write time series for this voxel */
01004   int ts_length,           /* length of time series array */
01005   float * ts_array,        /* time series data for voxel #iv */
01006   short ** d_array         /* output time series volume */
01007 )
01008 
01009 {
01010   int it;                      /* time index */
01011 
01012   for (it = 0;  it < ts_length;  it++)
01013     printf ("ts_array[%d] = %9.3f   d_array[%d][%d] = %d    \n", 
01014             it, ts_array[it], it, iv+1, d_array[it][iv]);
01015       
01016 }
01017 
01018 
01019 /*---------------------------------------------------------------------------*/
01020 /*
01021   Write parameters for this voxel.
01022 */
01023 
01024 void write_parameters 
01025 (
01026   int iv,                  /* current voxel number */
01027   char * nname,            /* name of noise model */
01028   char * sname,            /* name of signal model */
01029   int r,                   /* number of parameters in the noise model */
01030   int p,                   /* number of parameters in the signal model */
01031   char ** npname,          /* noise parameter names */
01032   char ** spname,          /* signal parameter names */
01033   float * par_full         /* estimated parameters for the full model */
01034 )
01035 
01036 {
01037   int ip;                  /* parameter index */
01038 
01039 
01040   if (iv < 0)
01041     printf ("\n\nVoxel Results: \n");
01042   else
01043     printf ("\n\nVoxel #%d\n", iv+1);
01044       
01045 
01046   /*----- write full model parameters -----*/  
01047   printf ("\nFull (%s + %s) Model: \n", nname, sname);
01048   
01049   for (ip = 0;  ip < r;  ip++)
01050     printf ("gn[%d]=%12.6f  %s \n", ip, par_full[ip], npname[ip]);
01051     
01052   for (ip = 0;  ip < p;  ip++)
01053     printf ("gs[%d]=%12.6f  %s \n", ip, par_full[ip+r], spname[ip]);                
01054 }
01055 
01056 
01057 /*---------------------------------------------------------------------------*/
01058 /*
01059   Save parameters for output later.
01060 */
01061 
01062 void save_parameters 
01063 (
01064   int iv,                  /* current voxel number */
01065   int r,                   /* number of parameters in the noise model */
01066   int p,                   /* number of parameters in the signal model */
01067   float * par_full,        /* estimated parameters for the full model */
01068 
01069   float ** ncoef_vol,      /* volume of noise model parameters */
01070   float ** scoef_vol       /* volume of signal model parameters */
01071 )
01072 
01073 {
01074   int ip;                  /* parameter index */
01075 
01076 
01077   /*----- save noise parameter estimates -----*/
01078   for (ip = 0;  ip < r;  ip++)
01079     {
01080       ncoef_vol[ip][iv] = par_full[ip];
01081     }
01082       
01083   /*----- save signal parameter estimates -----*/
01084   for (ip = 0;  ip < p;  ip++)
01085     {
01086       scoef_vol[ip][iv] = par_full[ip+r];
01087     }
01088 
01089 }
01090 
01091 
01092 /*---------------------------------------------------------------------------*/
01093 /*
01094   Convert one volume to another type, autoscaling:
01095      nxy   = # voxels
01096      itype = input datum type
01097      ivol  = pointer to input volume
01098      otype = output datum type
01099      ovol  = pointer to output volume (again, must be pre-malloc-ed)
01100   Return value is the scaling factor used (0.0 --> no scaling).
01101 */
01102 
01103 
01104 float EDIT_coerce_autoscale_new( int nxyz ,
01105                                  int itype,void *ivol , int otype,void *ovol )
01106 {
01107   float fac=0.0 , top ;
01108   
01109   if( MRI_IS_INT_TYPE(otype) ){
01110     top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
01111     if (top == 0.0)  fac = 0.0;
01112     else  fac = MRI_TYPE_maxval[otype]/top;
01113   }
01114   
01115   EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
01116   return ( fac );
01117 }
01118 
01119 
01120 /*---------------------------------------------------------------------------*/
01121 /*
01122   Routine to write one AFNI 3d+time data set. 
01123 */
01124 
01125 
01126 void output_ts_array 
01127 (
01128   short ** d_array,                   /* output time series volume data */
01129   int  ts_length,                     /* length of time series data */  
01130   char * input_filename,              /* input afni data set file name */
01131   char * filename                     /* output afni data set file name */
01132 )
01133 
01134 {
01135   THD_3dim_dataset * dset = NULL;     /* input afni data set pointer */
01136   THD_3dim_dataset * new_dset = NULL; /* output afni data set pointer */
01137   int ib;                             /* sub-brick index */ 
01138   int ierror;                         /* number of errors in editing data */
01139   float * fbuf;                       /* float buffer */
01140   float fimfac;                       /* scale factor for short data */
01141   char label[80];                     /* label for output file history */ 
01142   
01143      
01144   /*----- allocate memory -----*/
01145   fbuf = (float *)  malloc (sizeof(float) * ts_length);   MTEST (fbuf);
01146   for (ib = 0;  ib < ts_length;  ib++)    fbuf[ib] = 0.0;
01147 
01148   
01149   /*-- make an empty copy of the prototype dataset, for eventual output --*/
01150   dset = THD_open_one_dataset (input_filename);
01151   new_dset = EDIT_empty_copy (dset);
01152 
01153 
01154   /*----- Record history of dataset -----*/
01155 
01156   sprintf (label, "Output prefix: %s", filename);
01157   if( commandline != NULL )
01158      tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
01159   else
01160      tross_Append_History ( new_dset, label);
01161 
01162 
01163   /*----- delete prototype dataset -----*/
01164   THD_delete_3dim_dataset (dset, False);  dset = NULL ;
01165   
01166 
01167   ierror = EDIT_dset_items( new_dset ,
01168                             ADN_prefix , filename ,
01169                             ADN_label1 , filename ,
01170                             ADN_self_name , filename ,
01171                                     ADN_malloc_type, DATABLOCK_MEM_MALLOC ,  
01172                             ADN_nvals , ts_length ,
01173                             ADN_none ) ;
01174   
01175   if( ierror > 0 ){
01176     fprintf(stderr,
01177           "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01178     exit(1) ;
01179   }
01180   
01181   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01182     fprintf(stderr,
01183             "*** Output dataset file %s already exists--cannot continue!\a\n",
01184             new_dset->dblk->diskptr->header_name ) ;
01185     exit(1) ;
01186   }
01187 
01188   
01189   /*----- attach bricks to new data set -----*/
01190   for (ib = 0;  ib < ts_length;  ib++)
01191     mri_fix_data_pointer (d_array[ib], DSET_BRICK(new_dset,ib)); 
01192   
01193   
01194   /*----- write afni data set -----*/
01195 
01196   (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01197   
01198   THD_load_statistics( new_dset ) ;
01199   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01200   fprintf(stderr,"--- Wrote combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
01201 
01202   
01203   /*----- deallocate memory -----*/   
01204   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01205   free (fbuf);   fbuf = NULL;
01206 
01207 }
01208 
01209 
01210 /*---------------------------------------------------------------------------*/
01211 /*
01212   Routine to write one AFNI data set.  This data set must be of type 'fim'.
01213   The intensity data is in ffim.
01214 */
01215 
01216 
01217 void write_afni_data (char * input_filename, int nxyz, char * filename,  
01218                       float * ffim)
01219 {
01220   int ii;                             /* voxel index */
01221   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
01222   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
01223   int iv;                             /* sub-brick index */ 
01224   int ierror;                         /* number of errors in editing data */
01225   int ibuf[32];                       /* integer buffer */
01226   float fbuf[MAX_STAT_AUX];           /* float buffer */
01227   float fimfac;                       /* scale factor for short data */
01228   int output_datum;                   /* data type for output data */
01229   void  * vdif;                       /* 1st sub-brick data pointer */
01230   int func_type;                      /* afni data set type */
01231   float top, func_scale_short;        /* parameters for scaling data */
01232   char label[80];                     /* label for output file history */ 
01233   
01234     
01235   /*----- read input dataset -----*/
01236   dset = THD_open_one_dataset (input_filename) ;
01237   if( ! ISVALID_3DIM_DATASET(dset) ){
01238     fprintf(stderr,"*** Unable to open dataset file %s\n",
01239             input_filename);
01240     exit(1) ;
01241   }
01242   
01243   /*-- make an empty copy of this dataset, for eventual output --*/
01244   new_dset = EDIT_empty_copy( dset ) ;
01245 
01246   
01247   /*----- Record history of dataset -----*/
01248   if( commandline != NULL )
01249      tross_Append_History( new_dset , commandline ) ;
01250   sprintf (label, "Output prefix: %s", filename);
01251   tross_Append_History ( new_dset, label);
01252 
01253   
01254   iv = DSET_PRINCIPAL_VALUE(dset) ;
01255   output_datum = DSET_BRICK_TYPE(dset,iv) ;
01256   if( output_datum == MRI_byte ) output_datum = MRI_short ;
01257   
01258   
01259   ibuf[0] = output_datum ; 
01260   
01261   func_type = FUNC_FIM_TYPE;
01262   
01263   ierror = EDIT_dset_items( new_dset ,
01264                             ADN_prefix , filename ,
01265                             ADN_label1 , filename ,
01266                             ADN_self_name , filename ,
01267                             ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : 
01268                                                       GEN_FUNC_TYPE ,
01269                             ADN_func_type , func_type ,
01270                             ADN_nvals , FUNC_nvals[func_type] ,
01271                             ADN_datum_array , ibuf ,
01272                             ADN_ntt , 0 ,
01273                             ADN_malloc_type, DATABLOCK_MEM_MALLOC ,  
01274                             ADN_none ) ;
01275   
01276   if( ierror > 0 ){
01277     fprintf(stderr,
01278           "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01279     exit(1) ;
01280   }
01281   
01282   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01283     fprintf(stderr,
01284             "*** Output dataset file %s already exists--cannot continue!\a\n",
01285             new_dset->dblk->diskptr->header_name ) ;
01286     exit(1) ;
01287   }
01288 
01289   
01290   /*----- deleting exemplar dataset -----*/ 
01291   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
01292   
01293   
01294   /*----- allocate memory for output data -----*/
01295   vdif = (void *)  malloc( mri_datum_size(output_datum) * nxyz );
01296   if (vdif == NULL)   NLfit_error ("Unable to allocate memory for vdif");
01297 
01298   
01299   /*----- attach bricks to new data set -----*/
01300   mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0)); 
01301   
01302   
01303   /*----- convert data type to output specification -----*/
01304   fimfac = EDIT_coerce_autoscale_new (nxyz, 
01305                                       MRI_float, ffim, 
01306                                       output_datum, vdif);
01307   if (fimfac != 0.0)  fimfac = 1.0 / fimfac;
01308     
01309 
01310   /*----- write afni data set -----*/
01311   printf("--- Writing combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
01312   
01313   for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
01314   (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
01315   
01316   fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
01317   (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01318   
01319   THD_load_statistics( new_dset ) ;
01320   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01321   
01322   /*----- deallocate memory -----*/   
01323   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01324 
01325 }
01326 
01327 
01328 /*---------------------------------------------------------------------------*/
01329 /*
01330   Routine to write one bucket data set.
01331 */
01332 
01333 void write_bucket_data 
01334 (
01335   int q,                  /* number of parameters in the noise model */
01336   int p,                  /* number of parameters in the signal model */
01337   int  nxyz,              /* number of voxels in image */
01338   int  n,                 /* length of time series data */  
01339 
01340   float ** ncoef_vol,     /* volume of noise model parameters */
01341   float ** scoef_vol,     /* volume of signal model parameters */
01342 
01343   char * input_filename,
01344 
01345   NL_options * option_data     /* user input options */
01346 )
01347 
01348 {
01349   const float EPSILON = 1.0e-10;
01350 
01351   THD_3dim_dataset * old_dset = NULL;    /* prototype dataset */
01352   THD_3dim_dataset * new_dset = NULL;    /* output bucket dataset */
01353   char * output_prefix;     /* prefix name for bucket dataset */
01354   char * output_session;    /* directory for bucket dataset */
01355   int nbricks, ib;          /* number of sub-bricks in bucket dataset */
01356   short ** bar = NULL;      /* bar[ib] points to data for sub-brick #ib */
01357   float factor;             /* factor is new scale factor for sub-brick #ib */
01358   int brick_type;           /* indicates statistical type of sub-brick */
01359   int brick_coef;           /* regression coefficient index for sub-brick */
01360   char * brick_label;       /* character string label for sub-brick */
01361   int ierror;               /* number of errors in editing data */
01362   float * volume;           /* volume of floating point data */
01363   int dimension;            /* dimension of full model = p + q */
01364   char label[80];           /* label for output file history */ 
01365 
01366     
01367   /*----- initialize local variables -----*/
01368   nbricks = option_data->numbricks;
01369   output_prefix = option_data->bucket_filename;
01370   output_session = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
01371   strcpy (output_session, "./");
01372   dimension = p + q;
01373   
01374 
01375   /*----- allocate memory -----*/
01376   bar  = (short **) malloc (sizeof(short *) * nbricks);
01377   MTEST (bar);
01378 
01379  
01380   /*----- read first dataset -----*/
01381   old_dset = THD_open_one_dataset (input_filename);
01382   
01383 
01384   /*-- make an empty copy of this dataset, for eventual output --*/
01385   new_dset = EDIT_empty_copy (old_dset);
01386 
01387   
01388   /*----- Record history of dataset -----*/
01389   if( commandline != NULL )
01390      tross_Append_History( new_dset , commandline ) ;
01391   sprintf (label, "Output prefix: %s", output_prefix);
01392   tross_Append_History ( new_dset, label);
01393 
01394 
01395   /*----- Modify some structural properties.  Note that the nbricks
01396           just make empty sub-bricks, without any data attached. -----*/
01397   ierror = EDIT_dset_items (new_dset,
01398                             ADN_prefix,          output_prefix,
01399                             ADN_directory_name,  output_session,
01400                             ADN_type,            HEAD_FUNC_TYPE,
01401                             ADN_func_type,       FUNC_BUCK_TYPE,
01402                             ADN_ntt,             0,               /* no time */
01403                             ADN_nvals,           nbricks,
01404                             ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,  
01405                             ADN_none ) ;
01406   
01407   if( ierror > 0 )
01408     {
01409       fprintf(stderr, 
01410               "*** %d errors in attempting to create output dataset!\n", 
01411               ierror);
01412       exit(1);
01413     }
01414   
01415   if (THD_is_file(DSET_HEADNAME(new_dset))) 
01416     {
01417       fprintf(stderr,
01418               "*** Output dataset file %s already exists--cannot continue!\n",
01419               DSET_HEADNAME(new_dset));
01420       exit(1);
01421     }
01422   
01423 
01424   /*----- deleting exemplar dataset -----*/ 
01425   THD_delete_3dim_dataset( old_dset , False );  old_dset = NULL ;
01426   
01427 
01428   /*----- loop over new sub-brick index, attach data array with 
01429           EDIT_substitute_brick then put some strings into the labels and 
01430           keywords, and modify the sub-brick scaling factor -----*/
01431   for (ib = 0;  ib < nbricks;  ib++)
01432     {
01433       /*----- get information about this sub-brick -----*/
01434       brick_type  = option_data->brick_type[ib];
01435       brick_coef  = option_data->brick_coef[ib];
01436       brick_label = option_data->brick_label[ib];
01437 
01438       if (brick_type == FUNC_FIM_TYPE)
01439         {       
01440           if (brick_coef < q)
01441             volume = ncoef_vol[brick_coef];
01442           else if (brick_coef < p+q)
01443             volume = scoef_vol[brick_coef-q];
01444         }
01445 
01446       /*----- allocate memory for output sub-brick -----*/
01447       bar[ib]  = (short *) malloc (sizeof(short) * nxyz);
01448       MTEST (bar[ib]);
01449       factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
01450                                           MRI_short, bar[ib]);
01451 
01452       if (factor < EPSILON)  factor = 0.0;
01453       else factor = 1.0 / factor;
01454 
01455       /*----- edit the sub-brick -----*/
01456       EDIT_BRICK_LABEL (new_dset, ib, brick_label);
01457       EDIT_BRICK_FACTOR (new_dset, ib, factor);
01458 
01459       
01460       /*----- attach bar[ib] to be sub-brick #ib -----*/
01461       EDIT_substitute_brick (new_dset, ib, MRI_short, bar[ib]);
01462 
01463     }
01464 
01465 
01466   /*----- write bucket data set -----*/
01467   THD_load_statistics (new_dset);
01468   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01469   fprintf(stderr,"Wrote bucket dataset ");
01470   fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
01471 
01472   
01473   /*----- deallocate memory -----*/   
01474   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01475 
01476 }
01477 
01478 
01479 /*---------------------------------------------------------------------------*/
01480 /*
01481   Write out the parameter files.
01482 */
01483 
01484 void output_parameters
01485 (
01486   int r,                  /* number of parameters in the noise model */
01487   int p,                  /* number of parameters in the signal model */
01488   int  nxyz,              /* number of voxels in image */
01489   int  ts_length,         /* length of time series data */  
01490 
01491   float ** ncoef_vol,     /* volume of noise model parameters */
01492   float ** scoef_vol,     /* volume of signal model parameters */
01493 
01494   char * input_filename,     /* file name of prototype 3d+time dataset */
01495   char ** ncoef_filename,    /* file name for noise model parameters */
01496   char ** scoef_filename,    /* file name for signal model parameters */
01497 
01498   NL_options * option_data   /* user input options */
01499 )
01500 
01501 {
01502   int ip;                 /* parameter index */
01503   int dimension;          /* dimension of full model = r + p  */
01504   int numdof, dendof;     /* numerator and denominator degrees of freedom */
01505 
01506 
01507   dimension = r + p;
01508 
01509 
01510   /*----- write the bucket dataset -----*/
01511   if (option_data->numbricks > 0)
01512     write_bucket_data (r, p, nxyz, ts_length, ncoef_vol, scoef_vol,       
01513                        input_filename, option_data);
01514 
01515 
01516   /*----- write noise model parameters -----*/
01517   for (ip = 0;  ip < r;  ip++)
01518     {
01519       if (ncoef_filename[ip] != NULL)
01520         {
01521           write_afni_data (input_filename, nxyz, ncoef_filename[ip], 
01522                            ncoef_vol[ip]); 
01523         }
01524     }
01525 
01526   /*----- write signal model parameters -----*/
01527   for (ip = 0;  ip < p;  ip++)
01528     {
01529       if (scoef_filename[ip] != NULL)
01530         {
01531           write_afni_data (input_filename, nxyz, scoef_filename[ip], 
01532                            scoef_vol[ip]); 
01533         }
01534     }
01535 }
01536 
01537 
01538 /*---------------------------------------------------------------------------*/
01539 /*
01540   Release all allocated memory space.
01541 */
01542 
01543 void terminate_program 
01544 (
01545   int r,                       /* number of parameters in the noise model */
01546   int p,                       /* number of parameters in the signal model */
01547   int ts_length,               /* length of time series data */
01548   float *** x_array,           /* independent variable matrix */
01549   float ** ts_array,           /* input time series array */
01550   short *** d_array,           /* output time series volume data */
01551   char  ** nname,         /* noise model name */
01552   char  *** npname,       /* noise parameter labels */
01553   float ** min_nconstr,   /* min parameter constraints for noise model */
01554   float ** max_nconstr,   /* max parameter constraints for noise model */
01555   char ** sname,          /* signal model name */
01556   char *** spname,        /* signal parameter labels */
01557   float ** par_full,      /* estimated parameters for the full model */
01558   float ** min_sconstr,   /* min parameter constraints for signal model */
01559   float ** max_sconstr,   /* max parameter constraints for signal model */
01560   float *** ncoef_vol,        /* noise model parameters volume data */
01561   float *** scoef_vol,        /* signal model parameters volume data */
01562   char ** input_filename,        /* file name of prototype 3d+time dataset */
01563   char ** output_filename,       /* file name for output 3d+time dataset */
01564   char *** ncoef_filename,       /* file name for noise model parameters */
01565   char *** scoef_filename        /* file name for signal model parameters */
01566 )
01567  
01568 {
01569   int ip;                        /* parameter index */
01570   int it;                        /* time index */
01571 
01572 
01573   /*----- deallocate space for model names and parameters labels -----*/
01574   if (*nname != NULL)
01575     { free (*nname);  *nname = NULL; }
01576   if (*sname != NULL)
01577     { free (*sname);  *sname = NULL; }
01578   for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
01579     {
01580       if ((*npname)[ip] != NULL)
01581         { free ((*npname)[ip]);  (*npname)[ip] = NULL; }
01582       if ((*spname)[ip] != NULL)
01583         { free ((*spname)[ip]);  (*spname)[ip] = NULL; }
01584     }
01585 
01586   if (*npname != NULL)
01587     { free (*npname);  *npname = NULL; }
01588   if (*spname != NULL)
01589     { free (*spname);  *spname = NULL; }
01590 
01591 
01592   /*----- deallocate memory for parameter constraints -----*/
01593   if (*min_nconstr != NULL)  { free (*min_nconstr);  *min_nconstr = NULL; }
01594   if (*max_nconstr != NULL)  { free (*max_nconstr);  *max_nconstr = NULL; }
01595   if (*min_sconstr != NULL)  { free (*min_sconstr);  *min_sconstr = NULL; }
01596   if (*max_sconstr != NULL)  { free (*max_sconstr);  *max_sconstr = NULL; }
01597 
01598   
01599   /*----- deallocate memory space for filenames -----*/
01600   if (*input_filename != NULL)  
01601     { free (*input_filename);  *input_filename = NULL; }
01602   if (*output_filename != NULL)
01603     { free (*output_filename);   *output_filename = NULL; }
01604 
01605   if (*ncoef_filename != NULL)
01606     {
01607       for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
01608         {
01609           if ((*ncoef_filename)[ip] != NULL)
01610             { free ((*ncoef_filename)[ip]);  (*ncoef_filename)[ip] = NULL; } 
01611         }
01612       free (*ncoef_filename);  *ncoef_filename = NULL;
01613     }
01614 
01615   if (*scoef_filename != NULL)
01616     {
01617       for (ip = 0;  ip < MAX_PARAMETERS;  ip++)
01618         {
01619           if ((*scoef_filename)[ip] != NULL)
01620             { free ((*scoef_filename)[ip]);  (*scoef_filename)[ip] = NULL; } 
01621         }
01622       free (*scoef_filename);  *scoef_filename = NULL;
01623     }
01624 
01625 
01626   /*----- deallocate space for input time series -----*/
01627   if (*ts_array != NULL)    { free (*ts_array);    *ts_array = NULL; }
01628 
01629 
01630   /*----- deallocate space for independent variable matrix -----*/
01631   if (*x_array != NULL) 
01632     {
01633       for (it = 0;  it < ts_length;  it++)
01634         if ((*x_array)[it] != NULL)
01635           { free ((*x_array)[it]);  (*x_array)[it] = NULL; }
01636       free (*x_array);  *x_array = NULL; 
01637     }
01638  
01639 
01640   /*----- deallocate space for parameters -----*/
01641   if (*par_full != NULL)    { free (*par_full);    *par_full = NULL; }
01642 
01643 
01644   /*----- deallocate space for volume data -----*/
01645   if (*ncoef_vol != NULL)
01646     {
01647       for (ip = 0;  ip < r;  ip++)
01648         {
01649           if ((*ncoef_vol)[ip] != NULL)
01650             { free ((*ncoef_vol)[ip]);  (*ncoef_vol)[ip] = NULL; }
01651         }
01652       free (*ncoef_vol);   *ncoef_vol = NULL; 
01653     }
01654   
01655   if (*scoef_vol != NULL)
01656     {
01657       for (ip = 0;  ip < p;  ip++)
01658         {
01659           if ((*scoef_vol)[ip] != NULL)
01660             { free ((*scoef_vol)[ip]);  (*scoef_vol)[ip] = NULL; }
01661         }
01662       free (*scoef_vol);   *scoef_vol = NULL; 
01663     }
01664   
01665 }
01666 
01667 
01668 /*---------------------------------------------------------------------------*/
01669 
01670 int main 
01671 (
01672   int argc,                /* number of input arguments */
01673   char ** argv             /* array of input arguments */ 
01674 )
01675 
01676 {
01677   NL_options option_data;  /* bucket dataset options */
01678 
01679   /*----- declare time series variables -----*/
01680   int ts_length;                       /* length of time series data */
01681   float ** x_array = NULL;             /* independent variable matrix */
01682   float * ts_array = NULL;             /* generated time series array */
01683   int nxyz;                            /* number of voxels in image */
01684   int iv;                              /* voxel counter */
01685   int nvoxel;                          /* screen output for voxel #nvoxel */
01686   short ** d_array = NULL;             /* output time series volume */
01687 
01688   /*----- declare reduced (noise) model variables -----*/
01689   char * nname = NULL;         /* noise model name */
01690   vfp nmodel;                  /* pointer to noise model */
01691   int r;                       /* number of parameters in the noise model */
01692   char ** npname = NULL;       /* noise parameter labels */
01693   float * min_nconstr = NULL;  /* min parameter constraints for noise model */
01694   float * max_nconstr = NULL;  /* max parameter constraints for noise model */
01695 
01696   /*----- declare full (signal+noise) model variables -----*/
01697   char * sname = NULL;         /* signal model name */
01698   vfp smodel;                  /* pointer to signal model */
01699   int p;                       /* number of parameters in the signal model */
01700   float * par_full = NULL;     /* estimated parameters for the full model */
01701   char ** spname = NULL;       /* signal parameter labels */
01702   float * min_sconstr = NULL;  /* min parameter constraints for signal model */
01703   float * max_sconstr = NULL;  /* max parameter constraints for signal model */
01704 
01705   /*----- declare statistical test variables -----*/
01706   float sigma;             /* std. dev. of additive Gaussian noise */
01707 
01708   /*----- declare output volume data -----*/
01709   float ** ncoef_vol = NULL;    /* noise model parameters volume data */
01710   float ** scoef_vol = NULL;    /* signal model parameters volume data */
01711 
01712   /*----- declare file name variables -----*/
01713   char * input_filename = NULL;    /* file name of prototype 3d+time dataset */
01714   char * output_filename = NULL;   /* file name for output 3d+time dataset */
01715   char ** ncoef_filename = NULL;   /* file name for noise model parameters */
01716   char ** scoef_filename = NULL;   /* file name for signal model parameters */  
01717   
01718   /*----- Identify software -----*/
01719   printf ("\n\n");
01720   printf ("Program: %s \n", PROGRAM_NAME);
01721   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
01722   printf ("Date:    %s \n", PROGRAM_DATE);
01723   printf ("\n");
01724 
01725    
01726   /*----- program initialization -----*/
01727   initialize_program (argc, argv, 
01728                       &nname, &sname, &nmodel, &smodel, 
01729                       &r, &p, &npname, &spname,
01730                       &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr,
01731                       &sigma, &nvoxel,
01732                       &input_filename, &output_filename, &ncoef_filename,
01733                       &scoef_filename,
01734                       &nxyz, &ts_length, &x_array, &ts_array, &d_array, 
01735                       &par_full, 
01736                       &ncoef_vol, &scoef_vol, &option_data);
01737 
01738 
01739   /*----- loop over voxels in the data set -----*/
01740   for (iv = 0;  iv < nxyz;  iv++)
01741     {
01742 
01743       /*----- generate artificial time series for voxel iv -----*/
01744       generate_ts_array (nmodel, smodel, r, p,
01745                          min_nconstr, max_nconstr, min_sconstr, max_sconstr,
01746                          ts_length, x_array, sigma,
01747                          par_full, ts_array);
01748 
01749 
01750       /*----- save time series into data set -----*/
01751       save_ts_array (iv, ts_length, ts_array, d_array);
01752 
01753 
01754       /*----- write time series for this voxel -----*/
01755       if (iv == nvoxel-1)
01756         write_ts_array (iv, ts_length, ts_array, d_array);
01757       
01758 
01759       /*----- save parameters for this voxel into volume data -----*/
01760       save_parameters (iv, r, p, par_full, ncoef_vol, scoef_vol);
01761       
01762 
01763       /*----- write parameters for this voxel -----*/
01764       if (iv == nvoxel-1)
01765         write_parameters (iv, nname, sname, r, p, npname, spname, par_full);
01766 
01767     }
01768 
01769 
01770   /*----- save time series into data set -----*/
01771   output_ts_array (d_array, ts_length, input_filename, output_filename);
01772 
01773 
01774   /*----- output the parameter files -----*/
01775   output_parameters (r, p, nxyz, ts_length, ncoef_vol, scoef_vol,
01776                   input_filename,
01777                   ncoef_filename, scoef_filename, &option_data);
01778 
01779                  
01780   /*----- end of program -----*/
01781   terminate_program (r, p, ts_length, &x_array, &ts_array, &d_array, 
01782                      &nname, &npname, &min_nconstr, &max_nconstr, 
01783                      &sname, &spname, &par_full, &min_sconstr, &max_sconstr, 
01784                      &ncoef_vol, &scoef_vol,
01785                      &input_filename, &output_filename, 
01786                      &ncoef_filename, &scoef_filename);
01787 }
01788 
01789 
01790 
01791