/***************************************************************************** Details for this model are in prf_common_circular.c, #included. This file now just holds functionality that is specific to the basic 4-parameter method. R. Reynolds June, 2014 ******************************************************************************/ #include "NLfit_model.h" static char * g_model_ver = "model_conv_PRF, version 2.4, 19 Jun, 2018"; /* prototypes */ static int signal_model( float * , int , float ** , float *, int ); static int model_help(void); /* suck in all of the common functions */ #include "prf_common_circular.c" /*-----------------------------------------------------------------------*/ DEFINE_MODEL_PROTOTYPE MODEL_interface * initialize_model () { MODEL_interface * mi = NULL; /*----- first, see if the user wants help -----*/ if ( AFNI_yesenv("AFNI_MODEL_HELP_CONV_PRF") || AFNI_yesenv("AFNI_MODEL_HELP_ALL") ) model_help(); /*----- allocate memory space for model interface -----*/ mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface)); /*----- name of this model -----*/ strcpy (mi->label, "Conv_PRF"); /*----- this is a signal model -----*/ mi->model_type = MODEL_SIGNAL_TYPE; /*----- number of parameters in the model -----*/ mi->params = 4; /*----- parameter labels -----*/ strcpy (mi->plabel[0], "Amp"); strcpy (mi->plabel[1], "X"); strcpy (mi->plabel[2], "Y"); strcpy (mi->plabel[3], "Sigma"); /*----- minimum and maximum parameter constraints -----*/ /* amplitude, x/y ranges, and sigma range */ mi->min_constr[0] = -10.0; mi->max_constr[0] = 10.0; mi->min_constr[1] = -1.0; mi->max_constr[1] = 1.0; mi->min_constr[2] = -1.0; mi->max_constr[2] = 1.0; mi->min_constr[3] = 0.0; mi->max_constr[3] = 1.0; /*----- function which implements the model -----*/ mi->call_func = conv_model; return (mi); } /*----------------------------------------------------------------------*/ /* Routine to calculate the time series to (hopefully) fit the data. Definition of model parameters (gs[2] > 0) gs[0] = Amp = amplitude gs[1] = x0 = x-coordinate of gaussian center gs[2] = y0 = y-coordinate of gaussian center gs[3] = sigma = "width" of gaussian curve For each TR, integrate g(x,y) over stim aperture dset. g(x,y) = e^-[((x-x0)^2+(y-y0)^2)/(2*sigma^2)] The resulting returned time series will be convolved in the parent function. */ static int signal_model ( float * gs, /* parameters for signal model */ int ts_length, /* length of time series data */ float ** x_array, /* independent variable matrix */ float * ts_array, /* estimated signal model time series */ int debug /* make some noise */ ) { float A, x, y, sigma;/* model params */ int maxind; /* largest dimension */ int tmpmax; /* assign parameters */ A = gs[0]; x = gs[1]; y = gs[2]; sigma = gs[3]; if( debug ) { fprintf(stderr, "-d model_conv_PRF parameters: " "A = %f, x = %f, y = %f, sigma = %f\n" " nx = %d, nz = %d, nvals = %d, ts_len = %d\n", A, x, y, sigma, DSET_NX(g_saset), DSET_NZ(g_saset), DSET_NVALS(g_saset), ts_length); show_malloc_stats("signal model"); } if( ! ISVALID_3DIM_DATASET(g_saset) ) return 0; /* possibly restrict the length, via nx if on grid, else nt */ /* (was just NX) 30 Aug 2017 */ maxind = ts_length; if( genv_on_grid ) tmpmax = DSET_NX(g_saset); else tmpmax = DSET_NVALS(g_saset); if( maxind > tmpmax ) maxind = tmpmax; if( maxind == 0 ) return 0; if( debug ) fprintf( stderr,"-d NT orig=%d, applied=%d\n", ts_length, maxind); /* time array must be ordered according to stim dset */ if( genv_on_grid ) /* scale data directly from grid */ get_signal_from_grid(ts_array, maxind, g_saset, x, y, sigma, A, debug); else get_signal_computed(ts_array, maxind, g_saset, x, y, sigma, A, debug); if( debug ) disp_floats("+d signal model result : ", ts_array, ts_length); return maxind; } /*----------------------------------------------------------------------*/ static int model_help(void) { printf( "----------------------------------------------------------------------\n" "PRF - population receptive field (in visual cortex)\n" "\n" " This model is from the paper:\n" "\n" " Population receptive field estimates in human visual cortex\n" " NeuroImage 39 (2008) 647-660\n" " Serge O. Dumoulin, Brian A. Wandell\n" "\n" " The model is made from parameters A, x0, y0, sigma, and from stimulus\n" " time series input (visual field masks over time) by:\n" "\n" " 1. compute a Gaussian curve centered at x0, y0 of with spread sigma\n" " g(x,y) = e^-( [(x-x0)^2+(y-y0)^2] / (2*sigma^2) )\n" " 2. multiply this 2-D image by each 2-D stimulus mask image\n" " 3. convolve the result with an ideal HRF\n" " 4. scale by the amplitude A\n" "\n" " Currently, x0, y0, and sigma are limited to [-1,1], which the stimulus\n" " images are evaluated on. This use may be altered in the future.\n" "\n" "--------------------------------------------------\n" "To use this model function:\n" "\n" " 1. Generate the stimulus time series (currently, images must be square).\n" "\n" " This should be a 2D+time dataset of visual stimuli over time. They\n" " are viewed as binary masks by the model function.\n" "\n" " * If results are computed on a restricted grid (which is much faster\n" " and is the default (see AFNI_MODEL_PRF_ON_GRID)), the resolution of\n" " those X,Y results will come directly from this stimulus dataset.\n" " It might be reasonable to have this be 100 or 200 (or 101 or 201)\n" " voxels on a side.\n" "\n" " * The amount of memory used for the precomputation should be the size\n" " of this dataset (in float format) times AFNI_MODEL_PRF_SIGMA_NSTEPS.\n" " It is converted to floats because it is blurred internally.\n" " The default AFNI_MODEL_PRF_SIGMA_NSTEPS is 100.\n" "\n" " 2. Scale and demean the input EPI time series data.\n" "\n" " Scaling is done to put the amplitude values into a reasonable (i.e.\n" " expected) range, such as by scaling it to a fraction of the mean\n" " (or maybe twice that).\n" "\n" " Setting the mean to zero is done so that no baseline modeling is\n" " needed (though it might be good to model drifts in the future).\n" "\n" " 3. Generate a convolution reference time series, such as one for GAM.\n" " This should be on the same TR grid, which is 2 in this example.\n" "\n" " 3dDeconvolve -nodata 10 2 -polort -1 \\\n" " -num_stimts 1 -stim_times 1 '1D:0' GAM \\\n" " -x1D conv.ref.GAM.1D\n" "\n" " 4. Set up environment variables to control execution:\n" "\n" " setenv AFNI_CONVMODEL_REF conv.ref.GAM.1D\n" " setenv AFNI_MODEL_PRF_STIM_DSET stim.144.bmask.resam+tlrc\n" "\n" " 5. And execute:\n" "\n" " 3dNLfim -input epi.scale.demean+tlrc \\\n" " -noise Zero \\\n" " -signal Conv_PRF \\\n" " -sconstr 0 -10.0 10.0 \\\n" " -sconstr 1 -1.0 1.0 \\\n" " -sconstr 2 -1.0 1.0 \\\n" " -sconstr 3 0.0 1.0 \\\n" " -BOTH \\\n" " -nrand 10000 \\\n" " -nbest 5 \\\n" " -bucket 0 buck.PRF \\\n" " -snfit snfit.PRF\n" "\n" "--------------------------------------------------\n" "environment variables:\n" "\n" " -----------------------------------\n" " required:\n" "\n" " AFNI_CONVMODEL_REF : specify convolution reference file\n" "\n" " e.g. setenv AFNI_CONVMODEL_REF conv.ref.GAM.1D\n" "\n" " The file this specifies should contain a (short?) impulse\n" " response function, such as made in step #3, above.\n" "\n" " AFNI_MODEL_PRF_STIM_DSET : specify visual stimulus dataset\n" "\n" " e.g. setenv AFNI_MODEL_PRF_STIM_DSETstim.144.bmask.resam+tlrc\n" "\n" " This should be a 2D+time dataset of stimulus images over TRs.\n" " It will be converted to a byte mask over the visual field.\n" "\n" " -----------------------------------\n" " optional (for use with pre-computed grid):\n" "\n" " AFNI_MODEL_PRF_ON_GRID : Y/N - use pre-computed solutions\n" "\n" " e.g. setenv AFNI_MODEL_PRF_ON_GRID NO\n" " e.g. default YES\n" "\n" " Recommended.\n" "\n" " When set, the model function will actually pre-compute all possible\n" " (unscaled) fit solutions on the first pass. Since all of these\n" " parameters have a smooth effect on the result, this method should\n" " be sufficient.\n" "\n" " Note that the resolution of x0, y0 parameters comes directly from\n" " the stimulus dataset (AFNI_MODEL_PRF_STIM_DSET), while the sigma\n" " resolution comes from the maximum (AFNI_MODEL_PRF_SIGMA_MAX) and\n" " the number of computed values (AFNI_MODEL_PRF_SIGMA_NSTEPS).\n" "\n" " The more voxels to solve for in the input EPI, the more useful this\n" " is. For a single voxel, it is slow. For a large dataset, it can\n" " speed up the solution by a factor of 1000.\n" "\n" " AFNI_MODEL_PRF_SIGMA_MAX : specify maximum allowable sigma\n" "\n" " e.g. setenv AFNI_MODEL_PRF_SIGMA_MAX 2.0\n" " e.g. default 1.0\n" "\n" " Applies directly to AFNI_MODEL_PRF_ON_GRID.\n" "\n" " Use this variable to set the maximum pre-computed sigma.\n" " This should probably match the sconstr value for sigma.\n" "\n" " AFNI_MODEL_PRF_SIGMA_NSTEPS : specify number of pre-computed sigmas\n" "\n" " e.g. setenv AFNI_MODEL_PRF_SIGMA_NSTEPS 50\n" " e.g. default 100\n" "\n" " Applies directly to AFNI_MODEL_PRF_ON_GRID.\n" "\n" " Use this variable to set the number of pre-computed sigma values.\n" " Note that the resolution of pre-computed sigma values will be the\n" " ratio: AFNI_MODEL_PRF_SIGMA_MAX/AFNI_MODEL_PRF_SIGMA_NSTEPS.\n" "\n" " AFNI_MODEL_PRF_RAM_STATS VAL : request display of RAM usage\n" "\n" " e.g. setenv AFNI_MODEL_PRF_RAM_STATS Y\n" " e.g. default N\n" "\n" " Use this variable to control display of RAM usage. By default,\n" " is it off. VAL can be one of:\n" "\n" " Y : yes, show all information\n" " N : no [default], show no information\n" " MALLOC : show only MALLOC information\n" " PS : show only PS information\n" " ALL : same as Y\n" " WAIT : same as Y, and wait after output\n" "\n" " -----------------------------------\n" " helpful:\n" "\n" " AFNI_MODEL_HELP_CONV_PRF : Y/N - output this help\n" "\n" " e.g. setenv AFNI_MODEL_HELP_CONV_PRF YES\n" "\n" " When set, the model initialization function will output this help.\n" "\n" " Consider:\n" "\n" " 3dNLfim -signal Conv_PRF\n" "\n" " or more directly (without setenv):\n" "\n" " 3dNLfim -DAFNI_MODEL_HELP_CONV_PRF=Y -signal Conv_PRF\n" "\n" " AFNI_MODEL_DEBUG : specify debug/verbosity level\n" "\n" " e.g. setenv AFNI_MODEL_DEBUG 2\n" "\n" " Be more verbose. Valid levels are from 0 to 3, currently.\n" "\n" " AFNI_MODEL_DITER : specify debug iteration\n" "\n" " e.g. setenv AFNI_MODEL_DITER 999\n" "\n" " Get extra debug info at some iteration.\n" "\n" "----------------------------------------------------------------------\n" " Written for E Silson and C Baker.\n" "\n" " R. Reynolds 27 June, 2014\n" "----------------------------------------------------------------------\n" ); return 0 ; }