/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 2000-2001, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ /*---------------------------------------------------------------------------*/ /* Program to perform wavelet analysis of an FMRI 3d+time dataset. File: 3dWavelets.c Author: B. Douglas Ward Date: 28 March 2000 Mod: Added call to AFNI_logger. Date: 15 August 2001 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME. Date: 02 December 2002 Mod: Added -datum option. Date: 07 February 2006 [rickr] */ /*---------------------------------------------------------------------------*/ /* Software identification. */ #define PROGRAM_NAME "3dWavelets" /* name of this program */ #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */ #define PROGRAM_INITIAL "28 March 2000" /* date of initial program release */ #define PROGRAM_LATEST "02 December 2002" /* date of last program revision */ /*---------------------------------------------------------------------------*/ /* Include header files and source code. */ #include "mrilib.h" #include "Wavelets.h" #include "Wavelets.c" /*---------------------------------------------------------------------------*/ /* Global constants. */ #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */ #define MAX_FILTERS 20 /* max. number of filter specifications */ /*---------------------------------------------------------------------------*/ /* Data structure for operator inputs. */ typedef struct WA_options { int NFirst; /* first image from input 3d+time dataset to use */ int NLast; /* last image from input 3d+time dataset to use */ int N; /* number of usable data points from input data */ int f; /* number of parameters removed by the filter */ int p; /* number of parameters in the full model */ int q; /* number of parameters in the baseline model */ float fdisp; /* minimum F-statistic for display */ char * input_filename; /* input 3d+time dataset */ char * mask_filename; /* input mask dataset */ char * input1D_filename; /* input fMRI measurement time series */ int wavelet_type; /* indicates type of wavelet */ int num_stop_filters; /* number of wavelet stop filters */ int stop_band[MAX_FILTERS]; /* wavelet filter stop band */ int stop_mintr[MAX_FILTERS]; /* min. time for stop band */ int stop_maxtr[MAX_FILTERS]; /* max. time for stop band */ float * stop_filter; /* select wavelet coefs. to stop */ int num_base_filters; /* number of wavelet base filters */ int base_band[MAX_FILTERS]; /* wavelet filter base band */ int base_mintr[MAX_FILTERS]; /* min. time for base band */ int base_maxtr[MAX_FILTERS]; /* max. time for base band */ float * base_filter; /* select wavelet coefs. for baseline */ int num_sgnl_filters; /* number of wavelet signal filters */ int sgnl_band[MAX_FILTERS]; /* wavelet filter signal band */ int sgnl_mintr[MAX_FILTERS]; /* min. time for signal band */ int sgnl_maxtr[MAX_FILTERS]; /* max. time for signal band */ float * sgnl_filter; /* select wavelet coefs. for signal */ char * bucket_filename; /* bucket dataset file name */ char * coefts_filename; /* forward wavelet transform coefficients output */ char * fitts_filename; /* full model time series 3d+time output */ char * sgnlts_filename; /* signal model time series 3d+time output */ char * errts_filename; /* residual error time series 3d+time output */ int fout; /* flag to output F-statistics */ int rout; /* flag to output R^2 statistics */ int cout; /* flag to output wavelet coefficients */ int vout; /* flag to output variance map */ int stat_first; /* flag to output full model stats first */ int datum; /* data type for output bucket and time_series */ } WA_options; /*---------------------------------------------------------------------------*/ /* Routine to display 3dWavelets help menu. */ void display_help_menu() { printf ( "Program to perform wavelet analysis of an FMRI 3d+time dataset. \n" " \n" "Usage: \n" "3dWavelets \n" "-type wname wname = name of wavelet to use for the analysis \n" " At present, there are only two choices for wname: \n" " Haar --> Haar wavelets \n" " Daub --> Daubechies wavelets \n" "-input fname fname = filename of 3d+time input dataset \n" "[-input1D dname] dname = filename of single (fMRI) .1D time series \n" "[-mask mname] mname = filename of 3d mask dataset \n" "[-nfirst fnum] fnum = number of first dataset image to use in \n" " the wavelet analysis. (default = 0) \n" "[-nlast lnum] lnum = number of last dataset image to use in \n" " the wavelet analysis. (default = last) \n" "[-fdisp fval] Write (to screen) results for those voxels \n" " whose F-statistic is >= fval \n" " \n" "Filter options: \n" " \n" "[-filt_stop band mintr maxtr] Specify wavelet coefs. to set to zero \n" "[-filt_base band mintr maxtr] Specify wavelet coefs. for baseline model\n" "[-filt_sgnl band mintr maxtr] Specify wavelet coefs. for signal model \n" " where band = frequency band \n" " mintr = min. value for time window (in TR) \n" " maxtr = max. value for time window (in TR) \n" " \n" "Output options: \n" " \n" "[-datum DTYPE] Coerce the output data to be stored as type DTYPE, \n" " which may be short or float. (default = short) \n" " \n" "[-coefts cprefix] cprefix = prefix of 3d+time output dataset which \n" " will contain the forward wavelet transform \n" " coefficients \n" " \n" "[-fitts fprefix] fprefix = prefix of 3d+time output dataset which \n" " will contain the full model time series fit \n" " to the input data \n" " \n" "[-sgnlts sprefix] sprefix = prefix of 3d+time output dataset which \n" " will contain the signal model time series fit \n" " to the input data \n" " \n" "[-errts eprefix] eprefix = prefix of 3d+time output dataset which \n" " will contain the residual error time series \n" " from the full model fit to the input data \n" " \n" "The following options control the contents of the bucket dataset: \n" " \n" "[-fout] Flag to output the F-statistics \n" "[-rout] Flag to output the R^2 statistics \n" "[-cout] Flag to output the full model wavelet coefficients \n" "[-vout] Flag to output the sample variance (MSE) map \n" " \n" "[-stat_first] Flag to specify that the full model statistics will \n" " appear prior to the wavelet coefficients in the \n" " bucket dataset output \n" " \n" "[-bucket bprefix] bprefix = prefix of AFNI 'bucket' dataset containing\n" " parameters of interest, such as the F-statistic \n" " for significance of the wavelet signal model. \n" ); exit(0); } /*---------------------------------------------------------------------------*/ /* Routine to initialize the input options. */ void initialize_options ( WA_options * option_data /* wavelet analysis program options */ ) { int is; /* filter index */ /*----- Initialize default values -----*/ option_data->NFirst = 0; option_data->NLast = 32767; option_data->N = 0; option_data->f = 0; option_data->p = 0; option_data->q = 0; option_data->fdisp = -1.0; option_data->wavelet_type = 0; /*----- Initialize filter specifications -----*/ option_data->num_stop_filters = 0; option_data->num_base_filters = 0; option_data->num_sgnl_filters = 0; for (is = 0; is < MAX_FILTERS; is++) { option_data->stop_band[is] = 0; option_data->stop_mintr[is] = 0; option_data->stop_maxtr[is] = 0; option_data->base_band[is] = 0; option_data->base_mintr[is] = 0; option_data->base_maxtr[is] = 0; option_data->sgnl_band[is] = 0; option_data->sgnl_mintr[is] = 0; option_data->sgnl_maxtr[is] = 0; } option_data->stop_filter = NULL; option_data->base_filter = NULL; option_data->sgnl_filter = NULL; /*----- Initialize output flags -----*/ option_data->fout = 0; option_data->rout = 0; option_data->cout = 0; option_data->vout = 0; option_data->stat_first = 0; option_data->datum = MRI_short; /*----- Initialize file name character strings -----*/ option_data->input_filename = NULL; option_data->mask_filename = NULL; option_data->input1D_filename = NULL; option_data->bucket_filename = NULL; option_data->coefts_filename = NULL; option_data->fitts_filename = NULL; option_data->sgnlts_filename = NULL; option_data->errts_filename = NULL; } /*---------------------------------------------------------------------------*/ /* Routine to get user specified input options. */ void get_options ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ WA_options * option_data /* wavelet program options */ ) { int nopt = 1; /* input option argument counter */ int ival; /* integer input */ float fval; /* float input */ char message[MAX_NAME_LENGTH]; /* error message */ int ifilter; /* filter index */ /*-- addto the arglist, if user wants to --*/ 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 ; } } /*----- does user request help menu? -----*/ if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu(); /*----- Add to program log -----*/ AFNI_logger (PROGRAM_NAME,argc,argv); /*----- initialize the input options -----*/ initialize_options (option_data); /*----- main loop over input options -----*/ while (nopt < argc ) { /*----- -type wname -----*/ if (strncmp(argv[nopt], "-type", 5) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -type "); for (ival = 0; ival < MAX_WAVELET_TYPE; ival++) if (strncmp(argv[nopt], WAVELET_TYPE_name[ival], 4) == 0) break; if (ival >= MAX_WAVELET_TYPE) { sprintf(message,"Unrecognized wavelet type: %s\n", argv[nopt]); WA_error (message); } option_data->wavelet_type = ival; nopt++; continue; } /*----- -input filename -----*/ if (strncmp(argv[nopt], "-input", 8) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -input "); option_data->input_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->input_filename); strcpy (option_data->input_filename, argv[nopt]); nopt++; continue; } /*----- -input1D filename -----*/ if (strncmp(argv[nopt], "-input1D", 8) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -input1D "); option_data->input1D_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->input1D_filename); strcpy (option_data->input1D_filename, argv[nopt]); nopt++; continue; } /*----- -mask filename -----*/ if (strncmp(argv[nopt], "-mask", 5) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -mask "); option_data->mask_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->mask_filename); strcpy (option_data->mask_filename, argv[nopt]); nopt++; continue; } /*----- -nfirst num -----*/ if (strncmp(argv[nopt], "-nfirst", 7) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -nfirst "); sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("illegal argument after -nfirst "); option_data->NFirst = ival; nopt++; continue; } /*----- -nlast num -----*/ if (strncmp(argv[nopt], "-nlast", 6) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -nlast "); sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("illegal argument after -nlast "); option_data->NLast = ival; nopt++; continue; } /*----- -fdisp fval -----*/ if (strncmp(argv[nopt], "-fdisp", 6) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -fdisp "); sscanf (argv[nopt], "%f", &fval); option_data->fdisp = fval; nopt++; continue; } /*----- -filt_stop band mintr maxtr -----*/ if (strncmp(argv[nopt], "-filt_stop", 10) == 0) { nopt++; if (nopt+2 >= argc) WA_error ("need 3 arguments after -filt_stop"); ifilter = option_data->num_stop_filters; sscanf (argv[nopt], "%d", &ival); if (ival < -1) WA_error ("-filt_stop band mintr maxtr Require: -1 <= band"); option_data->stop_band[ifilter] = ival; nopt++; sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("-filt_stop band mintr maxtr Require: 0 <= mintr"); option_data->stop_mintr[ifilter] = ival; nopt++; sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("-filt_stop band mintr maxtr Require: 0 <= maxtr"); option_data->stop_maxtr[ifilter] = ival; nopt++; option_data->num_stop_filters++; continue; } /*----- -filt_base band mintr maxtr -----*/ if (strncmp(argv[nopt], "-filt_base", 10) == 0) { nopt++; if (nopt+2 >= argc) WA_error ("need 3 arguments after -filt_base"); ifilter = option_data->num_base_filters; sscanf (argv[nopt], "%d", &ival); if (ival < -1) WA_error ("-filt_base band mintr maxtr Require: -1 <= band"); option_data->base_band[ifilter] = ival; nopt++; sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("-filt_base band mintr maxtr Require: 0 <= mintr"); option_data->base_mintr[ifilter] = ival; nopt++; sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("-filt_base band mintr maxtr Require: 0 <= maxtr"); option_data->base_maxtr[ifilter] = ival; nopt++; option_data->num_base_filters++; continue; } /*----- -filt_sgnl band mintr maxtr -----*/ if (strncmp(argv[nopt], "-filt_sgnl", 10) == 0) { nopt++; if (nopt+2 >= argc) WA_error ("need 3 arguments after -filt_sgnl"); ifilter = option_data->num_sgnl_filters; sscanf (argv[nopt], "%d", &ival); if (ival < -1) WA_error ("-filt_sgnl band mintr maxtr Require: -1 <= band"); option_data->sgnl_band[ifilter] = ival; nopt++; sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("-filt_sgnl band mintr maxtr Require: 0 <= mintr"); option_data->sgnl_mintr[ifilter] = ival; nopt++; sscanf (argv[nopt], "%d", &ival); if (ival < 0) WA_error ("-filt_sgnl band mintr maxtr Require: 0 <= maxtr"); option_data->sgnl_maxtr[ifilter] = ival; nopt++; option_data->num_sgnl_filters++; continue; } /*----- -fout -----*/ if (strncmp(argv[nopt], "-fout", 5) == 0) { option_data->fout = 1; nopt++; continue; } /*----- -rout -----*/ if (strncmp(argv[nopt], "-rout", 5) == 0) { option_data->rout = 1; nopt++; continue; } /*----- -cout -----*/ if (strncmp(argv[nopt], "-cout", 5) == 0) { option_data->cout = 1; nopt++; continue; } /*----- -vout -----*/ if (strncmp(argv[nopt], "-vout", 5) == 0) { option_data->vout = 1; nopt++; continue; } /*----- -stat_first -----*/ if (strncmp(argv[nopt], "-stat_first", 11) == 0) { option_data->stat_first = 1; nopt++; continue; } /*----- -datum DTYPE -----*/ if (strncmp(argv[nopt], "-datum", 6) == 0) { nopt++; if (nopt >= argc) WA_error ("need argument after -datum "); if (!strcmp(argv[nopt],"short")) option_data->datum = MRI_short; else if (!strcmp(argv[nopt],"float")) option_data->datum = MRI_float; else WA_error ("illegal argument after -datum "); nopt++; continue; } /*----- -bucket filename -----*/ if (strncmp(argv[nopt], "-bucket", 6) == 0) { nopt++; if (nopt >= argc) WA_error ("need file prefixname after -bucket "); option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->bucket_filename); strcpy (option_data->bucket_filename, argv[nopt]); nopt++; continue; } /*----- -coefts filename -----*/ if (strncmp(argv[nopt], "-coefts", 8) == 0) { nopt++; if (nopt >= argc) WA_error ("need file prefixname after -coefts "); option_data->coefts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->coefts_filename); strcpy (option_data->coefts_filename, argv[nopt]); nopt++; continue; } /*----- -fitts filename -----*/ if (strncmp(argv[nopt], "-fitts", 7) == 0) { nopt++; if (nopt >= argc) WA_error ("need file prefixname after -fitts "); option_data->fitts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->fitts_filename); strcpy (option_data->fitts_filename, argv[nopt]); nopt++; continue; } /*----- -sgnlts filename -----*/ if (strncmp(argv[nopt], "-sgnlts", 7) == 0) { nopt++; if (nopt >= argc) WA_error ("need file prefixname after -sgnlts "); option_data->sgnlts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->sgnlts_filename); strcpy (option_data->sgnlts_filename, argv[nopt]); nopt++; continue; } /*----- -errts filename -----*/ if (strncmp(argv[nopt], "-errts", 6) == 0) { nopt++; if (nopt >= argc) WA_error ("need file prefixname after -errts "); option_data->errts_filename = malloc (sizeof(char)*MAX_NAME_LENGTH); MTEST (option_data->errts_filename); strcpy (option_data->errts_filename, argv[nopt]); nopt++; continue; } /*----- unknown command -----*/ sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]); WA_error (message); } } /*---------------------------------------------------------------------------*/ /* Read time series from specified file name. This file name may have a column selector attached. */ float * read_time_series ( char * ts_filename, /* time series file name (plus column index) */ int * ts_length /* output value for time series length */ ) { char message[MAX_NAME_LENGTH]; /* error message */ MRI_IMAGE * flim; /* pointer to image structures -- used to read 1D ASCII */ float * far; /* pointer to MRI_IMAGE floating point data */ int nx; /* number of time points in time series */ int ny; /* number of columns in time series file */ int iy; /* time series file column index */ int ipt; /* time point index */ float * ts_data = NULL; /* input time series data */ /*----- Read the time series file -----*/ flim = mri_read_1D(ts_filename) ; if (flim == NULL) { sprintf (message, "Unable to read time series file: %s", ts_filename); WA_error (message); } /*----- Set pointer to data, and set dimensions -----*/ far = MRI_FLOAT_PTR(flim); nx = flim->nx; ny = flim->ny; iy = 0 ; if( ny > 1 ){ fprintf(stderr,"WARNING: time series %s has %d columns\n",ts_filename,ny); } /*----- Save the time series data -----*/ *ts_length = nx; ts_data = (float *) malloc (sizeof(float) * nx); MTEST (ts_data); for (ipt = 0; ipt < nx; ipt++) ts_data[ipt] = far[ipt + iy*nx]; mri_free (flim); flim = NULL; return (ts_data); } /*---------------------------------------------------------------------------*/ /* Read the input data files. */ void read_input_data ( WA_options * option_data, /* wavelet program options */ THD_3dim_dataset ** dset_time, /* input 3d+time data set */ THD_3dim_dataset ** mask_dset, /* input mask data set */ float ** fmri_data, /* input fMRI time series data */ int * fmri_length /* length of fMRI time series */ ) { char message[MAX_NAME_LENGTH]; /* error message */ /*----- Read the input fMRI measurement data -----*/ if (option_data->input1D_filename != NULL) { /*----- Read the input fMRI 1D time series -----*/ *fmri_data = read_time_series (option_data->input1D_filename, fmri_length); if (*fmri_data == NULL) { sprintf (message, "Unable to read time series file: %s", option_data->input1D_filename); WA_error (message); } *dset_time = NULL; } else if (option_data->input_filename != NULL) { /*----- Read the input 3d+time dataset -----*/ *dset_time = THD_open_one_dataset (option_data->input_filename); if (!ISVALID_3DIM_DATASET(*dset_time)) { sprintf (message, "Unable to open data file: %s", option_data->input_filename); WA_error (message); } DSET_load(*dset_time) ; CHECK_LOAD_ERROR(*dset_time) ; if (option_data->mask_filename != NULL) { /*----- Read the input mask dataset -----*/ *mask_dset = THD_open_dataset (option_data->mask_filename); if (!ISVALID_3DIM_DATASET(*mask_dset)) { sprintf (message, "Unable to open mask file: %s", option_data->mask_filename); WA_error (message); } DSET_load(*mask_dset) ; } } else WA_error ("Must specify input data"); } /*---------------------------------------------------------------------------*/ /* 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 */ ) { char message[MAX_NAME_LENGTH]; /* error message */ THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ int ierror; /* number of errors in editing data */ /*----- 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 ) { sprintf (message, "*** %d errors in attempting to create output dataset!\n", ierror); WA_error (message); } if( THD_is_file(new_dset->dblk->diskptr->header_name) ) { sprintf (message, "Output dataset file %s already exists " " -- cannot continue!\a\n", new_dset->dblk->diskptr->header_name); WA_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 ( WA_options * option_data, /* wavelet analysis program options */ THD_3dim_dataset * dset_time /* input 3d+time data set */ ) { if (option_data->bucket_filename != NULL) check_one_output_file (dset_time, option_data->bucket_filename); if (option_data->coefts_filename != NULL) check_one_output_file (dset_time, option_data->coefts_filename); if (option_data->fitts_filename != NULL) check_one_output_file (dset_time, option_data->fitts_filename); if (option_data->sgnlts_filename != NULL) check_one_output_file (dset_time, option_data->sgnlts_filename); if (option_data->errts_filename != NULL) check_one_output_file (dset_time, option_data->errts_filename); } /*---------------------------------------------------------------------------*/ /* Routine to check for valid inputs. */ void check_for_valid_inputs ( WA_options * option_data, /* wavelet analysis program options */ THD_3dim_dataset * dset_time, /* input 3d+time data set */ THD_3dim_dataset * mask_dset /* input mask data set */ ) { char message[MAX_NAME_LENGTH]; /* error message */ /*----- If mask is used, check for compatible dimensions -----*/ if (mask_dset != NULL) { if ( (DSET_NX(dset_time) != DSET_NX(mask_dset)) || (DSET_NY(dset_time) != DSET_NY(mask_dset)) || (DSET_NZ(dset_time) != DSET_NZ(mask_dset)) ) { sprintf (message, "%s and %s have incompatible dimensions", option_data->input_filename, option_data->mask_filename); WA_error (message); } if (DSET_NVALS(mask_dset) != 1 ) WA_error ("Must specify 1 sub-brick from mask dataset"); } /*----- Check whether any of the output files already exist -----*/ if( THD_deathcon() ) check_output_files (option_data, dset_time); } /*---------------------------------------------------------------------------*/ /* Routine to initialize the filters and control variables. */ void initialize_filters ( WA_options * option_data, /* wavelet analysis program options */ THD_3dim_dataset * dset_time, /* input 3d+time data set */ THD_3dim_dataset * mask_dset, /* input mask data set */ int fmri_length /* length of input fMRI time series */ ) { int nt; /* number of images in input 3d+time dataset */ int NFirst; /* first image from input 3d+time dataset to use */ int NLast; /* last image from input 3d+time dataset to use */ int N; /* number of usable time points */ int i; /* filter coefficient index */ int f, q, p; /* numbers of model parameters */ /*----- Initialize local variables -----*/ if (option_data->input1D_filename != NULL) nt = fmri_length; else nt = DSET_NUM_TIMES (dset_time); /*----- Determine the number of usable time points -----*/ NFirst = option_data->NFirst; NLast = option_data->NLast; if (NLast > nt-1) NLast = nt-1; N = NLast - NFirst + 1; N = powerof2(my_log2(N)); NLast = N + NFirst + 1; option_data->N = N; option_data->NLast = NLast; /*----- Initialize the filters -----*/ option_data->stop_filter = FWT_1d_stop_filter (option_data->num_stop_filters, option_data->stop_band, option_data->stop_mintr, option_data->stop_maxtr, option_data->NFirst, option_data->N); option_data->base_filter = FWT_1d_pass_filter (option_data->num_base_filters, option_data->base_band, option_data->base_mintr, option_data->base_maxtr, option_data->NFirst, option_data->N); option_data->sgnl_filter = FWT_1d_pass_filter (option_data->num_sgnl_filters, option_data->sgnl_band, option_data->sgnl_mintr, option_data->sgnl_maxtr, option_data->NFirst, option_data->N); /*----- Count numbers of model parameters -----*/ f = 0; for (i = 0; i < N; i++) if (option_data->stop_filter[i] == 0.0) { f++; option_data->base_filter[i] = 0.0; option_data->sgnl_filter[i] = 0.0; } option_data->f = f; q = 0; for (i = 0; i < N; i++) if (option_data->base_filter[i] == 1.0) { q++; option_data->sgnl_filter[i] = 1.0; } option_data->q = q; p = 0; for (i = 0; i < N; i++) if (option_data->sgnl_filter[i] == 1.0) { p++; } option_data->p = p; return; } /*---------------------------------------------------------------------------*/ /* Allocate memory for output volumes. */ void allocate_memory ( WA_options * option_data, /* wavelet analysis options */ THD_3dim_dataset * dset, /* input 3d+time data set */ float *** coef_vol, /* array of volumes of signal model parameters */ float ** mse_vol, /* volume of full model mean square error */ float ** ffull_vol, /* volume of full model F-statistics */ float ** rfull_vol, /* volume of full model R^2 stats. */ float *** coefts_vol, /* volumes for forward wavelet transform coefs. */ float *** fitts_vol, /* volumes for wavelet filtered time series data */ float *** sgnlts_vol, /* volumes for signal model fit to input data */ float *** errts_vol /* volumes for residual errors */ ) { int nxyz; /* total number of voxels */ int ixyz; /* voxel index */ int it; /* time point index */ int N; /* number of usable data points */ int ip; /* parameter index */ int p; /* number of parameters in the full model */ /*----- Initialize local variables -----*/ nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; N = option_data->N; p = option_data->p; /*----- Allocate memory space for volume data -----*/ *coef_vol = (float **) malloc (sizeof(float *) * p); MTEST(*coef_vol); for (ip = 0; ip < p; ip++) { (*coef_vol)[ip] = (float *) malloc (sizeof(float) * nxyz); MTEST((*coef_vol)[ip]); for (ixyz = 0; ixyz < nxyz; ixyz++) { (*coef_vol)[ip][ixyz] = 0.0; } } *mse_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*mse_vol); *ffull_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*ffull_vol); *rfull_vol = (float *) malloc (sizeof(float) * nxyz); MTEST(*rfull_vol); for (ixyz = 0; ixyz < nxyz; ixyz++) { (*mse_vol)[ixyz] = 0.0; (*ffull_vol)[ixyz] = 0.0; (*rfull_vol)[ixyz] = 0.0; } /*----- Allocate memory for forward wavelet transform coefficients -----*/ if (option_data->coefts_filename != NULL) { *coefts_vol = (float **) malloc (sizeof(float **) * N); MTEST (*coefts_vol); for (it = 0; it < N; it++) { (*coefts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); MTEST ((*coefts_vol)[it]); for (ixyz = 0; ixyz < nxyz; ixyz++) (*coefts_vol)[it][ixyz] = 0.0; } } /*----- Allocate memory for filtered time series -----*/ if (option_data->fitts_filename != NULL) { *fitts_vol = (float **) malloc (sizeof(float **) * N); MTEST (*fitts_vol); for (it = 0; it < N; it++) { (*fitts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); MTEST ((*fitts_vol)[it]); for (ixyz = 0; ixyz < nxyz; ixyz++) (*fitts_vol)[it][ixyz] = 0.0; } } /*----- Allocate memory for signal model time series -----*/ if (option_data->sgnlts_filename != NULL) { *sgnlts_vol = (float **) malloc (sizeof(float **) * N); MTEST (*sgnlts_vol); for (it = 0; it < N; it++) { (*sgnlts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); MTEST ((*sgnlts_vol)[it]); for (ixyz = 0; ixyz < nxyz; ixyz++) (*sgnlts_vol)[it][ixyz] = 0.0; } } /*----- Allocate memory for residual time series -----*/ if (option_data->errts_filename != NULL) { *errts_vol = (float **) malloc (sizeof(float **) * N); MTEST (*errts_vol); for (it = 0; it < N; it++) { (*errts_vol)[it] = (float *) malloc (sizeof(float *) * nxyz); MTEST ((*errts_vol)[it]); for (ixyz = 0; ixyz < nxyz; ixyz++) (*errts_vol)[it][ixyz] = 0.0; } } } /*---------------------------------------------------------------------------*/ /* Perform all program initialization. */ void initialize_program ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ WA_options ** option_data, /* wavelet analysis options */ THD_3dim_dataset ** dset_time, /* input 3d+time data set */ THD_3dim_dataset ** mask_dset, /* input mask data set */ float ** fmri_data, /* input fMRI time series data */ int * fmri_length, /* length of fMRI time series */ float *** coef_vol, /* array of volumes of signal model parameters */ float ** mse_vol, /* volume of full model mean square error */ float ** ffull_vol, /* volume of full model F-statistics */ float ** rfull_vol, /* volume of full model R^2 stats. */ float *** coefts_vol, /* volumes for forward wavelet transform coefs. */ float *** fitts_vol, /* volumes for wavelet filtered time series */ float *** sgnlts_vol, /* volumes for signal model fit to input data */ float *** errts_vol /* volumes for residual errors */ ) { /*----- Allocate memory -----*/ *option_data = (WA_options *) malloc (sizeof(WA_options)); /*----- Get command line inputs -----*/ get_options (argc, argv, *option_data); /*----- Read input data -----*/ read_input_data (*option_data, dset_time, mask_dset, fmri_data, fmri_length); /*----- Check for valid inputs -----*/ check_for_valid_inputs (*option_data, *dset_time, *mask_dset); /*----- Initialize the filters and control variables -----*/ initialize_filters (*option_data, *dset_time, *mask_dset, *fmri_length); /*----- Allocate memory for output volumes -----*/ if ((*option_data)->input1D_filename == NULL) allocate_memory (*option_data, *dset_time, coef_vol, mse_vol, ffull_vol, rfull_vol, coefts_vol, fitts_vol, sgnlts_vol, errts_vol); } /*---------------------------------------------------------------------------*/ /* Get the time series for one voxel from the AFNI 3d+time data set. */ void extract_ts_array ( THD_3dim_dataset * dset_time, /* input 3d+time dataset */ int iv, /* get time series for this voxel */ float * ts_array /* time series data for voxel #iv */ ) { MRI_IMAGE * im = NULL; /* intermediate float data */ float * ar = NULL; /* pointer to float data */ int ts_length; /* length of input 3d+time data set */ 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) WA_error ("Unable to extract data from 3d+time dataset"); /*----- Now extract time series from MRI_IMAGE -----*/ ts_length = DSET_NUM_TIMES (dset_time); ar = MRI_FLOAT_PTR (im); for (it = 0; it < ts_length; it++) { ts_array[it] = ar[it]; } /*----- Release memory -----*/ mri_free (im); im = NULL; } /*---------------------------------------------------------------------------*/ /* Save results for this voxel. */ void save_voxel ( WA_options * option_data, /* wavelet analysis options */ int iv, /* current voxel index */ float * coef, /* regression parameters */ float mse, /* full model mean square error */ float ffull, /* full model F-statistic */ float rfull, /* full model R^2 stat. */ float * coefts, /* forward wavelet transform coefficients */ float * fitts, /* wavelet filtered time series */ float * sgnlts, /* signal model fitted time series */ float * errts, /* signal model residual error time series */ float ** coef_vol, /* array of volumes of signal model parameters */ float * mse_vol, /* volume of full model mean square error */ float * ffull_vol, /* volume of full model F-statistics */ float * rfull_vol, /* volume of full model R^2 stats. */ float ** coefts_vol, /* volumes for forward wavelet transform coefs. */ float ** fitts_vol, /* volumes for wavelet filtered time series data */ float ** sgnlts_vol, /* volumes for signal model fit to input data */ float ** errts_vol /* volumes for residual errors */ ) { int ip; /* parameter index */ int p; /* total number of parameters */ int it; /* time point index */ int N; /* number of usable data points */ /*----- Initialize local variables -----*/ p = option_data->p; N = option_data->N; /*----- Saved wavelet coefficients -----*/ for (ip = 0; ip < p; ip++) { coef_vol[ip][iv] = coef[ip]; } /*----- Save full model mean square error -----*/ mse_vol[iv] = mse; /*----- Save regression F-statistic -----*/ ffull_vol[iv] = ffull; /*----- Save R^2 values -----*/ rfull_vol[iv] = rfull; /*----- Save the forward wavelet transform coefficients -----*/ if (coefts_vol != NULL) for (it = 0; it < N; it++) coefts_vol[it][iv] = coefts[it]; /*----- Save the wavelet filtered time series -----*/ if (fitts_vol != NULL) for (it = 0; it < N; it++) fitts_vol[it][iv] = fitts[it]; /*----- Save the fitted signal model time series -----*/ if (sgnlts_vol != NULL) for (it = 0; it < N; it++) sgnlts_vol[it][iv] = sgnlts[it]; /*----- Save the residual errors time series -----*/ if (errts_vol != NULL) for (it = 0; it < N; it++) errts_vol[it][iv] = errts[it]; } /*---------------------------------------------------------------------------*/ /* Calculate the wavelet signal model and associated statistics. */ void calculate_results ( WA_options * option_data, /* wavelet analysis options */ THD_3dim_dataset * dset, /* input 3d+time data set */ THD_3dim_dataset * mask, /* input mask data set */ float * fmri_data, /* input fMRI time series data */ int fmri_length, /* length of fMRI time series */ float ** coef_vol, /* array of volumes of signal model parameters */ float * mse_vol, /* volume of full model mean square error */ float * ffull_vol, /* volume of F-statistic for the full model */ float * rfull_vol, /* volume of R^2 for the full model */ float ** coefts_vol, /* volumes for forward wavelet transform coefs. */ float ** fitts_vol, /* volumes for wavelet filtered time series */ float ** sgnlts_vol, /* volumes for full model fit to input data */ float ** errts_vol /* volumes for residual errors */ ) { float * ts_array = NULL; /* array of measured data for one voxel */ float mask_val[1]; /* value of mask at current voxel */ int f; int p; /* number of parameters in the full model */ int q; /* number of parameters in the baseline model */ float * coef = NULL; /* regression parameters */ float sse_base; /* baseline model error sum of squares */ float sse_full; /* full model error sum of squares */ float ffull; /* full model F-statistic */ float rfull; /* full model R^2 stat. */ int ixyz; /* voxel index */ int nxyz; /* number of voxels per dataset */ int nt; /* number of images in input 3d+time dataset */ int NFirst; /* first image from input 3d+time dataset to use */ int NLast; /* last image from input 3d+time dataset to use */ int N; /* number of usable data points */ int i; /* data point index */ float * coefts = NULL; /* filtered FWT coefficients */ float * fitts = NULL; /* filterd time series */ float * sgnlts = NULL; /* signal model fitted time series */ float * errts = NULL; /* residual error time series */ char * label = NULL; /* string containing stat. summary of results */ /*----- Initialize local variables -----*/ if (option_data->input1D_filename != NULL) { nxyz = 1; nt = fmri_length; } else { nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; nt = DSET_NUM_TIMES (dset); } NFirst = option_data->NFirst; NLast = option_data->NLast; N = option_data->N; f = option_data->f; q = option_data->q; p = option_data->p; /*----- Allocate memory for arrays -----*/ ts_array = (float *) malloc (sizeof(float) * nt); MTEST (ts_array); coef = (float *) malloc (sizeof(float) * p); MTEST (coef); coefts = (float *) malloc (sizeof(float) * N); MTEST (coefts); fitts = (float *) malloc (sizeof(float) * N); MTEST (fitts); sgnlts = (float *) malloc (sizeof(float) * N); MTEST (sgnlts); errts = (float *) malloc (sizeof(float) * N); MTEST (errts); /*----- Loop over all voxels -----*/ for (ixyz = 0; ixyz < nxyz; ixyz++) { /*----- Apply mask? -----*/ if (mask != NULL) { extract_ts_array (mask, ixyz, mask_val); if (mask_val[0] == 0.0) continue; } /*----- Extract Y-data for this voxel -----*/ if (option_data->input1D_filename != NULL) { for (i = 0; i < N; i++) ts_array[i] = fmri_data[i+NFirst]; } else extract_ts_array (dset, ixyz, ts_array); /*----- Perform the wavelet analysis for this voxel-----*/ wavelet_analysis (option_data->wavelet_type, f, option_data->stop_filter, q, option_data->base_filter, p, option_data->sgnl_filter, N, ts_array+NFirst, coef, &sse_base, &sse_full, &ffull, &rfull, coefts, fitts, sgnlts, errts); /*----- Save results for this voxel -----*/ if (option_data->input1D_filename == NULL) save_voxel (option_data, ixyz, coef, sse_full/(N-f-p), ffull, rfull, coefts, fitts, sgnlts, errts, coef_vol, mse_vol, ffull_vol, rfull_vol, coefts_vol, fitts_vol, sgnlts_vol, errts_vol); else { /*----- If only one voxel, save results as .1D files -----*/ ts_fprint ("WA.coefts.1D", N, coefts); ts_fprint ("WA.fitts.1D", N, fitts); ts_fprint ("WA.sgnlts.1D", N, sgnlts); ts_fprint ("WA.errts.1D", N, errts); } /*----- Report results for this voxel -----*/ if ( ((ffull > option_data->fdisp) && (option_data->fdisp >= 0.0)) || (option_data->input1D_filename != NULL) ) { printf ("\n\nResults for Voxel #%d: \n", ixyz); report_results (N, NFirst, f, p, q, option_data->base_filter, option_data->sgnl_filter, coef, sse_base, sse_full, ffull, rfull, &label); printf ("%s \n", label); } } /*----- Loop over voxels -----*/ /*----- Release memory -----*/ free (ts_array); ts_array = NULL; free (coef); coef = NULL; free (coefts); coefts = NULL; free (fitts); fitts = NULL; free (sgnlts); sgnlts = NULL; free (errts); errts = 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 3d+time data set. */ void write_ts_array ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ WA_options * option_data, /* wavelet analysis options */ 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 */ char * input_filename; /* input afni data set file name */ short ** bar = NULL; /* bar[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 */ /*----- Initialize local variables -----*/ input_filename = option_data->input_filename; dset = THD_open_one_dataset (input_filename); nxyz = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz; /*----- allocate memory -----*/ if( option_data->datum != MRI_float ) { bar = (short **) malloc (sizeof(short *) * ts_length); MTEST (bar); } fbuf = (float *) malloc (sizeof(float) * ts_length); MTEST (fbuf); /*-- 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 ) ; { char * commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ; sprintf (label, "Output prefix: %s", output_filename); if( commandline != NULL ) tross_multi_Append_History( new_dset , commandline,label,NULL ) ; else tross_Append_History ( new_dset, label); free(commandline) ; } /*----- 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, option_data->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) ; } /*----- attach bricks to new data set -----*/ for (ib = 0; ib < ts_length; ib++) { if( option_data->datum == MRI_float ) { fbuf[ib] = 1.0; /* factor array is all 1's for float output */ /*----- attach vol_array[ib] to be sub-brick #ib -----*/ mri_fix_data_pointer ((*vol_array)[ib], DSET_BRICK(new_dset,ib)); (*vol_array)[ib] = NULL; /* data is now owned by dataset */ } else { /*----- Set pointer to appropriate volume -----*/ volume = (*vol_array)[ib]; /*----- 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, MRI_short, 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)); } } if( option_data->datum == MRI_float ) *vol_array = NULL; /*----- 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 into %s\n",DSET_BRIKNAME(new_dset)) ; /*----- deallocate memory -----*/ THD_delete_3dim_dataset (new_dset, False); new_dset = NULL ; free (fbuf); fbuf = NULL; } /*---------------------------------------------------------------------------*/ /* Attach one sub-brick to output bucket data set. */ void attach_sub_brick ( THD_3dim_dataset * new_dset, /* output bucket dataset */ int ibrick, /* sub-brick indices */ float * volume, /* volume of floating point data */ int nxyz, /* total number of voxels */ int brick_type, /* indicates statistical type of sub-brick */ char * brick_label, /* character string label for sub-brick */ int dof, int ndof, int ddof, /* degrees of freedom */ short ** bar /* bar[ib] points to data for sub-brick #ib */ ) { const float EPSILON = 1.0e-10; float factor; /* factor is new scale factor for sub-brick #ib */ /*----- allocate memory for output sub-brick -----*/ bar[ibrick] = (short *) malloc (sizeof(short) * nxyz); MTEST (bar[ibrick]); factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume, MRI_short, bar[ibrick]); if (factor < EPSILON) factor = 0.0; else factor = 1.0 / factor; /*----- edit the sub-brick -----*/ EDIT_BRICK_LABEL (new_dset, ibrick, brick_label); EDIT_BRICK_FACTOR (new_dset, ibrick, factor); if (brick_type == FUNC_TT_TYPE) EDIT_BRICK_TO_FITT (new_dset, ibrick, dof); else if (brick_type == FUNC_FT_TYPE) EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof); /*----- attach bar[ib] to be sub-brick #ibrick -----*/ EDIT_substitute_brick (new_dset, ibrick, MRI_short, bar[ibrick]); } /*---------------------------------------------------------------------------*/ /* Attach one sub-brick to output bucket data set (of type MRI_float). */ void attach_float_brick ( THD_3dim_dataset * new_dset, /* output bucket dataset */ int ibrick, /* sub-brick indices */ float * volume, /* volume of floating point data */ int brick_type, /* indicates statistical type of sub-brick */ char * brick_label, /* character string label for sub-brick */ int dof, int ndof, int ddof /* degrees of freedom */ ) { /*----- edit the sub-brick -----*/ EDIT_BRICK_LABEL (new_dset, ibrick, brick_label); EDIT_BRICK_FACTOR (new_dset, ibrick, 1.0); /* all factors are 1.0 here */ if (brick_type == FUNC_TT_TYPE) EDIT_BRICK_TO_FITT (new_dset, ibrick, dof); else if (brick_type == FUNC_FT_TYPE) EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof); /*----- attach bar[ib] to be sub-brick #ibrick -----*/ EDIT_substitute_brick (new_dset, ibrick, MRI_float, volume); } /*---------------------------------------------------------------------------*/ /* Routine to write one bucket data set. All _vol variables have ptrs passed now, to free memory if it is used for the MRI_float dataset sub-bricks. 07 Feb 2006 [rickr] */ void write_bucket_data ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ WA_options * option_data, /* wavelet analysis options */ float *** coef_vol, /* array of volumes of signal model parameters */ float ** mse_vol, /* volume of full model mean square error */ float ** ffull_vol, /* volume of full model F-statistics */ float ** rfull_vol /* volume of full model R^2 statistics */ ) { THD_3dim_dataset * old_dset = NULL; /* prototype dataset */ THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */ char output_prefix[MAX_NAME_LENGTH]; /* prefix name for bucket dataset */ char output_session[MAX_NAME_LENGTH]; /* directory for bucket dataset */ int nbricks; /* number of sub-bricks in bucket dataset */ short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */ int brick_type; /* indicates statistical type of sub-brick */ char brick_label[MAX_NAME_LENGTH]; /* character string label for sub-brick */ int ierror; /* number of errors in editing data */ float * volume; /* volume of floating point data */ int N; /* number of usable data points */ int NFirst; int f; /* number of parameters removed by stop filter */ int p; /* number of parameters in the signal model */ int q; /* number of parameters in the rdcd model */ int nxyz; /* total number of voxels */ int nt; /* number of images in input 3d+time dataset */ int it; int icoef; /* coefficient index */ int ibrick; /* sub-brick index */ int ndof, ddof; /* degrees of freedom */ char label[MAX_NAME_LENGTH]; /* general label for sub-bricks */ /*----- read prototype dataset -----*/ old_dset = THD_open_one_dataset (option_data->input_filename); /*----- Initialize local variables -----*/ nxyz = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz; nt = DSET_NUM_TIMES (old_dset); f = option_data->f; q = option_data->q; p = option_data->p; N = option_data->N; NFirst = option_data->NFirst; printf ("f = %d q = %d p = %d N = %d \n", f, q, p, N); /*----- Calculate number of sub-bricks in the bucket -----*/ nbricks = p * option_data->cout + option_data->rout + option_data->fout + option_data->vout; printf ("nbricks = %d \n", nbricks); strcpy (output_prefix, option_data->bucket_filename); strcpy (output_session, "./"); /*----- allocate memory -----*/ if( option_data->datum != MRI_float ) { bar = (short **) malloc (sizeof(short *) * nbricks); MTEST (bar); } /*-- make an empty copy of prototype dataset, for eventual output --*/ new_dset = EDIT_empty_copy (old_dset); /*----- Record history of dataset -----*/ tross_Copy_History( old_dset , new_dset ) ; tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ; sprintf (label, "Output prefix: %s", output_prefix); tross_Append_History ( new_dset, label); /*----- delete prototype dataset -----*/ THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ; /*----- 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_datum_all, option_data->datum, 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 bucket 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); } /*----- Attach individual sub-bricks to the bucket dataset -----*/ ibrick = 0; /*----- Full model wavelet coefficients -----*/ if (option_data->cout) { /*----- User can choose to place full model stats first -----*/ if (option_data->stat_first) ibrick += option_data->vout + option_data->rout + option_data->fout; icoef = 0; for (it = 0; it < N; it++) { if (option_data->sgnl_filter[it]) { /*----- Wavelet coefficient -----*/ brick_type = FUNC_FIM_TYPE; { int band, mintr, maxtr; if (it == 0) { band = -1; mintr = 0; maxtr = N-1; } else { band = my_log2(it); mintr = (it - powerof2(band)) * powerof2(my_log2(N)-band); maxtr = mintr + powerof2(my_log2(N)-band) - 1; } mintr += NFirst; maxtr += NFirst; if (option_data->base_filter[it]) sprintf (brick_label, "B(%2d)[%3d,%3d]", band, mintr, maxtr); else sprintf (brick_label, "S(%2d)[%3d,%3d]", band, mintr, maxtr); } volume = (*coef_vol)[icoef]; if( option_data->datum == MRI_float ) { attach_float_brick (new_dset, ibrick, volume, brick_type, brick_label, 0, 0, 0); (*coef_vol)[icoef] = NULL; /* dataset owns this now */ } else attach_sub_brick (new_dset, ibrick, volume, nxyz, brick_type, brick_label, 0, 0, 0, bar); ibrick++; icoef++; } } /* End loop over Wavelet coefficients */ if( option_data->datum == MRI_float ) /* then free array, also */ { free( *coef_vol ); *coef_vol = NULL; } } /* if (option_data->cout) */ /*----- Statistics for full model -----*/ if (option_data->stat_first) ibrick = 0; /*----- Full model MSE -----*/ if (option_data->vout) { brick_type = FUNC_FIM_TYPE; sprintf (brick_label, "Full MSE"); volume = *mse_vol; if( option_data->datum == MRI_float ) { attach_float_brick (new_dset, ibrick, volume, brick_type, brick_label, 0, 0, 0); *mse_vol = NULL; /* since the dataset owns this now */ } else attach_sub_brick (new_dset, ibrick, volume, nxyz, brick_type, brick_label, 0, 0, 0, bar); ibrick++; } /*----- Full model R^2 -----*/ if (option_data->rout) { brick_type = FUNC_THR_TYPE; sprintf (brick_label, "Full R^2"); volume = *rfull_vol; if( option_data->datum == MRI_float ) { attach_float_brick (new_dset, ibrick, volume, brick_type, brick_label, 0, 0, 0); *rfull_vol = NULL; /* since the dataset owns this now */ } else attach_sub_brick (new_dset, ibrick, volume, nxyz, brick_type, brick_label, 0, 0, 0, bar); ibrick++; } /*----- Full model F-stat -----*/ if (option_data->fout) { brick_type = FUNC_FT_TYPE; ndof = p - q; ddof = N - f - p; sprintf (brick_label, "Full F-stat"); volume = *ffull_vol; if( option_data->datum == MRI_float ) { attach_float_brick (new_dset, ibrick, volume, brick_type, brick_label, 0, ndof, ddof); *ffull_vol = NULL; /* since the dataset owns this now */ } else attach_sub_brick (new_dset, ibrick, volume, nxyz, brick_type, brick_label, 0, ndof, ddof, bar); ibrick++; } /*----- write bucket data set -----*/ THD_load_statistics (new_dset); THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ; fprintf(stderr,"Writing `bucket' dataset "); fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset)); /*----- deallocate memory -----*/ THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; } /*---------------------------------------------------------------------------*/ /* Write out the user requested output files. Each item has address passed, so that data used for MRI_float datasets can have its pointer(s) cleared. 7 Feb 2006 [rickr] */ void output_results ( int argc, /* number of input arguments */ char ** argv, /* array of input arguments */ WA_options * option_data, /* wavelet analysis options */ float *** coef_vol, /* array of volumes of signal model parameters */ float ** mse_vol, /* volume of full model mean square error */ float ** ffull_vol, /* volume of full model F-statistics */ float ** rfull_vol, /* volume of full model R^2 statistics */ float *** coefts_vol, /* volumes for forward wavelet transform coefs. */ float *** fitts_vol, /* volumes for wavelet filtered time series */ float *** sgnlts_vol, /* volumes for signal model time series */ float *** errts_vol /* volumes for residual errors */ ) { int N; /* number of usable data points */ /*----- Initialize local variables -----*/ N = option_data->N; /*----- Write the bucket dataset -----*/ if (option_data->bucket_filename != NULL) write_bucket_data (argc, argv, option_data, coef_vol, mse_vol, ffull_vol, rfull_vol); /*----- Write the forward wavelet transform coefficients -----*/ if (option_data->coefts_filename != NULL) write_ts_array (argc, argv, option_data, N, coefts_vol, option_data->coefts_filename); /*----- Write the wavelet filtered 3d+time dataset -----*/ if (option_data->fitts_filename != NULL) write_ts_array (argc, argv, option_data, N, fitts_vol, option_data->fitts_filename); /*----- Write the signal model 3d+time dataset -----*/ if (option_data->sgnlts_filename != NULL) write_ts_array (argc, argv, option_data, N, sgnlts_vol, option_data->sgnlts_filename); /*----- Write the residual errors 3d+time dataset -----*/ if (option_data->errts_filename != NULL) write_ts_array (argc, argv, option_data, N, errts_vol, option_data->errts_filename); } /*---------------------------------------------------------------------------*/ void terminate_program ( WA_options ** option_data, /* wavelet analysis options */ float ** fmri_data, /* input fMRI time series data */ float *** coef_vol, /* array of volumes of signal model parameters */ float ** mse_vol, /* volume of full model mean square error */ float ** ffull_vol, /* volume of full model F-statistics */ float ** rfull_vol, /* volume of full model R^2 statistics */ float *** coefts_vol, /* volumes for forward wavelet transform coefs. */ float *** fitts_vol, /* volumes for wavelet filtered time series data */ float *** sgnlts_vol, /* volumes for signal model fit to input data */ float *** errts_vol /* volumes for residual errors */ ) { int p; /* number of parameters in full model */ int ip; /* parameter index */ int it; /* time index */ int N; /* number of usable data points */ /*----- Initialize local variables -----*/ p = (*option_data)->p; N = (*option_data)->N; /*----- Deallocate memory for option data -----*/ if ((*option_data)->stop_filter != NULL) free ((*option_data)->stop_filter); if ((*option_data)->base_filter != NULL) free ((*option_data)->base_filter); if ((*option_data)->sgnl_filter != NULL) free ((*option_data)->sgnl_filter); if ((*option_data)->input_filename != NULL) free ((*option_data)->input_filename); if ((*option_data)->mask_filename != NULL) free ((*option_data)->mask_filename); if ((*option_data)->input1D_filename != NULL) free ((*option_data)->input1D_filename); if ((*option_data)->bucket_filename != NULL) free ((*option_data)->bucket_filename); if ((*option_data)->coefts_filename != NULL) free ((*option_data)->coefts_filename); if ((*option_data)->fitts_filename != NULL) free ((*option_data)->fitts_filename); if ((*option_data)->sgnlts_filename != NULL) free ((*option_data)->sgnlts_filename); if ((*option_data)->errts_filename != NULL) free ((*option_data)->errts_filename); free (*option_data); *option_data = NULL; /*----- Deallocate memory for fMRI time series data -----*/ if (*fmri_data != NULL) { free (*fmri_data); *fmri_data = NULL; } /*----- Deallocate space for volume data -----*/ if (*coef_vol != NULL) { for (ip = 0; ip < p; ip++) { if ((*coef_vol)[ip] != NULL) { free ((*coef_vol)[ip]); (*coef_vol)[ip] = NULL; } } free (*coef_vol); *coef_vol = NULL; } if (*mse_vol != NULL) { free (*mse_vol); *mse_vol = NULL; } if (*ffull_vol != NULL) { free (*ffull_vol); *ffull_vol = NULL; } if (*rfull_vol != NULL) { free (*rfull_vol); *rfull_vol = NULL; } /*----- Deallocate space for forward wavelet transform coefficients -----*/ if (*coefts_vol != NULL) { for (it = 0; it < N; it++) { if ((*coefts_vol)[it] != NULL) { free ((*coefts_vol)[it]); (*coefts_vol)[it] = NULL; } } free (*coefts_vol); *coefts_vol = NULL; } /*----- Deallocate space for wavelet filtered time series -----*/ if (*fitts_vol != NULL) { for (it = 0; it < N; it++) { if ((*fitts_vol)[it] != NULL) { free ((*fitts_vol)[it]); (*fitts_vol)[it] = NULL; } } free (*fitts_vol); *fitts_vol = NULL; } /*----- Deallocate space for signal model time series -----*/ if (*sgnlts_vol != NULL) { for (it = 0; it < N; it++) { if ((*sgnlts_vol)[it] != NULL) { free ((*sgnlts_vol)[it]); (*sgnlts_vol)[it] = NULL; } } free (*sgnlts_vol); *sgnlts_vol = NULL; } /*----- Deallocate space for residual errors time series -----*/ if (*errts_vol != NULL) { for (it = 0; it < N; it++) { if ((*errts_vol)[it] != NULL) { free ((*errts_vol)[it]); (*errts_vol)[it] = NULL; } } free (*errts_vol); *errts_vol = NULL; } } /*---------------------------------------------------------------------------*/ int main ( int argc, /* number of input arguments */ char ** argv /* array of input arguments */ ) { WA_options * option_data; /* wavelet analysis options */ THD_3dim_dataset * dset_time = NULL; /* input 3d+time data set */ THD_3dim_dataset * mask_dset = NULL; /* input mask data set */ float * fmri_data = NULL; /* input fMRI time series data */ int fmri_length; /* length of fMRI time series */ float ** coef_vol = NULL; /* array of volumes for model parameters */ float * mse_vol = NULL; /* volume of full model mean square error */ float * ffull_vol = NULL; /* volume of full model F-statistics */ float * rfull_vol = NULL; /* volume of full model R^2 stats. */ float ** coefts_vol = NULL; /* volumes for wavelet transform coefs.*/ float ** fitts_vol = NULL; /* volumes for wavelet filtered time series */ float ** sgnlts_vol = NULL; /* volumes for signal model fit to input data */ float ** errts_vol = NULL; /* volumes for residual errors */ /*----- 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 PRINT_VERSION("3dWavelets") ; AUTHOR(PROGRAM_AUTHOR); mainENTRY("3dWavelets main") ; machdep() ; /*----- Program initialization -----*/ initialize_program (argc, argv, &option_data, &dset_time, &mask_dset, &fmri_data, &fmri_length, &coef_vol, &mse_vol, &ffull_vol, &rfull_vol, &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol); /*----- Perform wavelet analysis -----*/ calculate_results (option_data, dset_time, mask_dset, fmri_data, fmri_length, coef_vol, mse_vol, ffull_vol, rfull_vol, coefts_vol, fitts_vol, sgnlts_vol, errts_vol); /*----- Deallocate memory for input datasets -----*/ if (dset_time != NULL) { THD_delete_3dim_dataset (dset_time, False); dset_time = NULL; } if (mask_dset != NULL) { THD_delete_3dim_dataset (mask_dset, False); mask_dset = NULL; } /*----- Write requested output files -----*/ /* (pass addresses, in case data is stolen for datasets) 7 Feb 2006 [rickr] */ if (option_data->input1D_filename == NULL) output_results (argc, argv, option_data, &coef_vol, &mse_vol, &ffull_vol, &rfull_vol, &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol); /*----- Terminate program -----*/ terminate_program (&option_data, &fmri_data, &coef_vol, &mse_vol, &ffull_vol, &rfull_vol, &coefts_vol, &fitts_vol, &sgnlts_vol, &errts_vol); exit(0); }