1 /*****************************************************************************
2 Details for this model are in prf_common_circular.c, #included.
3
4 This file now just holds functionality that is specific to the
5 basic 4-parameter method.
6
7 R. Reynolds June, 2014
8 ******************************************************************************/
9
10 #include "NLfit_model.h"
11
12 static char * g_model_ver = "model_conv_PRF, version 2.4, 19 Jun, 2018";
13
14 /* prototypes */
15 static int signal_model( float * , int , float ** , float *, int );
16 static int model_help(void);
17
18 /* suck in all of the common functions */
19 #include "prf_common_circular.c"
20
21 /*-----------------------------------------------------------------------*/
22
23 DEFINE_MODEL_PROTOTYPE
24
initialize_model()25 MODEL_interface * initialize_model ()
26 {
27 MODEL_interface * mi = NULL;
28
29 /*----- first, see if the user wants help -----*/
30 if ( AFNI_yesenv("AFNI_MODEL_HELP_CONV_PRF") ||
31 AFNI_yesenv("AFNI_MODEL_HELP_ALL") ) model_help();
32
33 /*----- allocate memory space for model interface -----*/
34
35 mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));
36
37 /*----- name of this model -----*/
38
39 strcpy (mi->label, "Conv_PRF");
40
41 /*----- this is a signal model -----*/
42
43 mi->model_type = MODEL_SIGNAL_TYPE;
44
45 /*----- number of parameters in the model -----*/
46
47 mi->params = 4;
48
49 /*----- parameter labels -----*/
50
51 strcpy (mi->plabel[0], "Amp");
52 strcpy (mi->plabel[1], "X");
53 strcpy (mi->plabel[2], "Y");
54 strcpy (mi->plabel[3], "Sigma");
55
56 /*----- minimum and maximum parameter constraints -----*/
57
58 /* amplitude, x/y ranges, and sigma range */
59 mi->min_constr[0] = -10.0; mi->max_constr[0] = 10.0;
60
61 mi->min_constr[1] = -1.0; mi->max_constr[1] = 1.0;
62 mi->min_constr[2] = -1.0; mi->max_constr[2] = 1.0;
63
64 mi->min_constr[3] = 0.0; mi->max_constr[3] = 1.0;
65
66 /*----- function which implements the model -----*/
67 mi->call_func = conv_model;
68
69 return (mi);
70 }
71
72
73 /*----------------------------------------------------------------------*/
74 /*
75 Routine to calculate the time series to (hopefully) fit the data.
76
77 Definition of model parameters (gs[2] > 0)
78
79 gs[0] = Amp = amplitude
80 gs[1] = x0 = x-coordinate of gaussian center
81 gs[2] = y0 = y-coordinate of gaussian center
82 gs[3] = sigma = "width" of gaussian curve
83
84 For each TR, integrate g(x,y) over stim aperture dset.
85
86 g(x,y) = e^-[((x-x0)^2+(y-y0)^2)/(2*sigma^2)]
87
88 The resulting returned time series will be convolved in the
89 parent function.
90 */
signal_model(float * gs,int ts_length,float ** x_array,float * ts_array,int debug)91 static int signal_model
92 (
93 float * gs, /* parameters for signal model */
94 int ts_length, /* length of time series data */
95 float ** x_array, /* independent variable matrix */
96 float * ts_array, /* estimated signal model time series */
97 int debug /* make some noise */
98 )
99 {
100 float A, x, y, sigma;/* model params */
101 int maxind; /* largest dimension */
102 int tmpmax;
103
104
105 /* assign parameters */
106 A = gs[0];
107 x = gs[1]; y = gs[2]; sigma = gs[3];
108
109 if( debug ) {
110 fprintf(stderr, "-d model_conv_PRF parameters: "
111 "A = %f, x = %f, y = %f, sigma = %f\n"
112 " nx = %d, nz = %d, nvals = %d, ts_len = %d\n",
113 A, x, y, sigma,
114 DSET_NX(g_saset), DSET_NZ(g_saset), DSET_NVALS(g_saset),
115 ts_length);
116 show_malloc_stats("signal model");
117 }
118
119 if( ! ISVALID_3DIM_DATASET(g_saset) ) return 0;
120
121 /* possibly restrict the length, via nx if on grid, else nt */
122 /* (was just NX) 30 Aug 2017 */
123 maxind = ts_length;
124
125 if( genv_on_grid ) tmpmax = DSET_NX(g_saset);
126 else tmpmax = DSET_NVALS(g_saset);
127
128 if( maxind > tmpmax ) maxind = tmpmax;
129 if( maxind == 0 ) return 0;
130
131 if( debug )
132 fprintf( stderr,"-d NT orig=%d, applied=%d\n", ts_length, maxind);
133
134 /* time array must be ordered according to stim dset */
135 if( genv_on_grid ) /* scale data directly from grid */
136 get_signal_from_grid(ts_array, maxind, g_saset, x, y, sigma, A, debug);
137 else
138 get_signal_computed(ts_array, maxind, g_saset, x, y, sigma, A, debug);
139
140 if( debug )
141 disp_floats("+d signal model result : ", ts_array, ts_length);
142
143 return maxind;
144 }
145
146
147 /*----------------------------------------------------------------------*/
model_help(void)148 static int model_help(void)
149 {
150 printf(
151 "----------------------------------------------------------------------\n"
152 "PRF - population receptive field (in visual cortex)\n"
153 "\n"
154 " This model is from the paper:\n"
155 "\n"
156 " Population receptive field estimates in human visual cortex\n"
157 " NeuroImage 39 (2008) 647-660\n"
158 " Serge O. Dumoulin, Brian A. Wandell\n"
159 "\n"
160 " The model is made from parameters A, x0, y0, sigma, and from stimulus\n"
161 " time series input (visual field masks over time) by:\n"
162 "\n"
163 " 1. compute a Gaussian curve centered at x0, y0 of with spread sigma\n"
164 " g(x,y) = e^-( [(x-x0)^2+(y-y0)^2] / (2*sigma^2) )\n"
165 " 2. multiply this 2-D image by each 2-D stimulus mask image\n"
166 " 3. convolve the result with an ideal HRF\n"
167 " 4. scale by the amplitude A\n"
168 "\n"
169 " Currently, x0, y0, and sigma are limited to [-1,1], which the stimulus\n"
170 " images are evaluated on. This use may be altered in the future.\n"
171 "\n"
172 "--------------------------------------------------\n"
173 "To use this model function:\n"
174 "\n"
175 " 1. Generate the stimulus time series (currently, images must be square).\n"
176 "\n"
177 " This should be a 2D+time dataset of visual stimuli over time. They\n"
178 " are viewed as binary masks by the model function.\n"
179 "\n"
180 " * If results are computed on a restricted grid (which is much faster\n"
181 " and is the default (see AFNI_MODEL_PRF_ON_GRID)), the resolution of\n"
182 " those X,Y results will come directly from this stimulus dataset.\n"
183 " It might be reasonable to have this be 100 or 200 (or 101 or 201)\n"
184 " voxels on a side.\n"
185 "\n"
186 " * The amount of memory used for the precomputation should be the size\n"
187 " of this dataset (in float format) times AFNI_MODEL_PRF_SIGMA_NSTEPS.\n"
188 " It is converted to floats because it is blurred internally.\n"
189 " The default AFNI_MODEL_PRF_SIGMA_NSTEPS is 100.\n"
190 "\n"
191 " 2. Scale and demean the input EPI time series data.\n"
192 "\n"
193 " Scaling is done to put the amplitude values into a reasonable (i.e.\n"
194 " expected) range, such as by scaling it to a fraction of the mean\n"
195 " (or maybe twice that).\n"
196 "\n"
197 " Setting the mean to zero is done so that no baseline modeling is\n"
198 " needed (though it might be good to model drifts in the future).\n"
199 "\n"
200 " 3. Generate a convolution reference time series, such as one for GAM.\n"
201 " This should be on the same TR grid, which is 2 in this example.\n"
202 "\n"
203 " 3dDeconvolve -nodata 10 2 -polort -1 \\\n"
204 " -num_stimts 1 -stim_times 1 '1D:0' GAM \\\n"
205 " -x1D conv.ref.GAM.1D\n"
206 "\n"
207 " 4. Set up environment variables to control execution:\n"
208 "\n"
209 " setenv AFNI_CONVMODEL_REF conv.ref.GAM.1D\n"
210 " setenv AFNI_MODEL_PRF_STIM_DSET stim.144.bmask.resam+tlrc\n"
211 "\n"
212 " 5. And execute:\n"
213 "\n"
214 " 3dNLfim -input epi.scale.demean+tlrc \\\n"
215 " -noise Zero \\\n"
216 " -signal Conv_PRF \\\n"
217 " -sconstr 0 -10.0 10.0 \\\n"
218 " -sconstr 1 -1.0 1.0 \\\n"
219 " -sconstr 2 -1.0 1.0 \\\n"
220 " -sconstr 3 0.0 1.0 \\\n"
221 " -BOTH \\\n"
222 " -nrand 10000 \\\n"
223 " -nbest 5 \\\n"
224 " -bucket 0 buck.PRF \\\n"
225 " -snfit snfit.PRF\n"
226 "\n"
227 "--------------------------------------------------\n"
228 "environment variables:\n"
229 "\n"
230 " -----------------------------------\n"
231 " required:\n"
232 "\n"
233 " AFNI_CONVMODEL_REF : specify convolution reference file\n"
234 "\n"
235 " e.g. setenv AFNI_CONVMODEL_REF conv.ref.GAM.1D\n"
236 "\n"
237 " The file this specifies should contain a (short?) impulse\n"
238 " response function, such as made in step #3, above.\n"
239 "\n"
240 " AFNI_MODEL_PRF_STIM_DSET : specify visual stimulus dataset\n"
241 "\n"
242 " e.g. setenv AFNI_MODEL_PRF_STIM_DSETstim.144.bmask.resam+tlrc\n"
243 "\n"
244 " This should be a 2D+time dataset of stimulus images over TRs.\n"
245 " It will be converted to a byte mask over the visual field.\n"
246 "\n"
247 " -----------------------------------\n"
248 " optional (for use with pre-computed grid):\n"
249 "\n"
250 " AFNI_MODEL_PRF_ON_GRID : Y/N - use pre-computed solutions\n"
251 "\n"
252 " e.g. setenv AFNI_MODEL_PRF_ON_GRID NO\n"
253 " e.g. default YES\n"
254 "\n"
255 " Recommended.\n"
256 "\n"
257 " When set, the model function will actually pre-compute all possible\n"
258 " (unscaled) fit solutions on the first pass. Since all of these\n"
259 " parameters have a smooth effect on the result, this method should\n"
260 " be sufficient.\n"
261 "\n"
262 " Note that the resolution of x0, y0 parameters comes directly from\n"
263 " the stimulus dataset (AFNI_MODEL_PRF_STIM_DSET), while the sigma\n"
264 " resolution comes from the maximum (AFNI_MODEL_PRF_SIGMA_MAX) and\n"
265 " the number of computed values (AFNI_MODEL_PRF_SIGMA_NSTEPS).\n"
266 "\n"
267 " The more voxels to solve for in the input EPI, the more useful this\n"
268 " is. For a single voxel, it is slow. For a large dataset, it can\n"
269 " speed up the solution by a factor of 1000.\n"
270 "\n"
271 " AFNI_MODEL_PRF_SIGMA_MAX : specify maximum allowable sigma\n"
272 "\n"
273 " e.g. setenv AFNI_MODEL_PRF_SIGMA_MAX 2.0\n"
274 " e.g. default 1.0\n"
275 "\n"
276 " Applies directly to AFNI_MODEL_PRF_ON_GRID.\n"
277 "\n"
278 " Use this variable to set the maximum pre-computed sigma.\n"
279 " This should probably match the sconstr value for sigma.\n"
280 "\n"
281 " AFNI_MODEL_PRF_SIGMA_NSTEPS : specify number of pre-computed sigmas\n"
282 "\n"
283 " e.g. setenv AFNI_MODEL_PRF_SIGMA_NSTEPS 50\n"
284 " e.g. default 100\n"
285 "\n"
286 " Applies directly to AFNI_MODEL_PRF_ON_GRID.\n"
287 "\n"
288 " Use this variable to set the number of pre-computed sigma values.\n"
289 " Note that the resolution of pre-computed sigma values will be the\n"
290 " ratio: AFNI_MODEL_PRF_SIGMA_MAX/AFNI_MODEL_PRF_SIGMA_NSTEPS.\n"
291 "\n"
292 " AFNI_MODEL_PRF_RAM_STATS VAL : request display of RAM usage\n"
293 "\n"
294 " e.g. setenv AFNI_MODEL_PRF_RAM_STATS Y\n"
295 " e.g. default N\n"
296 "\n"
297 " Use this variable to control display of RAM usage. By default,\n"
298 " is it off. VAL can be one of:\n"
299 "\n"
300 " Y : yes, show all information\n"
301 " N : no [default], show no information\n"
302 " MALLOC : show only MALLOC information\n"
303 " PS : show only PS information\n"
304 " ALL : same as Y\n"
305 " WAIT : same as Y, and wait after output\n"
306 "\n"
307 " -----------------------------------\n"
308 " helpful:\n"
309 "\n"
310 " AFNI_MODEL_HELP_CONV_PRF : Y/N - output this help\n"
311 "\n"
312 " e.g. setenv AFNI_MODEL_HELP_CONV_PRF YES\n"
313 "\n"
314 " When set, the model initialization function will output this help.\n"
315 "\n"
316 " Consider:\n"
317 "\n"
318 " 3dNLfim -signal Conv_PRF\n"
319 "\n"
320 " or more directly (without setenv):\n"
321 "\n"
322 " 3dNLfim -DAFNI_MODEL_HELP_CONV_PRF=Y -signal Conv_PRF\n"
323 "\n"
324 " AFNI_MODEL_DEBUG : specify debug/verbosity level\n"
325 "\n"
326 " e.g. setenv AFNI_MODEL_DEBUG 2\n"
327 "\n"
328 " Be more verbose. Valid levels are from 0 to 3, currently.\n"
329 "\n"
330 " AFNI_MODEL_DITER : specify debug iteration\n"
331 "\n"
332 " e.g. setenv AFNI_MODEL_DITER 999\n"
333 "\n"
334 " Get extra debug info at some iteration.\n"
335 "\n"
336 "----------------------------------------------------------------------\n"
337 " Written for E Silson and C Baker.\n"
338 "\n"
339 " R. Reynolds 27 June, 2014\n"
340 "----------------------------------------------------------------------\n"
341 );
342
343 return 0 ;
344 }
345
346