/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1994-2001, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ #ifdef USE_SUNPERF /** for Solaris **/ # include #endif /* This program calculates a nonlinear regression for each voxel of the input AFNI 3d+time data set. The nonlinear regression is calculated by means of a least squares fit to the signal plus noise models which are specified by the user. File: 3dNLfim.c Author: B. Douglas Ward Date: 19 June 1997 Mod: Added external time reference capability (Rongyan Zhang) Date: 10 August 1997 Mod: Added option for absolute noise parameter constraints. Date: 22 August 1997 Mod: Added options for percent signal change above baseline, and calculation of area under the signal above baseline. Date: 26 November 1997 Mod: Print error message if user fails to specify the signal model or the noise model. Date: 26 December 1997 Mod: Extensive changes required to implement the 'bucket' dataset. Date: 14 January 1998 Mod: Added the -inTR option. 22 July 1998 -- RWCox Mod: Incorporated THD_extract_series routine. Date: 19 April 1999 Mod: Added -sfit and -snfit options to write out the signal and the signal+noise model time series fit for each voxel to a 3d+time dataset. Date: 08 July 1999 Mod: Added novar flag to eliminate unnecessary calculations. Date: 13 July 1999 Mod: Added changes for incorporating History notes. Date: 09 September 1999 Mod: Adjust F-statistics if parameter constraints force a parameter to be a constant. Date: 08 February 2000 Mod: Changes for output of R^2 (coefficient of multiple determination), and standard deviation of residuals from full model fit. Added global variable calc_tstats. Also, added screen display of p-values. Date: 10 May 2000 Mod: Corrected error in initializing of output data type (MRI_short) in routine write_3dtime. Date: 17 May 2000 Mod: Added -mask option. (Adapted from: 3dpc.c) Date: 18 May 2000 Mod: Changed "return" at end of program to exit(0). Also, increased maximum number of model parameters. Date: 08 August 2001 Mod: Added call to AFNI_logger. Date: 15 August 2001 Mod: Changes to allow -jobs option. Date: 07 May 2003 - RWCox. Mod: Added options -aux_name, -aux_fval and -voxel_count. Date: 25 Jan 2006 [rickr] Mod: Removed options -aux_name and -aux_fval, and the globals require linking to afni, too. Date: 30 Jan 2006 [rickr] Mod: Added NEWUOA stuff (mostly to simplex.c, actually). Date: 20 Jul 2006 [RWCox] Mod: Copied memmap code from 3dDeconvolve.c Date: 24 Oct 2006 [DRG] Mod: Limit reports to nth voxels via progress option Date: 25 Oct 2006 [DRG] Mod: Limit g_voxel_count reports to every 10th voxel. Date: 26 Oct 2006 [rickr] */ /*---------------------------------------------------------------------------*/ #define PROGRAM_NAME "3dNLfim" /* name of this program */ #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */ #define PROGRAM_INITIAL "19 June 1997" /* date of initial program release */ #define PROGRAM_LATEST "24 Oct 2006 - DRG" /* date of latest program revision */ /*---------------------------------------------------------------------------*/ #define DEFAULT_NRAND 19999 #define DEFAULT_NBEST 9 #define DEFAULT_FDISP 999.0 #define DEFAULT_PROGRESS 10000 #include #include #include #include "mrilib.h" /*---------------------------------------------------------------------------*/ /*--------- Global variables for multiple process execution - RWCox. --------*/ /*--------- All names start with "proc_", so search for that string. --------*/ #if !defined(DONT_USE_SHM) && !defined(DONT_USE_FORK) && !defined(CYGWIN) # include "thd_iochan.h" /* prototypes for shm stuff */ # define PROC_MAX 32 /* max num processes */ static int proc_numjob = 1 ; /* num processes */ static pid_t proc_pid[PROC_MAX] ; /* IDs of processes */ static int proc_shmid = 0 ; /* shared memory ID */ static char *proc_shmptr = NULL; /* pointer to shared memory */ static long long proc_shmsize = 0 ; /* total size of shared memory */ static int proc_shm_arnum = 0 ; /* num arrays in shared memory */ static float ***proc_shm_ar = NULL; /* *proc_shm_ar[i] = ptr to array #i */ static int *proc_shm_arsiz = NULL; /* proc_shm_arsiz[i] = floats in #i */ static int proc_vox_bot[PROC_MAX] ; /* 1st voxel to use in each process */ static int proc_vox_top[PROC_MAX] ; /* last voxel (+1) in each process */ static int proc_ind = 0 ; /* index of THIS job */ #else /* can't use multiple processes */ # define proc_numjob 1 /* flag that only 1 process is in use */ # define proc_ind 0 /* index of THIS job */ #endif /*---------------------------------------------------------------------------*/ #include "matrix.h" #include "simplex.h" #include "NLfit.h" #include "matrix.c" #include "simplex.c" /* 20 Jul 2006 - now includes NEWUOA variables [N_*] */ #include "NLfit.c" typedef struct NL_options { char * bucket_filename; /* file name for bucket dataset */ int numbricks; /* number of sub-bricks in bucket dataset */ int * brick_type; /* indicates type of sub-brick */ int * brick_coef; /* regression coefficient number for sub-brick*/ char ** brick_label; /* character string label for sub-brick */ } NL_options; /*---------------------------------------------------------------------------*/ /* Global data */ /***** 22 July 1998 -- RWCox: Modified to allow DELT to be set from the TR of the input file *****/ static float DELT = 1.0; /* default */ static int inTR = 0 ; /* set to 1 if -inTR option is used */ static float dsTR = 0.0 ; /* TR of the input file */ static int g_voxel_count = 0; /* display current voxel counter */ /* 25 Jan 2006 [rickr] */ static char * commandline = NULL ; /* command line for history notes */ static byte * mask_vol = NULL; /* mask volume */ static int mask_nvox = 0; /* number of voxels in mask volume */ static int output_datum = ILLEGAL_TYPE ; /*---------------------------------------------------------------------------*/ /* Routine to display 3dNLfim help menu. */ void display_help_menu() { printf ( "This program calculates a nonlinear regression for each voxel of the \n" "input AFNI 3d+time data set. The nonlinear regression is calculated \n" "by means of a least squares fit to the signal plus noise models which \n" "are specified by the user. \n" " \n" "Usage: \n" "3dNLfim \n" "-input fname fname = filename of 3d + time data file for input \n" "[-mask mset] Use the 0 sub-brick of dataset 'mset' as a mask \n" " to indicate which voxels to analyze (a sub-brick \n" " selector is allowed) [default = use all voxels] \n" "[-ignore num] num = skip this number of initial images in the \n" " time series for regresion analysis; default = 3 \n" "[-inTR] set delt = TR of the input 3d+time dataset \n" " [The default is to compute with delt = 1.0 ] \n" " [The model functions are calculated using a \n" " time grid of: 0, delt, 2*delt, 3*delt, ... ] \n" "[-time fname] fname = ASCII file containing each time point \n" " in the time series. Defaults to even spacing \n" " given by TR (this option overrides -inTR). \n" "-signal slabel slabel = name of (non-linear) signal model \n" "-noise nlabel nlabel = name of (linear) noise model \n" "-sconstr k c d constraints for kth signal parameter: \n" " c <= gs[k] <= d \n" "-nconstr k c d constraints for kth noise parameter: \n" " c+b[k] <= gn[k] <= d+b[k] \n" "[-nabs] use absolute constraints for noise parameters: \n" " c <= gn[k] <= d [default=relative, as above] \n" "[-nrand n] n = number of random test points [default=%d] \n" "[-nbest b] b = use b best test points to start [default=%d] \n" "[-rmsmin r] r = minimum rms error to reject reduced model \n" "[-fdisp fval] display (to screen) results for those voxels \n" " whose f-statistic is > fval [default=%.1f] \n" "[-progress ival] display (to screen) results for those voxels \n" " every ival number of voxels \n" "[-voxel_count] display (to screen) the current voxel index \n" " \n" "--- These options choose the least-square minimization algorithm --- \n" " \n" "[-SIMPLEX] use Nelder-Mead simplex method [default] \n" "[-POWELL] use Powell's NEWUOA method instead of the \n" " Nelder-Mead simplex method to find the \n" " nonlinear least-squares solution \n" " [slower; usually more accurate, but not always!] \n" "[-BOTH] use both Powell's and Nelder-Mead method \n" " [slowest, but should be most accurate] \n" " \n" "--- These options generate individual AFNI 2 sub-brick datasets --- \n" " \n" "[-freg fname] perform f-test for significance of the regression; \n" " output 'fift' is written to prefix filename fname\n" "[-frsqr fname] calculate R^2 (coef. of multiple determination); \n" " store along with f-test for regression; \n" " output 'fift' is written to prefix filename fname\n" "[-fsmax fname] estimate signed maximum of signal; store along \n" " with f-test for regression; output 'fift' is \n" " written to prefix filename fname \n" "[-ftmax fname] estimate time of signed maximum; store along \n" " with f-test for regression; output 'fift' is \n" " written to prefix filename fname \n" "[-fpsmax fname] calculate (signed) maximum percentage change of \n" " signal from baseline; output 'fift' is \n" " written to prefix filename fname \n" "[-farea fname] calculate area between signal and baseline; store \n" " with f-test for regression; output 'fift' is \n" " written to prefix filename fname \n" "[-fparea fname] percentage area of signal relative to baseline; \n" " store with f-test for regression; output 'fift' \n" " is written to prefix filename fname \n" "[-fscoef k fname] estimate kth signal parameter gs[k]; store along \n" " with f-test for regression; output 'fift' is \n" " written to prefix filename fname \n" "[-fncoef k fname] estimate kth noise parameter gn[k]; store along \n" " with f-test for regression; output 'fift' is \n" " written to prefix filename fname \n" "[-tscoef k fname] perform t-test for significance of the kth signal \n" " parameter gs[k]; output 'fitt' is written \n" " to prefix filename fname \n" "[-tncoef k fname] perform t-test for significance of the kth noise \n" " parameter gn[k]; output 'fitt' is written \n" " to prefix filename fname \n" " \n" "--- These options generate one AFNI 'bucket' type dataset --- \n" " \n" "[-bucket n prefixname] create one AFNI 'bucket' dataset containing \n" " n sub-bricks; n=0 creates default output; \n" " output 'bucket' is written to prefixname \n" "The mth sub-brick will contain: \n" "[-brick m scoef k label] kth signal parameter regression coefficient\n" "[-brick m ncoef k label] kth noise parameter regression coefficient \n" "[-brick m tmax label] time at max. abs. value of signal \n" "[-brick m smax label] signed max. value of signal \n" "[-brick m psmax label] signed max. value of signal as percent \n" " above baseline level \n" "[-brick m area label] area between signal and baseline \n" "[-brick m parea label] signed area between signal and baseline \n" " as percent of baseline area \n" "[-brick m tscoef k label] t-stat for kth signal parameter coefficient\n" "[-brick m tncoef k label] t-stat for kth noise parameter coefficient \n" "[-brick m resid label] std. dev. of the full model fit residuals \n" "[-brick m rsqr label] R^2 (coefficient of multiple determination)\n" "[-brick m fstat label] F-stat for significance of the regression \n" " \n" " --- These options write time series fit for --- \n" " --- each voxel to an AFNI 3d+time dataset --- \n" " \n" "[-sfit fname] fname = prefix for output 3d+time signal model fit \n" "[-snfit fname] fname = prefix for output 3d+time signal+noise fit \n" " \n" , DEFAULT_NRAND , DEFAULT_NBEST , DEFAULT_FDISP ); #ifdef PROC_MAX printf( "\n" " -jobs J Run the program with 'J' jobs (sub-processes).\n" " On a multi-CPU machine, this can speed the\n" " program up considerably. On a single CPU\n" " machine, using this option is silly.\n" " J should be a number from 1 up to the\n" " number of CPU sharing memory on the system.\n" " J=1 is normal (single process) operation.\n" " The maximum allowed value of J is %d.\n" " * For more information on parallelizing, see\n" " http://afni.nimh.nih.gov/afni/doc/misc/parallize.html\n" " * Use -mask to get more speed; cf. 3dAutomask.\n" , PROC_MAX ) ; #endif exit(0); } /*---------------------------------------------------------------------------*/ /** macro to test a malloc-ed pointer for validity **/ #define MTEST(ptr) \ if((ptr)==NULL) \ ( NLfit_error ("Cannot allocate memory") ) /*---------------------------------------------------------------------------*/ /* Routine to initialize the input options. */ void initialize_options ( int * ignore, /* delete this number of points from time series */ vfp * nmodel, /* pointer to noise model */ vfp * smodel, /* pointer to signal model */ int * r, /* number of parameters in the noise model */ int * p, /* number of parameters in the signal model */ char *** npname, /* noise parameter names */ char *** spname, /* signal parameter names */ float ** min_nconstr, /* minimum parameter constraints for noise model */ float ** max_nconstr, /* maximum parameter constraints for noise model */ float ** min_sconstr, /* minimum parameter constraints for signal model */ float ** max_sconstr, /* maximum parameter constraints for signal model */ int * nabs, /* use absolute constraints for noise parameters */ int * nrand, /* number of random vectors to generate */ int * nbest, /* number of random vectors to keep */ float * rms_min, /* minimum rms error to reject reduced model */ float * fdisp, /* minimum f-statistic for display */ int *progress, /* nth voxel to show report */ char ** input_filename, /* file name of input 3d+time dataset */ char ** tfilename, /* file name for time point series */ char ** freg_filename, /* file name for regression f-statistics */ char ** frsqr_filename, /* file name for R^2 statistics */ char *** fncoef_filename, /* file name for noise model parameters */ char *** fscoef_filename, /* file name for signal model parameters */ char *** tncoef_filename, /* file name for noise model t-statistics */ char *** tscoef_filename, /* file name for signal model t-statistics */ char ** sfit_filename, /* file name for 3d+time fitted signal model */ char ** snfit_filename, /* file name for 3d+time fitted signal+noise */ NL_options * option_data /* bucket dataset options */ ) { int ip; /* parameter index */ /*----- initialize default values -----*/ *ignore = 3; *nabs = 0; *nrand = DEFAULT_NRAND; *nbest = DEFAULT_NBEST; *rms_min = 0.0; *fdisp = DEFAULT_FDISP; *progress = DEFAULT_PROGRESS; *smodel = NULL; *nmodel = NULL; *r = -1; *p = -1; /*----- allocate memory for noise parameter names -----*/ *npname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS); MTEST (*npname); for (ip = 0; ip < MAX_PARAMETERS; ip++) { (*npname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST ((*npname)[ip]); } /*----- allocate memory for signal parameter names -----*/ *spname = (char **) malloc (sizeof(char *) * MAX_PARAMETERS); MTEST (*spname); for (ip = 0; ip < MAX_PARAMETERS; ip++) { (*spname)[ip] = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST ((*spname)[ip]); } /*----- allocate memory for parameter constraints -----*/ *min_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS); MTEST (*min_nconstr); *max_nconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS); MTEST (*max_nconstr); *min_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS); MTEST (*min_sconstr); *max_sconstr = (float *) malloc (sizeof(float) * MAX_PARAMETERS); MTEST (*max_sconstr); /*----- allocate memory space and initialize pointers for filenames -----*/ *input_filename = NULL; *tfilename = NULL; *freg_filename = NULL; *frsqr_filename = NULL; *fncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS); MTEST (*fncoef_filename); *fscoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS); MTEST (*fscoef_filename); *tncoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS); MTEST (*tncoef_filename); *tscoef_filename = (char **) malloc (sizeof(char *) * MAX_PARAMETERS); MTEST (*tscoef_filename); for (ip = 0; ip < MAX_PARAMETERS; ip++) { (*fncoef_filename)[ip] = NULL; (*fscoef_filename)[ip] = NULL; (*tncoef_filename)[ip] = NULL; (*tscoef_filename)[ip] = NULL; } *sfit_filename = NULL; *snfit_filename = NULL; /*----- initialize bucket dataset options -----*/ option_data->bucket_filename = NULL; option_data->numbricks = -1; option_data->brick_type = NULL; option_data->brick_coef = NULL; option_data->brick_label = NULL; } /*---------------------------------------------------------------------------*/ /* Routine to get user specified input options. */ void get_options ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ int * ignore, /* delete this number of points from time series */ char ** nname, /* name of noise model */ char ** sname, /* name of signal model */ vfp * nmodel, /* pointer to noise model */ vfp * smodel, /* pointer to signal model */ int * r, /* number of parameters in the noise model */ int * p, /* number of parameters in the signal model */ char *** npname, /* noise parameter names */ char *** spname, /* signal parameter names */ float ** min_nconstr, /* minimum parameter constraints for noise model */ float ** max_nconstr, /* maximum parameter constraints for noise model */ float ** min_sconstr, /* minimum parameter constraints for signal model */ float ** max_sconstr, /* maximum parameter constraints for signal model */ int * nabs, /* use absolute constraints for noise parameters */ int * nrand, /* number of random vectors to generate */ int * nbest, /* number of random vectors to keep */ float * rms_min, /* minimum rms error to reject reduced model */ float * fdisp, /* minimum f-statistic for display */ int * progress, /* progress report every nth voxel */ char ** input_filename, /* file name of input 3d+time dataset */ char ** tfilename, /* file name for time point series */ char ** freg_filename, /* file name for regression f-statistics */ char ** frsqr_filename, /* file name for R^2 statistics */ char ** fsmax_filename, /* file name for signal signed maximum */ char ** ftmax_filename, /* file name for time of signed maximum */ char ** fpmax_filename, /* file name for max. percentage change */ char ** farea_filename, /* file name for area under the signal */ char ** fparea_filename, /* file name for percent area under the signal */ char *** fncoef_filename, /* file name for noise model parameters */ char *** fscoef_filename, /* file name for signal model parameters */ char *** tncoef_filename, /* file name for noise model t-statistics */ char *** tscoef_filename, /* file name for signal model t-statistics */ char ** sfit_filename, /* file name for 3d+time fitted signal model */ char ** snfit_filename, /* file name for 3d+time fitted signal+noise */ THD_3dim_dataset ** dset_time, /* input 3d+time data set */ int * nxyz, /* number of voxels in image */ int * ts_length, /* length of time series data */ NL_options * option_data /* bucket dataset options */ ) { const MAX_BRICKS = 100; /* max. number of bricks in the bucket */ int nopt = 1; /* input option argument counter */ int ival, index; /* integer input */ float fval; /* float input */ char message[MAX_NAME_LENGTH]; /* error message */ int ok; /* boolean for specified model exists */ NLFIT_MODEL_array * model_array = NULL; /* array of SO models */ int im; /* model index */ int ibrick; /* sub-brick index */ int nbricks; /* number of bricks in the bucket */ /*----- does user request help menu? -----*/ if (argc < 2 || strcmp(argv[1], "-help") == 0) display_help_menu(); /*----- add to program log -----*/ AFNI_logger (PROGRAM_NAME,argc,argv); /*----- initialize the model array -----*/ model_array = NLFIT_get_many_MODELs (); if ((model_array == NULL) || (model_array->num == 0)) NLfit_error ("Unable to locate any models"); /*----- initialize the input options -----*/ initialize_options (ignore, nmodel, smodel, r, p, npname, spname, min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs, nrand, nbest, rms_min, fdisp, progress, input_filename, tfilename, freg_filename, frsqr_filename, fncoef_filename, fscoef_filename, tncoef_filename, tscoef_filename, sfit_filename, snfit_filename, option_data); /*----- main loop over input options -----*/ while (nopt < argc ) { /*----- -input filename -----*/ if (strcmp(argv[nopt], "-input") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -input "); *input_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*input_filename); strcpy (*input_filename, argv[nopt]); /* open input dataset - was THD_open_one_dataset, allow sub-bricks */ *dset_time = THD_open_dataset (*input_filename); if ((*dset_time) == NULL) { sprintf (message, "Unable to open data file: %s", *input_filename); NLfit_error (message); } DSET_UNMSEC( *dset_time ) ; /* RWCox */ DSET_load(*dset_time) ; CHECK_LOAD_ERROR(*dset_time) ; *nxyz = (*dset_time)->dblk->diskptr->dimsizes[0] * (*dset_time)->dblk->diskptr->dimsizes[1] * (*dset_time)->dblk->diskptr->dimsizes[2] ; *ts_length = DSET_NUM_TIMES(*dset_time); dsTR = DSET_TIMESTEP(*dset_time) ; if(output_datum==ILLEGAL_TYPE) { /* if output_datum type is not specified by user*/ output_datum = DSET_BRICK_TYPE(*dset_time,0); /* get datum type from dataset */ } nopt++; continue; } /**** -datum type ****/ if( strcmp(argv[nopt],"-datum") == 0 ){ if( ++nopt >= argc ) NLfit_error("need an argument after -datum!\n") ; if( strcmp(argv[nopt],"short") == 0 ){ output_datum = MRI_short ; } else if( strcmp(argv[nopt],"float") == 0 ){ output_datum = MRI_float ; } else { ERROR_exit("-datum of type '%s' not supported in 3dNLfim!\n",argv[nopt]) ; } nopt++ ; continue ; /* go to next arg */ } /**** -mask mset [18 May 2000] ****/ if( strcmp(argv[nopt],"-mask") == 0 ){ THD_3dim_dataset * mset ; int ii,mc ; nopt++ ; if (nopt >= argc) NLfit_error ("need argument after -mask!") ; mset = THD_open_dataset( argv[nopt] ) ; if (mset == NULL) NLfit_error ("can't open -mask dataset!") ; mask_vol = THD_makemask( mset , 0 , 1.0,0.0 ) ; mask_nvox = DSET_NVOX(mset) ; DSET_delete(mset) ; if (mask_vol == NULL ) NLfit_error ("can't use -mask dataset!") ; for( ii=mc=0 ; ii < mask_nvox ; ii++ ) if (mask_vol[ii]) mc++ ; if (mc == 0) NLfit_error ("mask is all zeros!") ; printf("++ %d voxels in mask %s\n",mc,argv[nopt]) ; nopt++ ; continue ; } /*----- 22 July 1998: the -inTR option -----*/ if( strcmp(argv[nopt],"-inTR") == 0 ){ inTR = 1 ; nopt++ ; continue ; } /*----- -ignore num -----*/ if (strcmp(argv[nopt], "-ignore") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -ignore "); sscanf (argv[nopt], "%d", &ival); if (ival < 0) NLfit_error ("illegal argument after -ignore "); *ignore = ival; nopt++; continue; } /*----- -time filename -----*/ if (strcmp(argv[nopt], "-time") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -time "); *tfilename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*tfilename); strcpy (*tfilename, argv[nopt]); nopt++; continue; } /*----- -signal slabel -----*/ if (strcmp(argv[nopt], "-signal") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -signal "); *sname = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*sname); strcpy (*sname, argv[nopt]); initialize_signal_model (model_array, *sname, smodel, p, *spname, *min_sconstr, *max_sconstr); nopt++; continue; } /*----- -noise nlabel -----*/ if (strcmp(argv[nopt], "-noise") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -noise "); *nname = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*nname); strcpy (*nname, argv[nopt]); initialize_noise_model (model_array, *nname, nmodel, r, *npname, *min_nconstr, *max_nconstr); nopt++; continue; } /*----- check that user has specified the signal and noise models -----*/ if ((*smodel) == NULL) NLfit_error ("Must specify signal model"); if ((*nmodel) == NULL) NLfit_error ("Must specify noise model"); /*----- -sconstr k min max -----*/ if (strcmp(argv[nopt], "-sconstr") == 0) { nopt++; if (nopt+2 >= argc) NLfit_error("need 3 arguments after -sconstr "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival >= *p)) NLfit_error ("illegal argument after -sconstr "); index = ival; nopt++; sscanf (argv[nopt], "%f", &fval); (*min_sconstr)[index] = fval; nopt++; sscanf (argv[nopt], "%f", &fval); (*max_sconstr)[index] = fval; nopt++; continue; } /*----- -nconstr k min max -----*/ if (strcmp(argv[nopt], "-nconstr") == 0) { nopt++; if (nopt+2 >= argc) NLfit_error("need 3 arguments after -nconstr "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival >= *r)) NLfit_error ("illegal argument after -nconstr "); index = ival; nopt++; sscanf (argv[nopt], "%f", &fval); (*min_nconstr)[index] = fval; nopt++; sscanf (argv[nopt], "%f", &fval); (*max_nconstr)[index] = fval; nopt++; continue; } /*----- -nabs -----*/ if (strcmp(argv[nopt], "-nabs") == 0) { *nabs = 1; nopt++; continue; } /*----- -POWELL -----*/ if( strcmp(argv[nopt],"-POWELL") == 0 ){ /* 20 Jul 2006 */ N_newuoa = 1 ; nopt++ ; continue ; } if( strcmp(argv[nopt],"-SIMPLEX") == 0 ){ N_newuoa = 0 ; nopt++ ; continue ; } if( strcmp(argv[nopt],"-BOTH") == 0 ){ N_newuoa = 2 ; nopt++ ; continue ; } /*----- -nrand n -----*/ if (strcmp(argv[nopt], "-nrand") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -nrand "); sscanf (argv[nopt], "%d", &ival); if (ival <= 0) NLfit_error ("illegal argument after -nrand "); *nrand = ival; nopt++; continue; } /*----- -nbest n -----*/ if (strcmp(argv[nopt], "-nbest") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -nbest "); sscanf (argv[nopt], "%d", &ival); if (ival <= 0) NLfit_error ("illegal argument after -nbest "); *nbest = ival; nopt++; continue; } /*----- -rmsmin r -----*/ if (strcmp(argv[nopt], "-rmsmin") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -rmsmin "); sscanf (argv[nopt], "%f", &fval); if (fval < 0.0) NLfit_error ("illegal argument after -rmsmin "); *rms_min = fval; nopt++; continue; } /*----- -fdisp fval -----*/ if (strcmp(argv[nopt], "-fdisp") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -fdisp "); sscanf (argv[nopt], "%f", &fval); *fdisp = fval; nopt++; continue; } /*----- -progress ival -----*/ if (strcmp(argv[nopt], "-progress") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -progress "); sscanf (argv[nopt], "%d", &ival); if (ival < 1) NLfit_error("illegal argument after -progress "); *progress = ival; nopt++; continue; } /*----- -voxel_count -----*/ if (strcmp(argv[nopt], "-voxel_count") == 0) { nopt++; g_voxel_count = 1; continue; } /*----- -freg filename -----*/ if (strcmp(argv[nopt], "-freg") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -freg "); *freg_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*freg_filename); strcpy (*freg_filename, argv[nopt]); nopt++; continue; } /*----- -frsqr filename -----*/ if (strcmp(argv[nopt], "-frsqr") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -frsqr "); *frsqr_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*frsqr_filename); strcpy (*frsqr_filename, argv[nopt]); nopt++; continue; } /*----- -fsmax filename -----*/ if (strcmp(argv[nopt], "-fsmax") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -fsmax "); *fsmax_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*fsmax_filename); strcpy (*fsmax_filename, argv[nopt]); nopt++; continue; } /*----- -ftmax filename -----*/ if (strcmp(argv[nopt], "-ftmax") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -ftmax "); *ftmax_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*ftmax_filename); strcpy (*ftmax_filename, argv[nopt]); nopt++; continue; } /*----- -fpmax filename -----*/ if (strcmp(argv[nopt], "-fpmax") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -fpmax "); *fpmax_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*fpmax_filename); strcpy (*fpmax_filename, argv[nopt]); nopt++; continue; } /*----- -farea filename -----*/ if (strcmp(argv[nopt], "-farea") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -farea "); *farea_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*farea_filename); strcpy (*farea_filename, argv[nopt]); nopt++; continue; } /*----- -fparea filename -----*/ if (strcmp(argv[nopt], "-fparea") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -fparea "); *fparea_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*fparea_filename); strcpy (*fparea_filename, argv[nopt]); nopt++; continue; } /*----- -fscoef k filename -----*/ if (strcmp(argv[nopt], "-fscoef") == 0) { nopt++; if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -fscoef "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival >= *p)) NLfit_error ("illegal argument after -fscoef "); index = ival; nopt++; (*fscoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST ((*fscoef_filename)[index]); strcpy ((*fscoef_filename)[index], argv[nopt]); nopt++; continue; } /*----- -fncoef k filename -----*/ if (strcmp(argv[nopt], "-fncoef") == 0) { nopt++; if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -fncoef "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival >= *r)) NLfit_error ("illegal argument after -fncoef "); index = ival; nopt++; (*fncoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST ((*fncoef_filename)[index]); strcpy ((*fncoef_filename)[index], argv[nopt]); nopt++; continue; } /*----- -tscoef k filename -----*/ if (strcmp(argv[nopt], "-tscoef") == 0) { nopt++; if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -tscoef "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival >= *p)) NLfit_error ("illegal argument after -tscoef "); index = ival; nopt++; (*tscoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST ((*tscoef_filename)[index]); strcpy ((*tscoef_filename)[index], argv[nopt]); calc_tstats = 1; nopt++; continue; } /*----- -tncoef k filename -----*/ if (strcmp(argv[nopt], "-tncoef") == 0) { nopt++; if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -tncoef "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival >= *r)) NLfit_error ("illegal argument after -tncoef "); index = ival; nopt++; (*tncoef_filename)[index] = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST ((*tncoef_filename)[index]); strcpy ((*tncoef_filename)[index], argv[nopt]); calc_tstats = 1; nopt++; continue; } /*----- -bucket n prefixname -----*/ if (strcmp(argv[nopt], "-bucket") == 0) { nopt++; if (nopt+1 >= argc) NLfit_error ("need 2 arguments after -bucket "); sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival > MAX_BRICKS)) NLfit_error ("illegal argument after -bucket "); nopt++; option_data->bucket_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (option_data->bucket_filename); strcpy (option_data->bucket_filename, argv[nopt]); /*----- set number of sub-bricks in the bucket -----*/ if (ival == 0) nbricks = (*p) + (*r) + 8; else nbricks = ival; option_data->numbricks = nbricks; /*----- allocate memory and initialize bucket dataset options -----*/ option_data->brick_type = malloc (sizeof(int) * nbricks); MTEST (option_data->brick_type); option_data->brick_coef = malloc (sizeof(int) * nbricks); MTEST (option_data->brick_coef); option_data->brick_label = malloc (sizeof(char *) * nbricks); MTEST (option_data->brick_label); for (ibrick = 0; ibrick < nbricks; ibrick++) { option_data->brick_type[ibrick] = -1; option_data->brick_coef[ibrick] = -1; option_data->brick_label[ibrick] = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (option_data->brick_label[ibrick]); } if (ival == 0) /*----- throw (almost) everything into the bucket -----*/ { for (ibrick = 0; ibrick < (*r); ibrick++) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], (*npname)[ibrick]); } for (ibrick = (*r); ibrick < (*p) + (*r); ibrick++) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], (*spname)[ibrick-(*r)]); } ibrick = (*p) + (*r); option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "Signal TMax"); ibrick++; option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "Signal SMax"); ibrick++; option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "Signal % SMax"); ibrick++; option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "Signal Area"); ibrick++; option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "Signal % Area"); ibrick++; option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "Sigma Resid"); ibrick++; option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "R^2"); ibrick++; option_data->brick_type[ibrick] = FUNC_FT_TYPE; option_data->brick_coef[ibrick] = ibrick; strcpy (option_data->brick_label[ibrick], "F-stat Regression"); } nopt++; continue; } /*----- -brick m type k label -----*/ if (strcmp(argv[nopt], "-brick") == 0) { nopt++; if (nopt+2 >= argc) NLfit_error ("need more arguments after -brick "); sscanf (argv[nopt], "%d", &ibrick); if ((ibrick < 0) || (ibrick >= option_data->numbricks)) NLfit_error ("illegal argument after -brick "); nopt++; if (strcmp(argv[nopt], "scoef") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; nopt++; sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival > (*p))) NLfit_error ("illegal argument after scoef "); option_data->brick_coef[ibrick] = ival + (*r); nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "ncoef") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; nopt++; sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival > (*r))) NLfit_error ("illegal argument after ncoef "); option_data->brick_coef[ibrick] = ival; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "tmax") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p); nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "smax") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 1; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "psmax") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 2; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "area") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 3; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "parea") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 4; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "resid") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 5; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "rsqr") == 0) { option_data->brick_type[ibrick] = FUNC_FIM_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 6; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "fstat") == 0) { option_data->brick_type[ibrick] = FUNC_FT_TYPE; option_data->brick_coef[ibrick] = (*r) + (*p) + 7; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); } else if (strcmp(argv[nopt], "tscoef") == 0) { option_data->brick_type[ibrick] = FUNC_TT_TYPE; nopt++; sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival > (*p))) NLfit_error ("illegal argument after tscoef "); option_data->brick_coef[ibrick] = ival + (*r); nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); calc_tstats = 1; } else if (strcmp(argv[nopt], "tncoef") == 0) { option_data->brick_type[ibrick] = FUNC_TT_TYPE; nopt++; sscanf (argv[nopt], "%d", &ival); if ((ival < 0) || (ival > (*r))) NLfit_error ("illegal argument after tncoef "); option_data->brick_coef[ibrick] = ival; nopt++; if (nopt >= argc) NLfit_error ("need more arguments after -brick "); strcpy (option_data->brick_label[ibrick], argv[nopt]); calc_tstats = 1; } else NLfit_error ("unable to interpret options after -brick "); nopt++; continue; } /*----- -sfit filename -----*/ if (strcmp(argv[nopt], "-sfit") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -sfit "); *sfit_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*sfit_filename); strcpy (*sfit_filename, argv[nopt]); nopt++; continue; } /*----- -snfit filename -----*/ if (strcmp(argv[nopt], "-snfit") == 0) { nopt++; if (nopt >= argc) NLfit_error ("need argument after -snfit "); *snfit_filename = malloc (sizeof(char) * MAX_NAME_LENGTH); MTEST (*snfit_filename); strcpy (*snfit_filename, argv[nopt]); nopt++; continue; } /*----- -jobs J -----*/ if( strcmp(argv[nopt],"-jobs") == 0 ){ /* RWCox */ nopt++ ; if (nopt >= argc) NLfit_error ("need J parameter after -jobs "); #ifdef PROC_MAX proc_numjob = strtol(argv[nopt],NULL,10) ; if( proc_numjob < 1 ){ fprintf(stderr,"** setting number of processes to 1!\n") ; proc_numjob = 1 ; } else if( proc_numjob > PROC_MAX ){ fprintf(stderr,"** setting number of processes to %d!\n",PROC_MAX); proc_numjob = PROC_MAX ; } #else fprintf(stderr,"** -jobs not supported in this version\n") ; #endif nopt++; continue; } /*----- unknown command -----*/ sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]); NLfit_error (message); } /*----- discard the model array -----*/ DESTROY_MODEL_ARRAY (model_array) ; } /*---------------------------------------------------------------------------*/ /* Routine to check whether one output file already exists. */ void check_one_output_file ( THD_3dim_dataset * dset_time, /* input 3d+time data set */ char * filename /* name of output file */ ) { THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ int ierror; /* number of errors in editing data */ char message[MAX_NAME_LENGTH]; /* error message */ /*----- make an empty copy of input dataset -----*/ new_dset = EDIT_empty_copy( dset_time ) ; ierror = EDIT_dset_items( new_dset , ADN_prefix , filename , ADN_label1 , filename , ADN_self_name , filename , ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE , ADN_none ) ; if( ierror > 0 ) { fprintf(stderr, "*** %d errors in attempting to create output dataset!\n", ierror ) ; exit(1) ; } if( THD_is_file(new_dset->dblk->diskptr->header_name) ) { sprintf (message, "Output dataset file %s already exists", new_dset->dblk->diskptr->header_name ) ; NLfit_error (message); } /*----- deallocate memory -----*/ THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; } /*---------------------------------------------------------------------------*/ /* Routine to check whether output files already exist. */ void check_output_files ( char * freg_filename, /* file name for regression f-statistics */ char * frsqr_filename, /* file name for R^2 statistics */ char * fsmax_filename, /* file name for signal signed maximum */ char * ftmax_filename, /* file name for epoch of signed maximum */ char * fpmax_filename, /* file name for max. percentage change */ char * farea_filename, /* file name for area under the signal */ char * fparea_filename, /* file name for % area under the signal */ char ** fncoef_filename, /* file name for noise model parameters */ char ** fscoef_filename, /* file name for signal model parameters */ char ** tncoef_filename, /* file name for noise model t-statistics */ char ** tscoef_filename, /* file name for signal model t-statistics */ char * bucket_filename, /* file name for bucket dataset */ char * sfit_filename, /* file name for 3d+time fitted signal model */ char * snfit_filename, /* file name for 3d+time fitted signal+noise */ THD_3dim_dataset * dset_time /* input 3d+time data set */ ) { int ip; /* parameter index */ if (freg_filename != NULL) check_one_output_file (dset_time, freg_filename); if (frsqr_filename != NULL) check_one_output_file (dset_time, frsqr_filename); if (fsmax_filename != NULL) check_one_output_file (dset_time, fsmax_filename); if (ftmax_filename != NULL) check_one_output_file (dset_time, ftmax_filename); if (fpmax_filename != NULL) check_one_output_file (dset_time, fpmax_filename); if (farea_filename != NULL) check_one_output_file (dset_time, farea_filename); if (fparea_filename != NULL) check_one_output_file (dset_time, fparea_filename); for (ip = 0; ip < MAX_PARAMETERS; ip++) { if (fncoef_filename[ip] != NULL) check_one_output_file (dset_time, fncoef_filename[ip]); if (fscoef_filename[ip] != NULL) check_one_output_file (dset_time, fscoef_filename[ip]); if (tncoef_filename[ip] != NULL) check_one_output_file (dset_time, tncoef_filename[ip]); if (tscoef_filename[ip] != NULL) check_one_output_file (dset_time, tscoef_filename[ip]); } if (bucket_filename != NULL) check_one_output_file (dset_time, bucket_filename); if (sfit_filename != NULL) check_one_output_file (dset_time, sfit_filename); if (snfit_filename != NULL) check_one_output_file (dset_time, snfit_filename); } /*---------------------------------------------------------------------------*/ /* Routine to check for valid inputs. */ void check_for_valid_inputs ( int nxyz, /* number of voxels in image */ int r, /* number of parameters in the noise model */ int p, /* number of parameters in the signal model */ float * min_nconstr, /* minimum parameter constraints for noise model */ float * max_nconstr, /* maximum parameter constraints for noise model */ float * min_sconstr, /* minimum parameter constraints for signal model */ float * max_sconstr, /* maximum parameter constraints for signal model */ int nrand, /* number of random vectors to generate */ int nbest, /* number of random vectors to keep */ char * freg_filename, /* file name for regression f-statistics */ char * frsqr_filename, /* file name for R^2 statistics */ char * fsmax_filename, /* file name for signal signed maximum */ char * ftmax_filename, /* file name for epoch of signed maximum */ char * fpmax_filename, /* file name for max. percentage change */ char * farea_filename, /* file name for area under the signal */ char * fparea_filename, /* file name for % area under the signal */ char ** fncoef_filename, /* file name for noise model parameters */ char ** fscoef_filename, /* file name for signal model parameters */ char ** tncoef_filename, /* file name for noise model t-statistics */ char ** tscoef_filename, /* file name for signal model t-statistics */ char * bucket_filename, /* file name for bucket dataset */ char * sfit_filename, /* file name for 3d+time fitted signal model */ char * snfit_filename, /* file name for 3d+time fitted signal+noise */ THD_3dim_dataset * dset_time /* input 3d+time data set */ ) { int ip; /* parameter index */ /*----- check if mask is right size -----*/ if (mask_vol != NULL) if (mask_nvox != nxyz) NLfit_error ("mask and input dataset bricks don't match in size"); /*----- check for valid constraints -----*/ for (ip = 0; ip < r; ip++) if (min_nconstr[ip] > max_nconstr[ip]) NLfit_error ("Must have minimum constraints <= maximum constraints"); for (ip = 0; ip < p; ip++) if (min_sconstr[ip] > max_sconstr[ip]) NLfit_error("Must have mininum constraints <= maximum constraints"); /*----- must have nbest <= nrand -----*/ if (nrand < nbest) NLfit_error ("Must have nbest <= nrand"); /*----- check whether any of the output files already exist -----*/ if( THD_deathcon() ) check_output_files (freg_filename, frsqr_filename, fsmax_filename, ftmax_filename, fpmax_filename, farea_filename, fparea_filename, fncoef_filename, fscoef_filename, tncoef_filename, tscoef_filename, bucket_filename, sfit_filename, snfit_filename, dset_time); } /*---------------------------------------------------------------------------*/ /* Allocate volume memory and fill with zeros. */ void zero_fill_volume (float ** fvol, int nxyz) { int ixyz; if( proc_numjob == 1 ){ /* 1 process ==> allocate locally */ *fvol = (float *) malloc (sizeof(float) * nxyz); MTEST(*fvol); for (ixyz = 0; ixyz < nxyz; ixyz++) (*fvol)[ixyz] = 0.0; } #ifdef PROC_MAX else { /* multiple processes ==> prepare for shared memory */ /* by remembering what to do */ proc_shm_arnum++ ; proc_shm_arsiz = (int *) realloc( proc_shm_arsiz , sizeof(int) *proc_shm_arnum ) ; proc_shm_ar = (float ***) realloc( proc_shm_ar , sizeof(float **)*proc_shm_arnum ) ; proc_shm_arsiz[proc_shm_arnum-1] = nxyz ; proc_shm_ar[proc_shm_arnum-1] = fvol ; /* actual allocation and pointer assignment (to *fvol) will take place in function proc_finalize_shm_volumes() */ } #endif } /*---------------------------------------------------------------------------*/ #ifdef PROC_MAX /*** signal handler for fatal errors -- make sure shared memory is deleted ***/ #include void proc_sigfunc(int sig) { char * sname ; int ii ; static volatile int fff=0 ; if( fff ) _exit(1); else fff=1 ; switch(sig){ default: sname = "unknown" ; break ; case SIGHUP: sname = "SIGHUP" ; break ; case SIGTERM: sname = "SIGTERM" ; break ; case SIGILL: sname = "SIGILL" ; break ; case SIGKILL: sname = "SIGKILL" ; break ; case SIGPIPE: sname = "SIGPIPE" ; break ; case SIGSEGV: sname = "SIGSEGV" ; break ; case SIGBUS: sname = "SIGBUS" ; break ; case SIGINT: sname = "SIGINT" ; break ; case SIGFPE: sname = "SIGFPE" ; break ; } if( proc_shmid > 0 ){ shmctl( proc_shmid , IPC_RMID , NULL ) ; proc_shmid = 0 ; } fprintf(stderr,"Fatal Signal %d (%s) received in job #%d\n", sig,sname,proc_ind) ; exit(1) ; } void proc_atexit( void ) /*** similarly - atexit handler ***/ { if( proc_shmid > 0 ) shmctl( proc_shmid , IPC_RMID , NULL ) ; } /*---------------------------------------------------------------*/ /*** This function is called to allocate all output volumes at once in shared memory, and set their pointers ***/ /** 24 Oct 2006: use mmap() instead of shared memory, if possible **/ #include #if !defined(MAP_ANON) && defined(MAP_ANONYMOUS) # define MAP_ANON MAP_ANONYMOUS #endif void proc_finalize_shm_volumes(void) { char kstr[32] ; int ii ; long long psum , twogig = 2ll * 1024ll * 1024ll * 1024ll ; if( proc_shm_arnum == 0 ) return ; /* should never happen */ /** 21 Oct 2005: check in big arithmetic how much mem we need **/ psum = 0 ; for( ii=0 ; ii < proc_shm_arnum ; ii++ ) psum += (long long)proc_shm_arsiz[ii] ; psum *= sizeof(float) ; #ifdef MAP_ANON if( psum >= twogig && ( sizeof(void *) < 8 || sizeof(size_t) < 8 ) ) /* too much for 32-bit */ #else if( psum >= twogig ) /* too much for shmem */ #endif ERROR_exit( "Total shared memory needed = %lld >= %lld (2 GB)\n" "** SUGGESTION: Use 3dZcutup to slice dataset into smaller pieces\n" "** and then 3dZcat to glue results back together.\n" "** SUGGESTION: Run on a 64-bit computer system, instead of 32-bit.\n" , psum,twogig) ; else INFO_message("total shared memory needed = %lld bytes",psum) ; proc_shmsize = psum ; /* global variable */ #if 0 /* old code */ if( proc_shm_arnum == 0 ) return ; /* should never happen */ proc_shmsize = 0 ; /* add up sizes of */ for( ii=0 ; ii < proc_shm_arnum ; ii++ ) /* all arrays for */ proc_shmsize += proc_shm_arsiz[ii] ; /* shared memory */ proc_shmsize *= sizeof(float) ; /* convert to byte count */ #endif /* create shared memory segment */ #ifdef MAP_ANON /** 24 Oct 2005: use mmap() instead of shmem **/ proc_shmptr = mmap( (void *)0 , (size_t)psum , PROT_READ | PROT_WRITE , MAP_ANON | MAP_SHARED , -1 , 0 ) ; if( proc_shmptr == NULL || proc_shmptr == (void *)(-1) ){ perror("** FATAL ERROR: Can't create shared mmap() segment\n" "** Unix message") ; exit(1) ; } #else /** old code: shared memory segment **/ UNIQ_idcode_fill( kstr ) ; /* unique string "key" */ proc_shmid = shm_create( kstr , proc_shmsize ) ; /* thd_iochan.c */ if( proc_shmid < 0 ){ fprintf(stderr,"\n** Can't create shared memory of size %d!\n" "** Try re-running without -jobs option!\n" , proc_shmsize ) ; /** if failed, print out some advice on how to tune SHMMAX **/ { char *cmd=NULL ; #if defined(LINUX) cmd = "cat /proc/sys/kernel/shmmax" ; #elif defined(SOLARIS) cmd = "/usr/sbin/sysdef | grep SHMMAX" ; #elif defined(DARWIN) cmd = "sysctl -n kern.sysv.shmmax" ; #endif if( cmd != NULL ){ FILE *fp = popen( cmd , "r" ) ; if( fp != NULL ){ long long smax=0 ; fscanf(fp,"%lld",&smax) ; pclose(fp) ; if( smax > 0 ) fprintf(stderr , "\n" "** POSSIBLY USEFUL ADVICE:\n" "** Current max shared memory size = %u bytes.\n" "** For information on how to change this, see\n" "** http://afni.nimh.nih.gov/afni/parallize.htm\n" "** and also contact your system administrator.\n" , smax ) ; } } } exit(1) ; } #endif /* MAP_ANON */ /* set a signal handler to catch most fatal errors and delete the shared memory segment if program crashes */ signal(SIGPIPE,proc_sigfunc) ; signal(SIGSEGV,proc_sigfunc) ; signal(SIGINT ,proc_sigfunc) ; signal(SIGFPE ,proc_sigfunc) ; signal(SIGBUS ,proc_sigfunc) ; signal(SIGHUP ,proc_sigfunc) ; signal(SIGTERM,proc_sigfunc) ; signal(SIGILL ,proc_sigfunc) ; signal(SIGKILL,proc_sigfunc) ; signal(SIGPIPE,proc_sigfunc) ; atexit(proc_atexit) ; /* or when the program exits */ #ifdef MAP_ANON fprintf(stderr , "++ mmap() memory allocated: %lld bytes\n" , proc_shmsize ) ; #else fprintf(stderr , "++ Shared memory allocated: %lld bytes at id=%d\n" , proc_shmsize , proc_shmid ) ; #endif /* get pointer to shared memory segment we just created */ #ifndef MAP_ANON if(proc_shmid > 0) { proc_shmptr = shm_attach( proc_shmid ) ; /* thd_iochan.c */ if( proc_shmptr == NULL ){ fprintf(stderr,"\n** Can't attach to shared memory!?\n" "** This is bizarre.\n" ) ; shmctl( proc_shmid , IPC_RMID , NULL ) ; exit(1) ; } } #endif /* clear all shared memory to zero*/ memset( proc_shmptr , 0 , proc_shmsize ) ; /* fix the local pointers to arrays in shared memory */ *proc_shm_ar[0] = (float *) proc_shmptr ; for( ii=1 ; ii < proc_shm_arnum ; ii++ ) *proc_shm_ar[ii] = *proc_shm_ar[ii-1] + proc_shm_arsiz[ii] ; return; } /*-------------------------------------------------------------*/ /*** This function replaces free(); it won't try to free things stored in the shared memory ***/ void proc_free( void *ptr ) { int ii ; if( ptr == NULL ) return ; /* from if ( proc_shmid == 0 ) 25 Oct 2006 DRG */ if( proc_shmptr == NULL ){ free(ptr); return; } /* no shm */ for( ii=0 ; ii < proc_shm_arnum ; ii++ ) if( ((float *)ptr) == *proc_shm_ar[ii] ) return; free(ptr); return; } #undef free /* replace use of library free() */ #define free proc_free /* with proc_free() function */ #endif /* PROC_MAX */ /*---------------------------------------------------------------------------*/ /* Perform all program initialization. */ void initialize_program ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ int * ignore, /* delete this number of points from time series */ char ** nname, /* name of noise model */ char ** sname, /* name of signal model */ vfp * nmodel, /* pointer to noise model */ vfp * smodel, /* pointer to signal model */ int *r, /* number of parameters in the noise model */ int *p, /* number of parameters in the signal model */ char *** npname, /* noise parameter names */ char *** spname, /* signal parameter names */ float ** min_nconstr, /* minimum parameter constraints for noise model */ float ** max_nconstr, /* maximum parameter constraints for noise model */ float ** min_sconstr, /* minimum parameter constraints for signal model */ float ** max_sconstr, /* maximum parameter constraints for signal model */ int * nabs, /* use absolute constraints for noise parameters */ int * nrand, /* number of random vectors to generate */ int * nbest, /* number of random vectors to keep */ float * rms_min, /* minimum rms error to reject reduced model */ float * fdisp, /* minimum f-statistic for display */ int *progress, /* progress report every nth voxel */ char ** input_filename, /* file name of input 3d+time dataset */ char ** tfilename, /* file name for time point series */ char ** freg_filename, /* file name for regression f-statistics */ char ** frsqr_filename, /* file name for R^2 statistics */ char ** fsmax_filename, /* file name for signal signed maximum */ char ** ftmax_filename, /* file name for time of signed maximum */ char ** fpmax_filename, /* file name for max. percentage change */ char ** farea_filename, /* file name for area under the signal */ char ** fparea_filename, /* file name for % area under the signal */ char *** fncoef_filename, /* file name for noise model parameters */ char *** fscoef_filename, /* file name for signal model parameters */ char *** tncoef_filename, /* file name for noise model t-statistics */ char *** tscoef_filename, /* file name for signal model t-statistics */ char ** sfit_filename, /* file name for 3d+time fitted signal model */ char ** snfit_filename, /* file name for 3d+time fitted signal+noise */ THD_3dim_dataset ** dset_time, /* input 3d+time data set */ int * nxyz, /* number of voxels in image */ int * ts_length, /* length of time series data */ float *** x_array, /* independent variable matrix */ float ** ts_array, /* input time series array */ float ** par_rdcd, /* estimated parameters for the reduced model */ float ** par_full, /* estimated parameters for the full model */ float ** tpar_full, /* t-statistic of the parameters in full model */ float ** rmsreg_vol, /* root-mean-square error for the full model */ float ** freg_vol, /* f-statistic for the full regression model */ float ** rsqr_vol, /* R^2 volume data */ float ** smax_vol, /* volume of signed maximum of signal */ float ** tmax_vol, /* volume of epoch of signed maximum of signal */ float ** pmax_vol, /* max. percentage change due to signal */ float ** area_vol, /* area between signal and baseline */ float ** parea_vol, /* percentage area between signal and baseline */ float *** ncoef_vol, /* volume of noise model parameters */ float *** scoef_vol, /* volume of signal model parameters */ float *** tncoef_vol, /* volume of noise model t-statistics */ float *** tscoef_vol, /* volume of signal model t-statistics */ float *** sfit_vol, /* voxelwise 3d+time fitted signal model */ float *** snfit_vol, /* voxelwise 3d+time fitted signal+noise model */ NL_options * option_data /* bucket dataset options */ ) { int dimension; /* dimension of full model */ int ip; /* parameter index */ int it; /* time index */ MRI_IMAGE * im, * flim; /* pointers to image structures -- used to read 1D ASCII */ int nt; /* number of points in 1D x data file */ float * tar; int ixyz; /* voxel index */ /*----- save command line for history notes -----*/ commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ; /*----- get command line inputs -----*/ get_options(argc, argv, ignore, nname, sname, nmodel, smodel, r, p, npname, spname, min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs, nrand, nbest, rms_min, fdisp, progress, input_filename, tfilename, freg_filename, frsqr_filename, fsmax_filename, ftmax_filename, fpmax_filename, farea_filename, fparea_filename, fncoef_filename, fscoef_filename, tncoef_filename, tscoef_filename, sfit_filename, snfit_filename, dset_time, nxyz, ts_length, option_data); /*----- check for valid inputs -----*/ check_for_valid_inputs (*nxyz, *r, *p, *min_nconstr, *max_nconstr, *min_sconstr, *max_sconstr, *nrand, *nbest, *freg_filename, *frsqr_filename, *fsmax_filename, *ftmax_filename, *fpmax_filename, *farea_filename, *fparea_filename, *fncoef_filename, *fscoef_filename, *tncoef_filename, *tscoef_filename, *sfit_filename, *snfit_filename, option_data->bucket_filename, *dset_time); /*----- allocate space for input time series -----*/ *ts_length = *ts_length - *ignore; *ts_array = (float *) malloc (sizeof(float) * (*ts_length)); MTEST (*ts_array); /*----- allocate space for independent variable matrix -----*/ *x_array = (float **) malloc (sizeof(float *) * (*ts_length)); MTEST (*x_array); for (it = 0; it < (*ts_length); it++) { (*x_array)[it] = (float *) malloc (sizeof(float) * 3); MTEST ((*x_array)[it]); } /*----- initialize independent variable matrix -----*/ if (*tfilename == NULL) { if( inTR && dsTR > 0.0 ){ /* 22 July 1998 */ DELT = dsTR ; fprintf(stderr,"--- computing with TR = %g\n",DELT) ; } for (it = 0; it < (*ts_length); it++) { (*x_array)[it][0] = 1.0; (*x_array)[it][1] = it * DELT; (*x_array)[it][2] = (it * DELT) * (it * DELT); } } else { flim = mri_read_1D(*tfilename) ; if (flim == NULL) NLfit_error ("Unable to read time reference file \n"); nt = flim -> nx; if (nt < (*ts_length)) NLfit_error ("Time reference array is too short"); tar = MRI_FLOAT_PTR(flim) ; for (it = 0; it < (*ts_length); it++) { (*x_array)[it][0] = 1.0; (*x_array)[it][1] = tar[it] ; (*x_array)[it][2] = tar[it] * tar[it]; } mri_free (flim); } /*--- 24 Jul 2006: special change to x_array[][2] for Linear+Ort [RWCox] ---*/ if( strcmp(*nname,"Linear+Ort") == 0 ){ char *fname ; MRI_IMAGE *fim ; int nx ; float *far ; fname = my_getenv("AFNI_ORTMODEL_REF") ; if( fname == NULL ) ERROR_exit("Linear+Ort model: 'AFNI_ORTMODEL_REF' not set") ; fim = mri_read_1D(fname) ; if( fim == NULL || fim->nx < 2 ) ERROR_exit( "Linear+Ort model: can't read file AFNI_ORTMODEL_REF='%s'",fname) ; nx = fim->nx ; far = MRI_FLOAT_PTR(fim) ; if( nx != (*ts_length) || fim->ny > 1 ) WARNING_message( "Linear+Ort model: AFNI_ORTMODEL_REF='%s' has nx=%d ny=%d; data length=%d", fname , nx , fim->ny , *ts_length ) ; else INFO_message( "Linear+Ort model: AFNI_ORTMODEL_REF='%s' loaded; length=%d" , fname , nx ) ; for( it=0 ; it < (*ts_length); it++) (*x_array)[it][2] = (it < nx) ? far[it] : 0.0f ; } /*----- allocate memory space for parameters -----*/ dimension = (*r) + (*p); *par_rdcd = (float *) malloc (sizeof(float) * dimension); MTEST (*par_rdcd); *par_full = (float *) malloc (sizeof(float) * dimension); MTEST (*par_full); *tpar_full = (float *) malloc (sizeof(float) * dimension); MTEST (*tpar_full); /*----- allocate memory space for volume data -----*/ zero_fill_volume( freg_vol , *nxyz ) ; if ((*freg_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( rmsreg_vol , *nxyz ) ; if ((*frsqr_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( rsqr_vol , *nxyz ) ; if ((*fsmax_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( smax_vol , *nxyz ) ; if ((*ftmax_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( tmax_vol , *nxyz ) ; if ((*fpmax_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( pmax_vol , *nxyz ) ; if ((*farea_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( area_vol , *nxyz ) ; if ((*fparea_filename != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( parea_vol , *nxyz ) ; *ncoef_vol = (float **) malloc (sizeof(float *) * (*r)); MTEST (*ncoef_vol); *tncoef_vol = (float **) malloc (sizeof(float *) * (*r)); MTEST (*tncoef_vol); for (ip = 0; ip < (*r); ip++) { if (((*fncoef_filename)[ip] != NULL) || ((*tncoef_filename)[ip] != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( &((*ncoef_vol)[ip]) , *nxyz ) ; else (*ncoef_vol)[ip] = NULL; if (((*tncoef_filename)[ip] != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( &((*tncoef_vol)[ip]) , *nxyz ) ; else (*tncoef_vol)[ip] = NULL; } *scoef_vol = (float **) malloc (sizeof(float *) * (*p)); MTEST (*scoef_vol); *tscoef_vol = (float **) malloc (sizeof(float *) * (*p)); MTEST (*tscoef_vol); for (ip = 0; ip < (*p); ip++) { if (((*fscoef_filename)[ip] != NULL) || ((*tscoef_filename)[ip] != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( &((*scoef_vol)[ip]) , *nxyz ) ; else (*scoef_vol)[ip] = NULL; if (((*tscoef_filename)[ip] != NULL) || (option_data->bucket_filename != NULL)) zero_fill_volume( &((*tscoef_vol)[ip]) , *nxyz ) ; else (*tscoef_vol)[ip] = NULL; } /*----- Allocate memory space for 3d+time fitted signal model -----*/ if (*sfit_filename != NULL) { *sfit_vol = (float **) malloc (sizeof(float *) * (*ts_length)); MTEST(*sfit_vol); for (it = 0; it < *ts_length; it++) zero_fill_volume( &((*sfit_vol)[it]) , *nxyz ) ; } /*----- Allocate memory space for 3d+time fitted signal+noise model -----*/ if (*snfit_filename != NULL) { *snfit_vol = (float **) malloc (sizeof(float *) * (*ts_length)); MTEST(*snfit_vol); for (it = 0; it < *ts_length; it++) zero_fill_volume( &((*snfit_vol)[it]) , *nxyz ) ; } #ifdef PROC_MAX if( proc_numjob > 1 ) proc_finalize_shm_volumes() ; /* RWCox */ #endif } /*---------------------------------------------------------------------------*/ /* Get the time series for one voxel from the AFNI 3d+time data set. */ void read_ts_array ( THD_3dim_dataset * dset_time, /* input 3d+time data set */ int iv, /* get time series for this voxel */ int ts_length, /* length of time series array */ int ignore, /* delete this number of points from time series */ float * ts_array /* input time series data for voxel #iv */ ) { MRI_IMAGE * im; /* intermediate float data */ float * ar; /* pointer to float data */ int it; /* time index */ /*----- Extract time series from 3d+time data set into MRI_IMAGE -----*/ im = THD_extract_series (iv, dset_time, 0); /*----- Verify extraction -----*/ if (im == NULL) NLfit_error ("Unable to extract data from 3d+time dataset"); /*----- Now extract time series from MRI_IMAGE -----*/ ar = MRI_FLOAT_PTR (im); for (it = 0; it < ts_length; it++) { ts_array[it] = ar[it+ignore]; } /*----- Release memory -----*/ mri_free (im); im = NULL; } /*---------------------------------------------------------------------------*/ /* Print out the given time series. */ void write_ts_array ( int ts_length, /* length of time series array */ float * ts_array /* time series data to be printed */ ) { int it; /* time index */ for (it = 0; it < ts_length; it++) printf ("%d %f \n", it, ts_array[it]); } /*---------------------------------------------------------------------------*/ /* Save results for output later. */ void save_results ( int iv, /* current voxel number */ vfp nmodel, /* pointer to noise model */ vfp smodel, /* pointer to signal model */ int r, /* number of parameters in the noise model */ int p, /* number of parameters in the signal model */ int novar, /* flag for insufficient variation in the data */ int ts_length, /* length of time series data */ float ** x_array, /* independent variable matrix */ float * par_full, /* estimated parameters for the full model */ float * tpar_full, /* t-statistic of the parameters in full model */ float rmsreg, /* root-mean-square error for the full model */ float freg, /* f-statistic for the full regression model */ float rsqr, /* R^2 (coef. of multiple determination) */ float smax, /* signed maximum of signal */ float tmax, /* epoch of signed maximum of signal */ float pmax, /* percentage change due to signal */ float area, /* area between signal and baseline */ float parea, /* percentage area between signal and baseline */ float * rmsreg_vol, /* root-mean-square for the full regression model */ float * freg_vol, /* f-statistic for the full regression model */ float * rsqr_vol, /* R^2 volume data */ float * smax_vol, /* signed maximum of signal volume data */ float * tmax_vol, /* epoch of signed maximum volume data */ float * pmax_vol, /* percentage change in signal volume data */ float * area_vol, /* area underneath signal volume data */ float * parea_vol, /* percentage area underneath signal volume data */ float ** ncoef_vol, /* volume of noise model parameters */ float ** scoef_vol, /* volume of signal model parameters */ float ** tncoef_vol, /* volume of noise model t-statistics */ float ** tscoef_vol, /* volume of signal model t-statistics */ float ** sfit_vol, /* voxelwise 3d+time fitted signal model */ float ** snfit_vol /* voxelwise 3d+time fitted signal+noise model */ ) { int ip; /* parameter index */ int it; /* time index */ float * s_array; /* fitted signal model time series */ float * n_array; /* fitted noise model time series */ /*----- save regression results into volume data -----*/ if (freg_vol != NULL) freg_vol[iv] = freg; if (rmsreg_vol != NULL) rmsreg_vol[iv] = rmsreg; if (rsqr_vol != NULL) rsqr_vol[iv] = rsqr; /*----- save signed maximum and epoch of signed maximum of signal -----*/ if (smax_vol != NULL) smax_vol[iv] = smax; if (tmax_vol != NULL) tmax_vol[iv] = tmax; /*----- save percentage change and area beneath the signal -----*/ if (pmax_vol != NULL) pmax_vol[iv] = pmax; if (area_vol != NULL) area_vol[iv] = area; if (parea_vol != NULL) parea_vol[iv] = parea; /*----- save noise parameter estimates -----*/ for (ip = 0; ip < r; ip++) { if (ncoef_vol[ip] != NULL) ncoef_vol[ip][iv] = par_full[ip]; if (tncoef_vol[ip] != NULL) tncoef_vol[ip][iv] = tpar_full[ip]; } /*----- save signal parameter estimates -----*/ for (ip = 0; ip < p; ip++) { if (scoef_vol[ip] != NULL) scoef_vol[ip][iv] = par_full[ip+r]; if (tscoef_vol[ip] != NULL) tscoef_vol[ip][iv] = tpar_full[ip+r]; } /*----- save fitted signal model time series -----*/ if (sfit_vol != NULL) { if (novar) { for (it = 0; it < ts_length; it++) sfit_vol[it][iv] = 0.0; } else { s_array = (float *) malloc (sizeof(float) * (ts_length)); MTEST (s_array); smodel (par_full + r, ts_length, x_array, s_array); for (it = 0; it < ts_length; it++) sfit_vol[it][iv] = s_array[it]; free (s_array); s_array = NULL; } } /*----- save fitted signal+noise model time series -----*/ if (snfit_vol != NULL) { n_array = (float *) malloc (sizeof(float) * (ts_length)); MTEST (n_array); nmodel (par_full, ts_length, x_array, n_array); for (it = 0; it < ts_length; it++) { snfit_vol[it][iv] = n_array[it]; } free (n_array); n_array = NULL; if (! novar) { s_array = (float *) malloc (sizeof(float) * (ts_length)); MTEST (s_array); smodel (par_full + r, ts_length, x_array, s_array); for (it = 0; it < ts_length; it++) { snfit_vol[it][iv] += s_array[it]; } free (s_array); s_array = NULL; } } } /*---------------------------------------------------------------------------*/ /* Convert one volume to another type, autoscaling: nxy = # voxels itype = input datum type ivol = pointer to input volume otype = output datum type ovol = pointer to output volume (again, must be pre-malloc-ed) Return value is the scaling factor used (0.0 --> no scaling). */ float EDIT_coerce_autoscale_new( int nxyz , int itype,void *ivol , int otype,void *ovol ) { float fac=0.0 , top ; if( MRI_IS_INT_TYPE(otype) ){ top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ; if (top == 0.0) fac = 0.0; else fac = MRI_TYPE_maxval[otype]/top; } EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ; return ( fac ); } /*---------------------------------------------------------------------------*/ /* Routine to write one AFNI data set. This data set may be either 'fitt' type (intensity + t-statistic) or 'fift' type (intensity + F-statistic). The intensity data is in ffim, and the corresponding statistic is in ftr. numdof = numerator degrees of freedom dendof = denominator degrees of freedom Note: if dendof = 0, then write 'fitt' type data set; if dendof > 0, then write 'fift' type data set. */ void write_afni_data (char * input_filename, int nxyz, char * filename, float * ffim, float * ftr, int numdof, int dendof) { int ii; /* voxel index */ THD_3dim_dataset * dset=NULL; /* input afni data set pointer */ THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ int iv; /* sub-brick index */ int ierror; /* number of errors in editing data */ int ibuf[32]; /* integer buffer */ float fbuf[MAX_STAT_AUX]; /* float buffer */ float fimfac; /* scale factor for short data */ /* int output_datum; */ /* data type for output data */ short * tsp; /* 2nd sub-brick data pointer */ void * vdif; /* 1st sub-brick data pointer */ int func_type; /* afni data set type */ float top, func_scale_short; /* parameters for scaling data */ char label[80]; /* label for output file history */ int nbad; /* number of bad voxels in volume */ /*----- read input dataset -----*/ dset = THD_open_dataset (input_filename) ; /* was THD_open_one...*/ if( ! ISVALID_3DIM_DATASET(dset) ){ fprintf(stderr,"*** Unable to open dataset file %s\n", input_filename); exit(1) ; } /*-- make an empty copy of this dataset, for eventual output --*/ new_dset = EDIT_empty_copy( dset ) ; /*----- Record history of dataset -----*/ tross_Copy_History( dset , new_dset ) ; sprintf (label, "Output prefix: %s", filename); if( commandline != NULL ) tross_multi_Append_History( new_dset , commandline,label,NULL ) ; else tross_Append_History ( new_dset, label); iv = DSET_PRINCIPAL_VALUE(dset) ; fprintf(stderr,"--output datum is %d\n", output_datum); /* output_datum = DSET_BRICK_TYPE(dset,iv) ; if( output_datum == MRI_byte ) output_datum = MRI_short ; */ ibuf[0] = output_datum ; ibuf[1] = MRI_short ; if (dendof == 0) func_type = FUNC_TT_TYPE; else func_type = FUNC_FT_TYPE; ierror = EDIT_dset_items( new_dset , ADN_prefix , filename , ADN_label1 , filename , ADN_self_name , filename , ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE , ADN_func_type , func_type , ADN_nvals , FUNC_nvals[func_type] , ADN_datum_array , ibuf , ADN_ntt , 0 , ADN_malloc_type, DATABLOCK_MEM_MALLOC , ADN_none ) ; if( ierror > 0 ){ fprintf(stderr, "*** %d errors in attempting to create output dataset!\n", ierror ) ; exit(1) ; } if( THD_is_file(new_dset->dblk->diskptr->header_name) ){ fprintf(stderr, "*** Output dataset file %s already exists--cannot continue!\a\n", new_dset->dblk->diskptr->header_name ) ; exit(1) ; } /*----- deleting exemplar dataset -----*/ THD_delete_3dim_dataset( dset , False ) ; dset = NULL ; nbad = thd_floatscan( nxyz , ffim ) ; if(nbad) fprintf(stderr,"++ %d bad floating point values in combined dataset\n",nbad); /*----- allocate memory for output data -----*/ vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz ); if (vdif == NULL) NLfit_error ("Unable to allocate memory for vdif"); tsp = (short *) malloc( sizeof(short) * nxyz ); if (tsp == NULL) NLfit_error ("Unable to allocate memory for tsp"); /*----- attach bricks to new data set -----*/ mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0)); mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1)); /*----- convert data type to output specification -----*/ fimfac = EDIT_coerce_autoscale_new (nxyz, MRI_float, ffim, output_datum, vdif); if (fimfac != 0.0) fimfac = 1.0 / fimfac; #define TOP_SS 32700 if (dendof == 0) /* t-statistic */ { top = TOP_SS/FUNC_TT_SCALE_SHORT; func_scale_short = FUNC_TT_SCALE_SHORT; } else /* F-statistic */ { top = TOP_SS/FUNC_FT_SCALE_SHORT; func_scale_short = FUNC_FT_SCALE_SHORT; } for (ii = 0; ii < nxyz; ii++) { if (ftr[ii] > top) tsp[ii] = TOP_SS; else if (ftr[ii] < -top) tsp[ii] = -TOP_SS; else if (ftr[ii] >= 0.0) tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5); else tsp[ii] = (short) (func_scale_short * ftr[ii] - 0.5); } /*----- write afni data set -----*/ printf("++ Writing combined dataset into %s\n",DSET_BRIKNAME(new_dset)) ; fbuf[0] = numdof; fbuf[1] = dendof; for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ; (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ; fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ; fbuf[1] = 1.0 / func_scale_short ; (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ; THD_load_statistics( new_dset ) ; THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ; /*----- deallocate memory -----*/ THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; } /*---------------------------------------------------------------------------*/ /* Routine to write one bucket data set. */ void write_bucket_data ( int q, /* number of parameters in the noise model */ int p, /* number of parameters in the signal model */ int numdof, /* numerator dof for F-statistic */ int dendof, /* denominator dof for F-statistic */ int nxyz, /* number of voxels in image */ int n, /* length of time series data */ float * rmsreg_vol, /* root-mean-square error for the full model */ float * freg_vol, /* f-statistic for the full regression model */ float * rsqr_vol, /* R^2 volume data */ float * smax_vol, /* signed maximum of signal volume data */ float * tmax_vol, /* epoch of signed maximum volume data */ float * pmax_vol, /* percentage change in signal volume data */ float * area_vol, /* area underneath signal volume data */ float * parea_vol, /* percentage area underneath signal volume data */ float ** ncoef_vol, /* volume of noise model parameters */ float ** scoef_vol, /* volume of signal model parameters */ float ** tncoef_vol, /* volume of noise model t-statistics */ float ** tscoef_vol, /* volume of signal model t-statistics */ char * input_filename, NL_options * option_data /* user input options */ ) { const float EPSILON = 1.0e-10; THD_3dim_dataset * old_dset = NULL; /* prototype dataset */ THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */ char * output_prefix; /* prefix name for bucket dataset */ char * output_session; /* directory for bucket dataset */ int nbricks, ib; /* number of sub-bricks in bucket dataset */ short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */ float ** far = NULL; /* far[ib] points to data for sub-brick #ib */ float factor; /* factor is new scale factor for sub-brick #ib */ int brick_type; /* indicates statistical type of sub-brick */ int brick_coef; /* regression coefficient index for sub-brick */ char * brick_label; /* character string label for sub-brick */ int ierror; /* number of errors in editing data */ float * volume; /* volume of floating point data */ int dimension; /* dimension of full model = p + q */ char label[80]; /* label for output file history */ void * imptr; /* pointer to volume in correct datum type to actually write out*/ int nbad; /* number of bad floating point values in volume */ /*----- initialize local variables -----*/ nbricks = option_data->numbricks; output_prefix = option_data->bucket_filename; output_session = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH); strcpy (output_session, "./"); dimension = p + q; /*----- allocate memory -----*/ if(output_datum==MRI_short) { bar = (short **) malloc (sizeof(short *) * nbricks); MTEST (bar); } else { far = (float **) malloc (sizeof(float *) * nbricks); MTEST (far); } /*----- read first dataset -----*/ old_dset = THD_open_dataset (input_filename); /* was THD_open_one...*/ /*-- make an empty copy of this dataset, for eventual output --*/ new_dset = EDIT_empty_copy (old_dset); /*----- Record history of dataset -----*/ tross_Copy_History( old_dset , new_dset ) ; if( commandline != NULL ) tross_Append_History( new_dset , commandline ) ; sprintf (label, "Output prefix: %s", output_prefix); tross_Append_History ( new_dset, label); /*----- Modify some structural properties. Note that the nbricks just make empty sub-bricks, without any data attached. -----*/ ierror = EDIT_dset_items (new_dset, ADN_prefix, output_prefix, ADN_directory_name, output_session, ADN_type, HEAD_FUNC_TYPE, ADN_func_type, FUNC_BUCK_TYPE, ADN_ntt, 0, /* no time */ ADN_nvals, nbricks, ADN_malloc_type, DATABLOCK_MEM_MALLOC , ADN_none ) ; if( ierror > 0 ) { fprintf(stderr, "*** %d errors in attempting to create output dataset!\n", ierror); exit(1); } if (THD_is_file(DSET_HEADNAME(new_dset))) { fprintf(stderr, "*** Output dataset file %s already exists--cannot continue!\n", DSET_HEADNAME(new_dset)); exit(1); } /*----- deleting exemplar dataset -----*/ THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ; nbad = 0; /*----- loop over new sub-brick index, attach data array with EDIT_substitute_brick then put some strings into the labels and keywords, and modify the sub-brick scaling factor -----*/ for (ib = 0; ib < nbricks; ib++) { /*----- get information about this sub-brick -----*/ brick_type = option_data->brick_type[ib]; brick_coef = option_data->brick_coef[ib]; brick_label = option_data->brick_label[ib]; if (brick_type == FUNC_FIM_TYPE) { if (brick_coef < q) volume = ncoef_vol[brick_coef]; else if (brick_coef < p+q) volume = scoef_vol[brick_coef-q]; else if (brick_coef == p+q) volume = tmax_vol; else if (brick_coef == p+q+1) volume = smax_vol; else if (brick_coef == p+q+2) volume = pmax_vol; else if (brick_coef == p+q+3) volume = area_vol; else if (brick_coef == p+q+4) volume = parea_vol; else if (brick_coef == p+q+5) volume = rmsreg_vol; else if (brick_coef == p+q+6) volume = rsqr_vol; } else if (brick_type == FUNC_TT_TYPE) { if (brick_coef < q) volume = tncoef_vol[brick_coef]; else if (brick_coef < p+q) volume = tscoef_vol[brick_coef-q]; EDIT_BRICK_TO_FITT (new_dset, ib, dendof); } else if (brick_type == FUNC_FT_TYPE) { volume = freg_vol; EDIT_BRICK_TO_FIFT (new_dset, ib, numdof, dendof); } nbad += thd_floatscan( nxyz , volume ) ; /*----- allocate memory for output sub-brick -----*/ if(output_datum==MRI_short) { bar[ib] = (short *) malloc (sizeof(short) * nxyz); MTEST (bar[ib]); factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, MRI_short, bar[ib]); imptr = bar[ib]; } else { far[ib] = (float *) malloc (sizeof(float) * nxyz); MTEST (far[ib]); factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, output_datum, far[ib]); imptr = far[ib]; } if (factor < EPSILON) factor = 0.0; else factor = 1.0 / factor; /*----- edit the sub-brick -----*/ EDIT_BRICK_LABEL (new_dset, ib, brick_label); EDIT_BRICK_FACTOR (new_dset, ib, factor); /*----- attach image pointer to be sub-brick #ib -----*/ EDIT_substitute_brick (new_dset, ib, output_datum, imptr); } if(nbad) fprintf(stderr,"++ %d bad floating point values in dataset\n",nbad); /*----- write bucket data set -----*/ THD_load_statistics (new_dset); THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ; fprintf(stderr,"++ Wrote bucket dataset: %s\n",DSET_BRIKNAME(new_dset)) ; /*----- deallocate memory -----*/ /* if(output_datum==MRI_short) {*/ THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; /*} */ } /*---------------------------------------------------------------------------*/ /* Routine to write one AFNI 3d+time data set. */ void write_3dtime ( char * input_filename, /* file name of input 3d+time dataset */ int ts_length, /* length of time series data */ float ** vol_array, /* output time series volume data */ char * output_filename /* output afni data set file name */ ) { const float EPSILON = 1.0e-10; THD_3dim_dataset * dset = NULL; /* input afni data set pointer */ THD_3dim_dataset * new_dset = NULL; /* output afni data set pointer */ int ib; /* sub-brick index */ int ierror; /* number of errors in editing data */ int nxyz; /* total number of voxels */ float factor; /* factor is new scale factor for sub-brick #ib */ short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */ float ** far = NULL; /* far[ib] points to data for sub-brick #ib */ float * fbuf; /* float buffer */ float * volume; /* pointer to volume of data */ char label[80]; /* label for output file history */ int nbad; /* number of voxels in volume with bad floating point values */ /*----- Initialize local variables -----*/ dset = THD_open_dataset (input_filename); /* was THD_open_one...*/ nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; /*----- allocate memory -----*/ if(output_datum==MRI_short) { bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar); } else { far = (float **) malloc (sizeof(float *) * ts_length); MTEST (far); } fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf); for (ib = 0; ib < ts_length; ib++) fbuf[ib] = 0.0; /*-- make an empty copy of the prototype dataset, for eventual output --*/ new_dset = EDIT_empty_copy (dset); /*----- Record history of dataset -----*/ tross_Copy_History( dset , new_dset ) ; if( commandline != NULL ) tross_Append_History( new_dset , commandline ) ; sprintf (label, "Output prefix: %s", output_filename); tross_Append_History ( new_dset, label); /*----- delete prototype dataset -----*/ THD_delete_3dim_dataset (dset, False); dset = NULL ; ierror = EDIT_dset_items (new_dset, ADN_prefix, output_filename, ADN_label1, output_filename, ADN_self_name, output_filename, ADN_malloc_type, DATABLOCK_MEM_MALLOC, ADN_datum_all, output_datum, ADN_nvals, ts_length, ADN_ntt, ts_length, ADN_none); if( ierror > 0 ){ fprintf(stderr, "*** %d errors in attempting to create output dataset!\n", ierror ) ; exit(1) ; } if( THD_is_file(new_dset->dblk->diskptr->header_name) ){ fprintf(stderr, "*** Output dataset file %s already exists--cannot continue!\a\n", new_dset->dblk->diskptr->header_name ) ; exit(1) ; } nbad = 0; /*----- attach bricks to new data set -----*/ for (ib = 0; ib < ts_length; ib++) { /*----- Set pointer to appropriate volume -----*/ volume = vol_array[ib]; nbad += thd_floatscan( nxyz , volume ) ; if(output_datum==MRI_short) { /*----- Allocate memory for output sub-brick -----*/ bar[ib] = (short *) malloc (sizeof(short) * nxyz); MTEST (bar[ib]); /*----- Convert data type to short for this sub-brick -----*/ factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, output_datum, bar[ib]); if (factor < EPSILON) factor = 0.0; else factor = 1.0 / factor; fbuf[ib] = factor; /*----- attach bar[ib] to be sub-brick #ib -----*/ mri_fix_data_pointer (bar[ib], DSET_BRICK(new_dset,ib)); } else { /*----- Allocate memory for output sub-brick -----*/ far[ib] = (float *) malloc (sizeof(float) * nxyz); MTEST (far[ib]); /*----- Convert data type to float for this sub-brick -----*/ factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, output_datum, far[ib]); if (factor < EPSILON) factor = 0.0; else factor = 1.0 / factor; fbuf[ib] = factor; /*----- attach far[ib] to be sub-brick #ib -----*/ mri_fix_data_pointer (far[ib], DSET_BRICK(new_dset,ib)); } } if(nbad) fprintf(stderr,"++ %d bad floating point values in dataset\n",nbad); /*----- write afni data set -----*/ (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ; THD_load_statistics (new_dset); THD_write_3dim_dataset (NULL, NULL, new_dset, True); fprintf (stderr,"++ Wrote 3D+time dataset %s\n",DSET_BRIKNAME(new_dset)) ; /*----- deallocate memory -----*/ /* if(output_datum==MRI_short) {*/ THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ; free (fbuf); fbuf = NULL; /* }*/ } /*---------------------------------------------------------------------------*/ /* Write out the user requested output files. */ void output_results ( int r, /* number of parameters in the noise model */ int p, /* number of parameters in the signal model */ float * min_nconstr, /* minimum parameter constraints for noise model */ float * max_nconstr, /* maximum parameter constraints for noise model */ float * min_sconstr, /* minimum parameter constraints for signal model */ float * max_sconstr, /* maximum parameter constraints for signal model */ int nxyz, /* number of voxels in image */ int ts_length, /* length of time series data */ float * rmsreg_vol, /* root-mean-square error for the full model */ float * freg_vol, /* f-statistic for the full regression model */ float * rsqr_vol, /* R^2 volume data */ float * smax_vol, /* signed maximum of signal volume data */ float * tmax_vol, /* epoch of signed maximum volume data */ float * pmax_vol, /* percentage change in signal volume data */ float * area_vol, /* area underneath signal volume data */ float * parea_vol, /* percentage area underneath signal volume data */ float ** ncoef_vol, /* volume of noise model parameters */ float ** scoef_vol, /* volume of signal model parameters */ float ** tncoef_vol, /* volume of noise model t-statistics */ float ** tscoef_vol, /* volume of signal model t-statistics */ float ** sfit_vol, /* voxelwise 3d+time fitted signal model */ float ** snfit_vol, /* voxelwise 3d+time fitted signal+noise model */ char * input_filename, /* file name of input 3d+time dataset */ char * freg_filename, /* file name for f-statistics */ char * frsqr_filename, /* file name for R^2 statistics */ char * fsmax_filename, /* file name for signal signed maximum */ char * ftmax_filename, /* file name for time of signed maximum */ char * fpmax_filename, /* file name for percentage signal change */ char * farea_filename, /* file name for area underneath signal */ char * fparea_filename, /* file name for % area underneath signal */ char ** fncoef_filename, /* file name for noise model parameters */ char ** fscoef_filename, /* file name for signal model parameters */ char ** tncoef_filename, /* file name for noise model t-statistics */ char ** tscoef_filename, /* file name for signal model t-statistics */ char * sfit_filename, /* file name for 3d+time fitted signal model */ char * snfit_filename, /* file name for 3d+time fitted signal+noise */ NL_options * option_data /* user input options */ ) { int ip; /* parameter index */ int dimension; /* dimension of full model = r + p */ int numdof, dendof; /* numerator and denominator degrees of freedom */ /*----- Initialization -----*/ dimension = r + p; numdof = p; dendof = ts_length - dimension; /*----- Adjust dof if constraints force a parameter to be a constant -----*/ for (ip = 0; ip < r; ip++) if (min_nconstr[ip] == max_nconstr[ip]) dendof++; for (ip = 0; ip < p; ip++) if (min_sconstr[ip] == max_sconstr[ip]) { numdof--; dendof++; } /*----- write the bucket dataset -----*/ if (option_data->numbricks > 0) write_bucket_data (r, p, numdof, dendof, nxyz, ts_length, rmsreg_vol, freg_vol, rsqr_vol, smax_vol, tmax_vol, pmax_vol, area_vol, parea_vol, ncoef_vol, scoef_vol, tncoef_vol, tscoef_vol, input_filename, option_data); /*----- write out the f-statistics for the regression -----*/ if (freg_filename != NULL) write_afni_data (input_filename, nxyz, freg_filename, rmsreg_vol, freg_vol, numdof, dendof); /*----- write out the R^2 (coefficient of multiple determination) -----*/ if (frsqr_filename != NULL) write_afni_data (input_filename, nxyz, frsqr_filename, rsqr_vol, freg_vol, numdof, dendof); /*----- write out the signed maximum of signal estimates -----*/ if (fsmax_filename != NULL) write_afni_data (input_filename, nxyz, fsmax_filename, smax_vol, freg_vol, numdof, dendof); /*----- write out the epoch of signed maximum of signal estimates -----*/ if (ftmax_filename != NULL) write_afni_data (input_filename, nxyz, ftmax_filename, tmax_vol, freg_vol, numdof, dendof); /*----- write out the max. percentage change due to signal -----*/ if (fpmax_filename != NULL) write_afni_data (input_filename, nxyz, fpmax_filename, pmax_vol, freg_vol, numdof, dendof); /*----- write out the area between the signal and baseline -----*/ if (farea_filename != NULL) write_afni_data (input_filename, nxyz, farea_filename, area_vol, freg_vol, numdof, dendof); /*----- write out the percentage area between the signal and baseline -----*/ if (fparea_filename != NULL) write_afni_data (input_filename, nxyz, fparea_filename, parea_vol, freg_vol, numdof, dendof); /*----- write noise model parameters -----*/ for (ip = 0; ip < r; ip++) { if (tncoef_filename[ip] != NULL) write_afni_data (input_filename, nxyz, tncoef_filename[ip], ncoef_vol[ip], tncoef_vol[ip], dendof, 0); if (fncoef_filename[ip] != NULL) write_afni_data (input_filename, nxyz, fncoef_filename[ip], ncoef_vol[ip], freg_vol, numdof, dendof); } /*----- write signal model parameters -----*/ for (ip = 0; ip < p; ip++) { if (tscoef_filename[ip] != NULL) write_afni_data (input_filename, nxyz, tscoef_filename[ip], scoef_vol[ip], tscoef_vol[ip], dendof, 0); if (fscoef_filename[ip] != NULL) write_afni_data (input_filename, nxyz, fscoef_filename[ip], scoef_vol[ip], freg_vol, numdof, dendof); } /*----- write the fitted 3d+time signal model -----*/ if (sfit_filename != NULL) { write_3dtime (input_filename, ts_length, sfit_vol, sfit_filename); } /*----- write the fitted 3d+time signal+noise model -----*/ if (snfit_filename != NULL) { write_3dtime (input_filename, ts_length, snfit_vol, snfit_filename); } } /*---------------------------------------------------------------------------*/ /* Release all allocated memory space. */ void terminate_program ( int r, /* number of parameters in the noise model */ int p, /* number of parameters in the signal model */ int ts_length, /* length of time series data */ float *** x_array, /* independent variable matrix */ float ** ts_array, /* input time series array */ char ** nname, /* noise model name */ char *** npname, /* noise parameter labels */ float ** par_rdcd, /* estimated parameters for the reduced model */ float ** min_nconstr, /* min parameter constraints for noise model */ float ** max_nconstr, /* max parameter constraints for noise model */ char ** sname, /* signal model name */ char *** spname, /* signal parameter labels */ float ** par_full, /* estimated parameters for the full model */ float ** tpar_full, /* t-statistic of parameters in full model */ float ** min_sconstr, /* min parameter constraints for signal model */ float ** max_sconstr, /* max parameter constraints for signal model */ float ** rmsreg_vol, /* rms error for the full regression model */ float ** freg_vol, /* f-statistic for the full regression model */ float ** rsqr_vol, /* R^2 volume data */ float ** smax_vol, /* signed max. of signal volume data */ float ** tmax_vol, /* epoch of signed max. volume data */ float ** pmax_vol, /* percentage change in signal volume data */ float ** area_vol, /* area underneath signal volume data */ float ** parea_vol, /* percent area underneath signal volume data */ float *** ncoef_vol, /* noise model parameters volume data */ float *** scoef_vol, /* signal model parameters volume data */ float *** tncoef_vol, /* noise model t-statistics volume data */ float *** tscoef_vol, /* signal model t-statistics volume data */ float *** sfit_vol, /* voxelwise 3d+time fitted signal model */ float *** snfit_vol, /* voxelwise 3d+time fitted signal+noise */ char ** input_filename, /* file name of input 3d+time dataset */ char ** freg_filename, /* file name for regression f-statistics */ char ** frsqr_filename, /* file name for R^2 statistics */ char ** fsmax_filename, /* file name for signal signed maximum */ char ** ftmax_filename, /* file name for epoch of signed maximum */ char ** fpmax_filename, /* file name for percentage signal change */ char ** farea_filename, /* file name for area underneath signal */ char ** fparea_filename, /* file name for % area underneath signal */ char *** fncoef_filename, /* file name for noise model parameters */ char *** fscoef_filename, /* file name for signal model parameters */ char *** tncoef_filename, /* file name for noise model t-statistics */ char *** tscoef_filename, /* file name for signal model t-statistics */ char ** sfit_filename, /* file name for 3d+time fitted signal model */ char ** snfit_filename /* file name for 3d+time fitted signal+noise */ ) { int ip; /* parameter index */ int it; /* time index */ /*----- deallocate space for model names and parameters labels -----*/ if (*nname != NULL) { free (*nname); *nname = NULL; } if (*sname != NULL) { free (*sname); *sname = NULL; } for (ip = 0; ip < MAX_PARAMETERS; ip++) { if ((*npname)[ip] != NULL) { free ((*npname)[ip]); (*npname)[ip] = NULL; } if ((*spname)[ip] != NULL) { free ((*spname)[ip]); (*spname)[ip] = NULL; } } if (*npname != NULL) { free (*npname); *npname = NULL; } if (*spname != NULL) { free (*spname); *spname = NULL; } /*----- deallocate memory for parameter constraints -----*/ if (*min_nconstr != NULL) { free (*min_nconstr); *min_nconstr = NULL; } if (*max_nconstr != NULL) { free (*max_nconstr); *max_nconstr = NULL; } if (*min_sconstr != NULL) { free (*min_sconstr); *min_sconstr = NULL; } if (*max_sconstr != NULL) { free (*max_sconstr); *max_sconstr = NULL; } /*----- deallocate memory space for filenames -----*/ if (*input_filename != NULL) { free (*input_filename); *input_filename = NULL; } if (*freg_filename != NULL) { free (*freg_filename); *freg_filename = NULL; } if (*frsqr_filename != NULL) { free (*frsqr_filename); *frsqr_filename = NULL; } if (*fsmax_filename != NULL) { free (*fsmax_filename); *fsmax_filename = NULL; } if (*ftmax_filename != NULL) { free (*ftmax_filename); *ftmax_filename = NULL; } if (*fpmax_filename != NULL) { free (*fpmax_filename); *fpmax_filename = NULL; } if (*farea_filename != NULL) { free (*farea_filename); *farea_filename = NULL; } if (*fparea_filename != NULL) { free (*fparea_filename); *fparea_filename = NULL; } for (ip = 0; ip < MAX_PARAMETERS; ip++) { if ((*fncoef_filename)[ip] != NULL) { free ((*fncoef_filename)[ip]); (*fncoef_filename)[ip] = NULL; } if ((*fscoef_filename)[ip] != NULL) { free ((*fscoef_filename)[ip]); (*fscoef_filename)[ip] = NULL; } if ((*tncoef_filename)[ip] != NULL) { free ((*tncoef_filename)[ip]); (*tncoef_filename)[ip] = NULL; } if ((*tscoef_filename)[ip] != NULL) { free ((*tscoef_filename)[ip]); (*tscoef_filename)[ip] = NULL; } } if (*fncoef_filename != NULL) { free (*fncoef_filename); *fncoef_filename = NULL; } if (*fscoef_filename != NULL) { free (*fscoef_filename); *fscoef_filename = NULL; } if (*tncoef_filename != NULL) { free (*tncoef_filename); *tncoef_filename = NULL; } if (*tscoef_filename != NULL) { free (*tscoef_filename); *tscoef_filename = NULL; } if (*sfit_filename != NULL) { free (*sfit_filename); *sfit_filename = NULL; } if (*snfit_filename != NULL) { free (*snfit_filename); *snfit_filename = NULL; } /*----- deallocate space for input time series -----*/ if (*ts_array != NULL) { free (*ts_array); *ts_array = NULL; } /*----- deallocate space for independent variable matrix -----*/ for (it = 0; it < ts_length; it++) if ((*x_array)[it] != NULL) { free ((*x_array)[it]); (*x_array)[it] = NULL; } if (*x_array != NULL) { free (*x_array); *x_array = NULL; } /*----- deallocate space for parameters -----*/ if (*par_rdcd != NULL) { free (*par_rdcd); *par_rdcd = NULL; } if (*par_full != NULL) { free (*par_full); *par_full = NULL; } if (*tpar_full != NULL) { free (*tpar_full); *tpar_full = NULL; } /*----- deallocate space for volume data -----*/ if (*rmsreg_vol != NULL) { free (*rmsreg_vol); *rmsreg_vol = NULL; } if (*freg_vol != NULL) { free (*freg_vol); *freg_vol = NULL; } if (*rsqr_vol != NULL) { free (*rsqr_vol); *rsqr_vol = NULL; } if (*smax_vol != NULL) { free (*smax_vol); *smax_vol = NULL; } if (*tmax_vol != NULL) { free (*tmax_vol); *tmax_vol = NULL; } if (*pmax_vol != NULL) { free (*pmax_vol); *pmax_vol = NULL; } if (*area_vol != NULL) { free (*area_vol); *area_vol = NULL; } if (*parea_vol != NULL) { free (*parea_vol); *parea_vol = NULL; } if (*ncoef_vol != NULL) { for (ip = 0; ip < r; ip++) if ((*ncoef_vol)[ip] != NULL) { free ((*ncoef_vol)[ip]); (*ncoef_vol)[ip] = NULL; } free (*ncoef_vol); *ncoef_vol = NULL; } if (*tncoef_vol != NULL) { for (ip = 0; ip < r; ip++) if ((*tncoef_vol)[ip] != NULL) { free ((*tncoef_vol)[ip]); (*tncoef_vol)[ip] = NULL; } free (*tncoef_vol); *tncoef_vol = NULL; } if (*scoef_vol != NULL) { for (ip = 0; ip < p; ip++) if ((*scoef_vol)[ip] != NULL) { free ((*scoef_vol)[ip]); (*scoef_vol)[ip] = NULL; } free (*scoef_vol); *scoef_vol = NULL; } if (*tscoef_vol != NULL) { for (ip = 0; ip < p; ip++) if ((*tscoef_vol)[ip] != NULL) { free ((*tscoef_vol)[ip]); (*tscoef_vol)[ip] = NULL; } free (*tscoef_vol); *tscoef_vol = NULL; } if (*sfit_vol != NULL) { for (it = 0; it < ts_length; it++) if ((*sfit_vol)[it] != NULL) { free ((*sfit_vol)[it]); (*sfit_vol)[it] = NULL; } free (*sfit_vol); *sfit_vol = NULL; } if (*snfit_vol != NULL) { for (it = 0; it < ts_length; it++) if ((*snfit_vol)[it] != NULL) { free ((*snfit_vol)[it]); (*snfit_vol)[it] = NULL; } free (*snfit_vol); *snfit_vol = NULL; } } /*---------------------------------------------------------------------------*/ int main ( int argc, /* number of input arguments */ char ** argv /* array of input arguments */ ) { /*----- declare input option variables -----*/ NL_options option_data; /* bucket dataset options */ int nabs; /* use absolute constraints for noise parameters */ int nrand; /* number of random vectors to generate */ int nbest; /* number of random vectors to keep */ float rms_min; /* minimum rms error to reject reduced model */ float fdisp; /* minimum f-statistic for display */ int progress; /* show progress report every n voxels */ /*----- declare time series variables -----*/ THD_3dim_dataset * dset_time = NULL; /* input 3d+time data set */ int ts_length; /* length of time series data */ int ignore; /* delete this number of points from time series */ float ** x_array = NULL; /* independent variable matrix */ float * ts_array = NULL; /* input time series array */ int nxyz; /* number of voxels in image */ int iv; /* voxel counter */ /*----- declare reduced (noise) model variables -----*/ char * nname = NULL; /* noise model name */ vfp nmodel; /* pointer to noise model */ int r; /* number of parameters in the noise model */ char ** npname = NULL; /* noise parameter labels */ float * par_rdcd = NULL; /* estimated parameters for the reduced model */ float sse_rdcd; /* error sum of squares for the reduced model */ float * min_nconstr = NULL; /* min parameter constraints for noise model */ float * max_nconstr = NULL; /* max parameter constraints for noise model */ /*----- declare full (signal+noise) model variables -----*/ char * sname = NULL; /* signal model name */ vfp smodel; /* pointer to signal model */ int p; /* number of parameters in the signal model */ char ** spname = NULL; /* signal parameter labels */ float * par_full = NULL; /* estimated parameters for the full model */ float sse_full; /* error sum of squares for the full model */ float * tpar_full = NULL; /* t-statistic of parameters in full model */ float freg; /* f-statistic for the full regression model */ float rmsreg; /* rms error for the full regression model */ float rsqr; /* R^2 (coef. of multiple determination) */ float smax; /* signed maximum of signal */ float tmax; /* epoch of signed maximum of signal */ float pmax; /* percentage change due to signal */ float area; /* area between signal and baseline */ float parea; /* percent area between signal and baseline */ float * min_sconstr = NULL; /* min parameter constraints for signal model */ float * max_sconstr = NULL; /* max parameter constraints for signal model */ /*----- declare output volume data -----*/ float * rmsreg_vol = NULL; /* rms error for the full regression model */ float * freg_vol = NULL; /* f-statistic for the full regression model */ float * rsqr_vol = NULL; /* R^2 volume data */ float * smax_vol = NULL; /* signed max. of signal volume data */ float * tmax_vol = NULL; /* epoch of signed max. volume data */ float * pmax_vol = NULL; /* max. percentage change due to signal */ float * area_vol = NULL; /* area between signal and baseline */ float * parea_vol = NULL; /* percent area between signal and baseline */ float ** ncoef_vol = NULL; /* noise model parameters volume data */ float ** scoef_vol = NULL; /* signal model parameters volume data */ float ** tncoef_vol = NULL; /* noise model t-statistics volume data */ float ** tscoef_vol = NULL; /* signal model t-statistics volume data */ float ** sfit_vol = NULL; /* voxelwise 3d+time fitted signal model */ float ** snfit_vol = NULL; /* voxelwise 3d+time fitted signal+noise */ /*----- declare file name variables -----*/ char * input_filename = NULL; /* file name of input 3d+time dataset */ char * tfilename = NULL; /* file name of time points */ char * freg_filename = NULL; /* file name for regression f-statistics */ char * frsqr_filename= NULL; /* file name for R^2 statistics */ char * fsmax_filename = NULL; /* file name for signal signed maximum */ char * ftmax_filename = NULL; /* file name for time of signed maximum */ char * fpmax_filename = NULL; /* file name for max. percentage change */ char * farea_filename = NULL; /* file name for area under the signal */ char * fparea_filename = NULL; /* file name for % area under the signal */ char ** fncoef_filename = NULL; /* file name for noise model parameters */ char ** fscoef_filename = NULL; /* file name for signal model parameters */ char ** tncoef_filename = NULL; /* file name for noise model t-statistics */ char ** tscoef_filename = NULL; /* file name for signal model t-statistics */ char * sfit_filename = NULL; /* file name for fitted signal model */ char * snfit_filename = NULL; /* file name for fitted signal+noise model */ char * label; /* report results for one voxel */ int novar; /* flag for insufficient variation in the data */ int ixyz_bot , ixyz_top ; int voxel_count_index = 0 ; /* for g_voxel_count output 26 Oct 2006 */ /*----- start the elapsed time counter -----*/ (void) COX_clock_time() ; /*----- Identify software -----*/ #if 0 printf ("\n\n"); printf ("Program: %s \n", PROGRAM_NAME); printf ("Author: %s \n", PROGRAM_AUTHOR); printf ("Initial Release: %s \n", PROGRAM_INITIAL); printf ("Latest Revision: %s \n", PROGRAM_LATEST); printf ("\n"); #endif /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/ PRINT_VERSION("3dNLfim") ; AUTHOR(PROGRAM_AUTHOR); mainENTRY("3dNLfim main") ; machdep() ; { int new_argc ; char ** new_argv ; addto_args( argc , argv , &new_argc , &new_argv ) ; if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; } } /*----- program initialization -----*/ initialize_program (argc, argv, &ignore, &nname, &sname, &nmodel, &smodel, &r, &p, &npname, &spname, &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr, &nabs, &nrand, &nbest, &rms_min, &fdisp,&progress, &input_filename, &tfilename, &freg_filename, &frsqr_filename, &fsmax_filename, &ftmax_filename, &fpmax_filename, &farea_filename, &fparea_filename, &fncoef_filename, &fscoef_filename, &tncoef_filename, &tscoef_filename, &sfit_filename, &snfit_filename, &dset_time, &nxyz, &ts_length, &x_array, &ts_array, &par_rdcd, &par_full, &tpar_full, &rmsreg_vol, &freg_vol, &rsqr_vol, &smax_vol, &tmax_vol, &pmax_vol, &area_vol, &parea_vol, &ncoef_vol, &scoef_vol, &tncoef_vol, &tscoef_vol, &sfit_vol, &snfit_vol, &option_data); #if 0 #ifdef SAVE_RAN RAN_setup( nmodel , smodel , r , p , nabs , min_nconstr, max_nconstr, min_sconstr, max_sconstr, ts_length, x_array, nrand ) ; #endif #endif INFO_message( "At each voxel, will use %d best of %d random parameter sets", nbest , nrand ) ; ixyz_bot = 0 ; ixyz_top = nxyz ; /* RWCox */ #ifdef PROC_MAX if( proc_numjob > 1 ){ /*---- set up multiple processes ----*/ int vv , nvox=nxyz , nper , pp , nv ; pid_t newpid ; /* count number of voxels to compute with into nvox */ if( mask_vol != NULL ){ for( vv=nvox=0 ; vv < nxyz ; vv++ ) if( mask_vol[vv] != 0 ) nvox++ ; } if( nvox < proc_numjob ){ /* too few voxels for multiple jobs? */ proc_numjob = 1 ; } else { /* prepare jobs */ /* split voxels between jobs evenly */ #if 0 if( g_voxel_count ){ g_voxel_count = 0 ; WARNING_message("-voxel_count disabled by -jobs") ; } #endif nper = nvox / proc_numjob ; /* # voxels per job */ if( mask_vol == NULL ){ proc_vox_bot[0] = 0 ; for( pp=0 ; pp < proc_numjob ; pp++ ){ proc_vox_top[pp] = proc_vox_bot[pp] + nper ; if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ; } proc_vox_top[proc_numjob-1] = nxyz ; } else { proc_vox_bot[0] = 0 ; for( pp=0 ; pp < proc_numjob ; pp++ ){ for( nv=0,vv=proc_vox_bot[pp] ; /* count ahead until */ nv < nper && vv < nxyz ; vv++ ){ /* find nper voxels */ if( mask_vol[vv] != 0 ) nv++ ; /* inside the mask */ } proc_vox_top[pp] = vv ; if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ; } proc_vox_top[proc_numjob-1] = nxyz ; } /* make sure dataset is in memory before forks */ DSET_load(dset_time) ; /* so dataset will be common */ /* start processes */ fprintf(stderr,"++ Voxels in dataset: %d\n",nxyz) ; if( nvox < nxyz ) fprintf(stderr,"++ Voxels in mask: %d\n",nvox) ; fprintf(stderr,"++ Voxels per job: %d\n",nper) ; for( pp=1 ; pp < proc_numjob ; pp++ ){ ixyz_bot = proc_vox_bot[pp] ; /* these 3 variables */ ixyz_top = proc_vox_top[pp] ; /* are for the process */ proc_ind = pp ; /* we're about to fork */ newpid = fork() ; if( newpid == -1 ){ fprintf(stderr,"** Can't fork job #%d! Error exit!\n",pp); exit(1) ; } if( newpid == 0 ) break ; /* I'm the child */ proc_pid[pp] = newpid ; /* I'm the parent */ iochan_sleep(10) ; } if( pp == proc_numjob ){ /* only in the parent */ ixyz_bot = proc_vox_bot[0] ; /* set the 3 control */ ixyz_top = proc_vox_top[0] ; /* variables needed */ proc_ind = 0 ; /* below */ } fprintf(stderr,"++ Job #%d: processing voxels %d to %d; elapsed time=%.3f\n", proc_ind,ixyz_bot,ixyz_top-1,COX_clock_time()) ; } } #endif /* PROC_MAX */ if( proc_numjob == 1 ) fprintf(stderr,"++ Calculations starting; elapsed time=%.3f\n", COX_clock_time()) ; /*----- loop over voxels in the data set -----*/ for (iv = ixyz_bot; iv < ixyz_top; iv++) { /*----- check for mask -----*/ if (mask_vol != NULL) if (mask_vol[iv] == 0) continue; /*----- display progress for user (1-based) -----*/ if (g_voxel_count && proc_ind == 0 ) { /* only print every 10th 26 Oct 2006 [rickr] */ if( voxel_count_index % 10 == 0 ) fprintf(stderr,"\r++ voxel count: %8d (of %d)", iv+1, ixyz_top); voxel_count_index++ ; } AFNI_store_dset_index( iv, -1); /* set voxel index 03 Nov 2006 */ /*----- read the time series for voxel iv -----*/ read_ts_array (dset_time, iv, ts_length, ignore, ts_array); /*----- calculate the reduced (noise) model -----*/ calc_reduced_model (ts_length, r, x_array, ts_array, par_rdcd, &sse_rdcd); /*----- calculate the full (signal+noise) model -----*/ calc_full_model (nmodel, smodel, r, p, min_nconstr, max_nconstr, min_sconstr, max_sconstr, ts_length, x_array, ts_array, par_rdcd, sse_rdcd, nabs, nrand, nbest, rms_min, par_full, &sse_full, &novar); /*----- calculate statistics for the full model -----*/ analyze_results (nmodel, smodel, r, p, novar, min_nconstr, max_nconstr, min_sconstr, max_sconstr, ts_length, x_array, par_rdcd, sse_rdcd, par_full, sse_full, &rmsreg, &freg, &rsqr, &smax, &tmax, &pmax, &area, &parea, tpar_full); /*----- report results for this voxel -----*/ if ((freg >= fdisp && proc_ind == 0 ) && (iv % progress == 0)) { report_results (nname, sname, r, p, npname, spname, ts_length, par_rdcd, sse_rdcd, par_full, sse_full, tpar_full, rmsreg, freg, rsqr, smax, tmax, pmax, area, parea, &label); printf ("\n\nVoxel #%d\n", iv); printf ("%s \n", label); } /*----- save results for this voxel into volume data -----*/ save_results (iv, nmodel, smodel, r, p, novar, ts_length, x_array, par_full, tpar_full, rmsreg, freg, rsqr, smax, tmax, pmax, area, parea, rmsreg_vol, freg_vol, rsqr_vol, smax_vol, tmax_vol, pmax_vol, area_vol, parea_vol,ncoef_vol, scoef_vol, tncoef_vol, tscoef_vol, sfit_vol, snfit_vol); } if(g_voxel_count) fputc('\n',stderr); /*-- if this is a child process, we're done. if this is the parent process, wait for the children --*/ #ifdef PROC_MAX if( proc_numjob > 1 ){ if( proc_ind > 0 ){ /* death of child */ fprintf(stderr,"++ Job #%d finished; elapsed time=%.3f\n", proc_ind,COX_clock_time()) ; _exit(0) ; } else { /* parent waits for children */ int pp ; fprintf(stderr,"++ Job #0 waiting for children to finish; elapsed time=%.3f\n",COX_clock_time()) ; for( pp=1 ; pp < proc_numjob ; pp++ ) waitpid( proc_pid[pp] , NULL , 0 ) ; fprintf(stderr,"++ Job #0 now finishing up; elapsed time=%.3f\n", COX_clock_time()) ; } /* when get to here, only parent process is left alive, and all the results are in the shared memory segment arrays */ } #endif if( proc_numjob == 1 ) fprintf(stderr,"++ Calculations finished; elapsed time=%.3f\n", COX_clock_time()) ; /*----- delete input dataset -----*/ THD_delete_3dim_dataset( dset_time , False ) ; dset_time = NULL ; /*----- write requested output files -----*/ output_results (r, p, min_nconstr, max_nconstr, min_sconstr, max_sconstr, nxyz, ts_length, rmsreg_vol, freg_vol, rsqr_vol, smax_vol, tmax_vol, pmax_vol, area_vol, parea_vol, ncoef_vol, scoef_vol, tncoef_vol, tscoef_vol, sfit_vol, snfit_vol, input_filename, freg_filename, frsqr_filename, fsmax_filename, ftmax_filename, fpmax_filename, farea_filename, fparea_filename, fncoef_filename, fscoef_filename, tncoef_filename, tscoef_filename, sfit_filename, snfit_filename, &option_data); /*----- end of program -----*/ terminate_program (r, p, ts_length, &x_array, &ts_array, &nname, &npname, &par_rdcd, &min_nconstr, &max_nconstr, &sname, &spname, &par_full, &tpar_full, &min_sconstr, &max_sconstr, &rmsreg_vol, &freg_vol, &rsqr_vol, &smax_vol, &tmax_vol, &pmax_vol, &area_vol, &parea_vol, &ncoef_vol, &scoef_vol, &tncoef_vol, &tscoef_vol, &sfit_vol, &snfit_vol, &input_filename, &freg_filename, &frsqr_filename, &fsmax_filename, &ftmax_filename, &fpmax_filename, &farea_filename, &fparea_filename, &fncoef_filename, &fscoef_filename, &tncoef_filename, &tscoef_filename, &sfit_filename, &snfit_filename); if( nwin_pow > 0 && (nwin_sim > 0 || nwin_stp > 0) ) INFO_message( "# best via POWELL=%d; via SIMPLEX=%d; SIMPLEX then POWELL=%d", nwin_pow , nwin_sim , nwin_stp ) ; INFO_message("Program finished; elapsed time=%.3f\n",COX_clock_time()) ; exit (0); }