1 /*****************************************************************************
2 Based on model_conv_PRF.c, but this is with 6 parameters.
3
4 The original (4 parameter) model is from:
5
6 Population receptive field estimates in human visual cortex
7 NeuroImage 39 (2008) 647-660
8 Serge O. Dumoulin, Brian A. Wandell
9
10 Given stimulus images over time s(x,y,t), find x0, y0, sigma, R and theta
11 values that produce a best fit of the model to the data. Here x0, y0 are
12 taken to be the center of the population receptive field, sigma is the
13 basic width of it, R is the ratio of axis lengths (sigma_y / sigma_x), and
14 theta is the rotation from the y-direction major axis (so zero is in the
15 positive y-direction).
16
17 domains:
18 x,y : [-1,1], scaled by the mask, itself
19 sigma : (0,1], where 1 means the mask radius
20 R : [1,inf), since sigma defines the smaller size
21 theta : [-PI/2, PI/2), since rotation by PI has no effect
22
23 The model function of x0, y0, sigma, R and theta is constructed as follows:
24
25 1. generate a 2-D Gaussian density function, centered at x0, y0,
26 with given sigma, R (=sigma_y/sigma_x), and theta:
27
28 -> pRF model g(x,y) = e^-(A(x-x0)^2 + B(x-x0)(y-y0) + C(y-y0)^2)
29
30 where A = cos^2(theta)
31 ------------
32
33
34 2. integrate (dot product) over binary stimulus image per time point
35 -> pRF response r(t)
36 3. convolve with HRF
37 -> Prediction p(t)
38
39 Then perform a simple regression fit, returning A (amplitude), the input
40 x0, y0, sigma values, and the error (sum squares of the residuals)
41
42 inputs
43
44 s(x,y,t) : 2d+time visual field bitmap timeseries
45 (e.g. rotating stimulus wedge)
46 hrf : HRF convolution curve - assume GAM for now
47 data : time series, per voxel
48
49 outputs
50
51 A : amplitude of fit
52 x, y, sigma : best fit coordinates of center of Gaussian kernel
53
54 R. Reynolds June, 2014
55 ******************************************************************************/
56
57 #include "NLfit_model.h"
58
59 static int refnum = 0; /* # pts in refts */
60 static int refnz = 0; /* # of nonzero pts */
61 static float * refts = NULL; /* reference time series */
62 static int * refin = NULL; /* indexes of nonzero pts */
63 static int g_iter = -1; /* iteration number */
64
65 static char * g_model_ver = "model_conv_PRF_6_BAD, version 1.3, 21 Jun, 2018";
66
67 /* exp variables, parameters */
68 static float g_exp_maxval = 8.0; /* max x in exp(-x) */
69 static int g_exp_ipieces = 1000; /* vals per unit length */
70
71 /* exp variables, implied */
72 static int g_exp_nvals = 0; /* maxval*ipieces + 1 */
73 static float * g_exp_ts = NULL; /* exp(-x) for x at index VAL*pieces */
74
75 static THD_3dim_dataset * g_saset=NULL; /* stimulus aperture dataset */
76
77 /* prototypes */
78 static void conv_set_ref( int num , float * ref );
79 static int signal_model( float * , int , float ** , float *, int );
80 static int reset_stim_aperture_dset(int);
81 static int reset_exp_time_series(void);
82
83 static THD_3dim_dataset * convert_to_blurred_masks(THD_3dim_dataset *);
84 static THD_3dim_dataset * THD_reorg_dset(THD_3dim_dataset * din);
85 static int convolve_dset(THD_3dim_dataset * tset);
86 static float * get_float_volume_copy(THD_3dim_dataset * dset, int index,int nz);
87 static int compute_e_x_grid(float * e, int nx, int ny, float x0, float y0,
88 float sigma, float sigrat, float theta);
89 static int fill_computed_farray(float * ts, int tslen, THD_3dim_dataset * dset,
90 float x0, float y0, float sigma, float sigrat,
91 float theta, float A, int debug);
92 static int fill_scaled_farray(float * fdest, int nt, THD_3dim_dataset * dsrc,
93 float x, float y, float sigma, float scale, int debug);
94 static int get_ABC(float sigma, float sigrat, float theta,
95 double * A, double * B, double * C);
96 static int inputs_to_coords(THD_3dim_dataset * dset, float x, float y,
97 float sigma, float sigrat, float theta);
98
99 static int disp_floats(char * mesg, float * p, int len);
100 static int write_gauss_file(char * fname, float * curve,
101 int nx, int ny, char *hist);
102 static int model_help(void);
103 static int convolve_by_ref(float *, int, float *, int, int, int);
104
105 static void conv_model( float * gs , int ts_length ,
106 float ** x_array , float * ts_array );
107
108 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
109
110 /* MATH uses exp(x) directly, while ARRAY uses pre-computed values */
111 #define ACCUM_EXP_MATH(val) do { if(val>=0) result += exp(-val); } while (0)
112 #define ACCUM_EXP_ARRAY(val) do { \
113 eind = (int)(val*g_exp_ipieces); /* truncate towards 0? */ \
114 if ( eind >= 0 && eind <= g_exp_nvals ) \
115 result += g_exp_ts[eind]; } while (0)
116
117 /* ---------------------------------------------------------------------- */
118 /* interface to the environment */
119 static char * genv_conv_ref = NULL; /* AFNI_CONVMODEL_REF */
120 static char * genv_prf_stim = NULL; /* AFNI_MODEL_PRF_STIM_DSET */
121 static char * genv_gauss_file = NULL; /* AFNI_MODEL_PRF_GAUSS_FILE */
122 static int genv_diter = -1; /* debug iteration */
123 static int genv_debug = 0; /* AFNI_MODEL_DEBUG */
124
125 static int genv_on_grid = 0; /* restrict computations to grid */
126 static float genv_sigma_max = 1.0; /* on_grid: maximum blur sigma */
127 static int genv_sigma_nsteps = 100; /* on_grid: number of blur steps */
128 static int genv_sigma_ratio_nsteps = 4; /* integers >= 1, or any >= 1 if 0 */
129 static int genv_theta_nsteps = 6; /* truncate [-PI/2,PI/2) to S steps
130 (if > 0) */
131
132 static int genv_get_help = 0; /* AFNI_MODEL_HELP_ALL/HELP_CONV_PRF_6 */
133
set_env_vars(void)134 static int set_env_vars(void)
135 {
136 genv_conv_ref = my_getenv("AFNI_CONVMODEL_REF"); /* reference file */
137 if( genv_conv_ref ) fprintf(stderr,"-- PRF: have REF %s\n", genv_conv_ref);
138 else fprintf(stderr,"** model PRF: AFNI_CONVMODEL_REF is not set\n");
139
140 genv_prf_stim = my_getenv("AFNI_MODEL_PRF_STIM_DSET"); /* visual stim */
141 if( genv_prf_stim ) fprintf(stderr,"-- PRF: have stim %s\n", genv_prf_stim);
142 else fprintf(stderr,"** model PRF: AFNI_MODEL_PRF_STIM_DSET is not set\n");
143
144 genv_diter = (int)AFNI_numenv_def("AFNI_MODEL_DITER", -1);
145 genv_debug = (int)AFNI_numenv_def("AFNI_MODEL_DEBUG", 0);
146 fprintf(stderr,"-- PRF: debug %d, iter %d\n", genv_debug, genv_diter);
147
148 /* on grid - default to no (yes->no 6 Jun 2013) */
149 genv_on_grid = AFNI_yesenv("AFNI_MODEL_PRF_ON_GRID"); /* flag */
150 fprintf(stderr,"-- PRF: results on grid: %s\n", genv_on_grid?"yes":"no");
151
152 genv_sigma_max = AFNI_numenv_def("AFNI_MODEL_PRF_SIGMA_MAX", genv_sigma_max);
153 genv_sigma_nsteps = (int)AFNI_numenv_def("AFNI_MODEL_PRF_SIGMA_NSTEPS",
154 genv_sigma_nsteps);
155 genv_sigma_ratio_nsteps =
156 (int)AFNI_numenv_def("AFNI_MODEL_PRF_SIGMA_RATIO_NSTEPS",
157 genv_sigma_ratio_nsteps);
158 genv_theta_nsteps = (int)AFNI_numenv_def("AFNI_MODEL_PRF_THETA_NSTEPS",
159 genv_theta_nsteps);
160 if( genv_on_grid )
161 fprintf(stderr, "-- PRF: sigma_max = %f, nsteps = %d,"
162 " ratio_steps = %d, theta_steps = %d\n",
163 genv_sigma_max, genv_sigma_nsteps,
164 genv_sigma_ratio_nsteps, genv_theta_nsteps);
165
166 /* help */
167 genv_get_help = AFNI_yesenv("AFNI_MODEL_HELP_CONV_PRF_6_BAD")
168 || AFNI_yesenv("AFNI_MODEL_HELP_ALL");
169
170 /* write a Gaussian mask? */
171 genv_gauss_file = my_getenv("AFNI_MODEL_PRF_GAUSS_FILE");
172
173 return 0;
174 }
175 /* ---------------------------------------------------------------------- */
176
177
178 /*----------------------------------------------------------------------
179 Function to set the reference time series, with which the
180 model function is convolved to produce the simulated data.
181 ------------------------------------------------------------------------*/
182
conv_set_ref(int num,float * ref)183 static void conv_set_ref( int num , float * ref )
184 {
185 if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
186 int ii ;
187
188 /* get rid of old data */
189
190 if(refts != NULL){ free(refts); refts = NULL; free(refin); refin = NULL; }
191
192 refnum = num ;
193 refts = (float *) malloc( sizeof(float) * num ) ;
194 refin = (int *) malloc( sizeof(int) * num ) ;
195 memcpy( refts , ref , sizeof(float) * num ) ;
196 for( ii=0,refnz=0 ; ii < num ; ii++ ) /* build list of nonzero */
197 if( refts[ii] != 0.0 ) refin[refnz++] = ii ; /* points in refts */
198 if( refnz == 0 )
199 ERREX("model_conv_PRF: All zero reference timeseries!") ;
200
201 if( genv_debug ) {
202 fprintf(stderr,"+d conv_set_ref: num=%d nonzero=%d\n",num,refnz) ;
203 if( genv_debug > 1 ) {
204 fprintf(stderr," TR locked stimuli :");
205 for( ii = 0; ii < refnz; ii++ )
206 fprintf(stderr," %f", refts[refin[ii]]);
207 fputc('\n',stderr);
208 }
209 }
210
211 return;
212
213 } else { /*** if no inputs, read it from AFNI_CONVMODEL_REF 1D file ***/
214
215 MRI_IMAGE * flim ;
216
217 if( genv_conv_ref == NULL )
218 ERREX("model_conv_PRF: need ref file as AFNI_CONVMODEL_REF") ;
219
220 flim = mri_read_1D(genv_conv_ref);
221 if( flim == NULL ){
222 fprintf(stderr,"** model_conv_PRF: Can't read timeseries file %s",
223 genv_conv_ref);
224 ERREX("failing...");
225 }
226
227 if( genv_debug ) fprintf(stderr,"+d conv_set_ref: refts=%s nx=%d\n",
228 genv_conv_ref, flim->ny);
229
230 conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ; /* recursion! */
231 mri_free(flim) ;
232 }
233 return;
234 }
235
236 /* any failure should leave g_saset == NULL
237 *
238 * return 0 on success */
reset_stim_aperture_dset(int needed_length)239 static int reset_stim_aperture_dset(int needed_length)
240 {
241 THD_3dim_dataset * sanew;
242 int errs=0;
243
244 /* free and reload saset */
245 if( g_saset ) DSET_delete(g_saset);
246 g_saset = THD_open_dataset(genv_prf_stim);
247 if( ! g_saset ) return 1;
248
249 /* check for square dataset and sufficient nt */
250 if( fabs(DSET_DX(g_saset) - DSET_DY(g_saset)) > 0.0001 ) {
251 fprintf(stderr,"** PRF: stimset voxels not square (%f, %f)\n",
252 DSET_DX(g_saset), DSET_DY(g_saset));
253 errs++;
254 }
255
256 if( DSET_NX(g_saset) != DSET_NY(g_saset) ) {
257 fprintf(stderr,"** PRF: stimset not square (%d, %d)\n",
258 DSET_NX(g_saset), DSET_NY(g_saset));
259 errs++;
260 }
261
262 if( DSET_NVALS(g_saset) < needed_length ) {
263 fprintf(stderr,"** PRF: dset nt = %d, tslen = %d\n",
264 DSET_NVALS(g_saset), needed_length);
265 errs++;
266 }
267
268 if( errs ) { DSET_delete(g_saset); g_saset = NULL; return 1; }
269
270
271 if( THD_dset_to_mask(g_saset, 1.0, 0.0) ) return 1;
272
273 if( genv_on_grid ) {
274 sanew = convert_to_blurred_masks(g_saset);
275 DSET_delete(g_saset);
276 g_saset = sanew;
277 //EDIT_dset_items(sanew, ADN_prefix, "sa.blur", ADN_none);
278 //fprintf(stderr, "++ writing blur dset...\n");
279 //DSET_write(sanew);
280
281 sanew = THD_reorg_dset(g_saset);
282 DSET_delete(g_saset);
283 g_saset = sanew;
284 //EDIT_dset_items(sanew, ADN_prefix, "sa.reorg", ADN_none);
285 //fprintf(stderr, "++ writing reorg dset...\n");
286 //DSET_write(sanew);
287
288 convolve_dset(g_saset);
289 //EDIT_dset_items(sanew, ADN_prefix, "sa.conv", ADN_none);
290 //fprintf(stderr, "++ writing convolve dset...\n");
291 //DSET_write(sanew);
292
293 if( ! g_saset ) return 1;
294 }
295
296 return 0;
297 }
298
299
300 /* x-axis is convolution time axis, so for each row, convolve */
convolve_dset(THD_3dim_dataset * tset)301 static int convolve_dset(THD_3dim_dataset * tset)
302 {
303 float * result, * volbase, * dbase;
304 int nx, ny, nz, nt, nxy;
305 int xind, yind, zind, tind;
306
307 if( ! tset ) { fprintf(stderr, "** no tset to convolve\n"); return 1; }
308
309 if(genv_debug>1)fprintf(stderr, "++ starting convolution at time %6.1f\n",
310 0.001*NI_clock_time());
311
312 /* allocate room for computation and result */
313 nx = DSET_NX(tset); ny = DSET_NY(tset); nz = DSET_NZ(tset);
314 nt = DSET_NVALS(tset); nxy = nx * ny;
315 result = (float *)malloc(nx * sizeof(float));
316 if( !result ) {
317 fprintf(stderr, "** PRF: conv_dset: failed malloc of %d floats\n", nx);
318 return 1;
319 }
320
321 for( tind = 0; tind < nt; tind++ ) {
322 volbase = (float *)DSET_ARRAY(tset, tind);
323
324 /* for each row (the signal), convolve */
325 for( zind = 0; zind < nz; zind++ ) {
326 for( yind = 0; yind < ny; yind++ ) {
327 /* result = convolved signal (initialize and demean) */
328 dbase = volbase + yind*nx + zind*nxy;
329 convolve_by_ref(result, nx, dbase, nx, 1, 1);
330 /* and put result back in */
331 for( xind = 0; xind < nx; xind++ ) dbase[xind] = result[xind];
332 }
333 }
334 }
335
336 if(genv_debug>1)fprintf(stderr, "-- finished convolution at time %6.1f\n",
337 0.001*NI_clock_time());
338
339 return 0;
340 }
341
342
343 /* Convert 2D+time byte binary mask dataset into 3D+time float
344 * dataset of slices (blur_nsteps) at different blur levels
345 * (multiples of blur_res).
346 *
347 * input: 1 binary mask slice per TR
348 * output: duplicate dset, but where each slice is converted to an array of
349 * slices at different blur levels
350 *
351 * Resample (reorganize) output so that time is the fastest dimension,
352 * then x, y and sigma (sub-brick direction).
353 *
354 * Convert to float array for now, too, just to waste RAM.
355 * At least that will make blurring more straightforward
356 * and avoid truncation issues.
357 *
358 * given blur_max (e.g. 0.1)
359 * blur_nsteps (e.g. 50)
360 * compute blur_res (e.g. 0.002)
361 *
362 *------------------------------------------------------------
363 * for each TR
364 * convert mask slice to float
365 * allocate for volume of floats
366 * duplicate slice across volume slices
367 * for each blur level
368 * EDIT_blur_volume(NX, NY, NZ=1, dx, dy, dz, MRI_float, data, SIGMA);
369 * attach to dset
370 * EDIT_dset_item: set datum, NZ
371 *------------------------------------------------------------
372 */
convert_to_blurred_masks(THD_3dim_dataset * dset)373 static THD_3dim_dataset * convert_to_blurred_masks(THD_3dim_dataset * dset)
374 {
375 THD_3dim_dataset * dnew;
376 THD_ivec3 iv_nxyz;
377 float * fdata, * foffset, sigma;
378 int nx, ny, nt;
379 int vind, sind;
380
381 /* test required inputs */
382 if( !dset ) return NULL;
383 if( genv_sigma_max <= 0.0 || genv_sigma_max > 1.0 || genv_sigma_nsteps<=1 ){
384 fprintf(stderr,"** PRF on grid: bad sigma max %f or nsteps %d\n",
385 genv_sigma_max, genv_sigma_nsteps);
386 return NULL;
387 }
388
389 if( DSET_NZ(dset) != 1 ) {
390 fprintf(stderr,"** invalid stim NZ = %d\n", DSET_NZ(dset));
391 return NULL;
392 }
393
394 /* create initial copy */
395 dnew = EDIT_empty_copy(dset);
396 if( !dnew ) return NULL;
397
398 nx = DSET_NX(dnew); ny = DSET_NY(dnew); nt = DSET_NVALS(dnew);
399 LOAD_IVEC3(iv_nxyz , nx, ny, genv_sigma_nsteps);
400
401 EDIT_dset_items(dnew,
402 ADN_prefix, "I_like_jello",
403 ADN_nxyz, iv_nxyz,
404 ADN_datum_all, MRI_float,
405 ADN_none);
406
407 if(genv_debug)
408 fprintf(stderr, "++ making blurred time series: %d x %d x %d x nt=%d\n",
409 DSET_NX(dnew), DSET_NY(dnew), DSET_NZ(dnew), DSET_NVALS(dnew));
410 if(genv_debug>1)fprintf(stderr, "++ starting blur at time %6.1f\n",
411 0.001*NI_clock_time());
412
413 for( vind = 0; vind < nt; vind++ ) {
414 if( genv_debug > 1 ) fputc('.', stderr);
415
416 fdata = get_float_volume_copy(dset, vind, genv_sigma_nsteps);
417 if( !fdata ) { DSET_delete(dnew); return NULL; }
418
419 for( sind = 0; sind < genv_sigma_nsteps; sind++ ) {
420 sigma = genv_sigma_max * ((sind + 1.0)/genv_sigma_nsteps);
421 foffset = fdata + sind*nx*ny;
422 FIR_blur_volume_3d(nx, ny, 1, 2.0/(nx-1), 2.0/(ny-1), 1.0,
423 foffset, sigma, sigma, 0.0);
424 }
425
426 mri_fix_data_pointer(fdata, DSET_BRICK(dnew, vind));
427 }
428
429 if(genv_debug>1)fprintf(stderr, "\n-- finished blur volume at time %6.1f\n",
430 0.001*NI_clock_time());
431
432 return dnew;
433 }
434
435 /* rotate list of axes and correspondingly shuffle data:
436 * x,y,z,t -> t,x,y,z
437 *
438 * The main point is to make the time axes the fast direction.
439 *
440 * require MRI_float for now */
THD_reorg_dset(THD_3dim_dataset * din)441 static THD_3dim_dataset * THD_reorg_dset(THD_3dim_dataset * din)
442 {
443 THD_3dim_dataset * dout;
444 THD_ivec3 iv_nxyz;
445 float * newvol, * inslice, * outbase, * inbase;
446 int in_nx, in_ny, in_nz, in_nt;
447 int out_nx, out_ny, out_nz, out_nt;
448 int in_nxyz, in_nxy, out_nxyz, out_nxy;
449 int xind, yind, zind, tind;
450
451 if( !din ) { fprintf(stderr,"** reorg_dset: missing input\n"); return NULL; }
452
453 dout = EDIT_empty_copy(din);
454
455 in_nx = DSET_NX(din); in_ny = DSET_NY(din); in_nz = DSET_NZ(din);
456 in_nt = DSET_NVALS(din);
457 in_nxy = in_nx * in_ny;
458 in_nxyz = in_nx * in_ny * in_nz;
459
460 out_nx = in_nt; out_ny = in_nx; out_nz = in_ny; out_nt = in_nz;
461 out_nxy = out_nx * out_ny;
462 out_nxyz = out_nx * out_ny * out_nz;
463
464 LOAD_IVEC3(iv_nxyz , out_nx, out_ny, out_nz);
465
466 EDIT_dset_items(dout,
467 ADN_nxyz, iv_nxyz, ADN_nvals, out_nt, ADN_ntt, out_nt,
468 ADN_none);
469
470 if(genv_debug)
471 fprintf(stderr,"-- reorg_dset: nxyzt (%d,%d,%d,%d) -> (%d,%d,%d,%d)\n",
472 in_nx, in_ny, in_nz, in_nt, out_nx, out_ny, out_nz, out_nt);
473 if(genv_debug>1)fprintf(stderr, "\n== reorg starting at %6.1f\n",
474 0.001*NI_clock_time());
475
476 /* create and attach each volume */
477 for( tind = 0; tind < out_nt; tind++ ) {
478 newvol = (float *)malloc(out_nxyz * sizeof(float));
479 if( !newvol ) {
480 fprintf(stderr,"** PRF reorg: failed to alloc volume %d\n", tind);
481 DSET_delete(dout);
482 return NULL;
483 }
484
485 for( xind = 0; xind < out_nx; xind++ ) {
486
487 /* copy one slice of input across y and z directions */
488 inslice = ((float *)DSET_ARRAY(din, xind)) + tind*in_nxy;
489
490 for( yind = 0; yind < out_ny; yind++ ) {
491 outbase = newvol + xind + yind*out_nx;
492 inbase = inslice + yind;
493 for( zind = 0; zind < out_nz; zind++ )
494 outbase[zind*out_nxy] = inbase[zind*in_nx];
495 }
496
497 }
498
499 mri_fix_data_pointer(newvol, DSET_BRICK(dout, tind));
500 }
501
502 if(genv_debug>1)fprintf(stderr, "\n== reorg finished at %6.1f\n",
503 0.001*NI_clock_time());
504
505 if( genv_debug > 2 ) {
506 MRI_IMAGE * im;
507 float * fp;
508 int i, j, k;
509
510 i = in_nx/3; j = in_ny/3; k = in_nz/3;
511
512 im = THD_extract_series(i+j*in_nx+k*in_nxy, din, 0);
513 disp_floats("== ARY: sig [nxyz/3]: ", MRI_FLOAT_PTR(im), in_nt);
514 fp = ((float *)DSET_ARRAY(dout,k)) + i*out_nx + j*out_nxy;
515 disp_floats("== ARY: reorg : ", fp, out_nx);
516
517 mri_free(im);
518 }
519
520 return dout;
521 }
522
523
524 /* duplicate dset slice at given index nz times (into volume of nz slices) */
get_float_volume_copy(THD_3dim_dataset * dset,int index,int nz)525 static float * get_float_volume_copy(THD_3dim_dataset * dset, int index, int nz)
526 {
527 MRI_IMAGE * fim;
528 float * fdata;
529 int nxy, sind;
530
531 /* make space for new data */
532 nxy = DSET_NXY(dset);
533 fdata = (float *)malloc(nxy * nz * sizeof(float));
534 if( ! fdata ) {
535 fprintf(stderr,"** PRF: failed to alloc %d floats at index %d\n",
536 nxy * nz, index);
537 return NULL;
538 }
539
540 /* get slice to duplicate */
541 fim = mri_to_float( DSET_BRICK(dset, index) );
542 if( ! fim ) {
543 fprintf(stderr,"** PRF: failed mri_to_float\n");
544 free(fdata);
545 return NULL;
546 }
547
548 /* dupe over slices */
549 for( sind = 0; sind < nz; sind++ )
550 memcpy(fdata+sind*nxy, MRI_FLOAT_PTR(fim), nxy * sizeof(float));
551
552 mri_free(fim);
553
554 return fdata;
555 }
556
557
558 /* want e^-x computed out to x=7 (c^-7 goes below 0.001)
559 *
560 * g_exp_maxval = 7
561 * g_exp_ipieces = 1000
562 * so step = 1/1000
563 * g_exp_nvals = g_exp_maxval * g_exp_ipieces + 1 (for 0), plus a few
564 */
reset_exp_time_series(void)565 static int reset_exp_time_series(void)
566 {
567 int ind;
568 float resol = 1.0/g_exp_ipieces;
569
570 g_exp_nvals = (int)(g_exp_maxval * g_exp_ipieces) + 1;
571
572 if(genv_debug) fprintf(stderr, "-- exp nvals = %d, max = %f, pieces = %d\n",
573 g_exp_nvals, g_exp_maxval, g_exp_ipieces);
574
575 if( g_exp_ts ) free(g_exp_ts);
576 g_exp_ts = (float *)malloc(g_exp_nvals * sizeof(float));
577 if( ! g_exp_ts ) {
578 fprintf(stderr,"** failed to alloc %d floats, buring...\n", g_exp_nvals);
579
580 if( g_saset ) DSET_delete(g_saset);
581 g_saset = NULL; /* this blocks analysis */
582 return 1;
583 }
584
585 for( ind = 0; ind < g_exp_nvals; ind++ )
586 g_exp_ts[ind] = exp(-resol*ind);
587
588 return 0;
589 }
590
591 /*-----------------------------------------------------------------------
592 Function to compute the simulated time series.
593 -------------------------------------------------------------------------*/
594
conv_model(float * gs,int ts_length,float ** x_array,float * ts_array)595 static void conv_model( float * gs , int ts_length ,
596 float ** x_array , float * ts_array )
597 {
598 int ii, cur_debug = 0, irfdur=0;
599 int iter_step;
600
601 static int nid = 0 ; /* number of pts in impulse */
602 static float * fid = NULL; /* impulse response function */
603
604 /*----- check for env vars -----*/
605
606 g_iter++ ;
607
608 if( g_iter == 0 ) {
609 set_env_vars(); /* process environment variables */
610 if(genv_debug && x_array) fprintf(stderr,"\n+d TR = %f\n",
611 x_array[1][1]-x_array[0][1]);
612 fprintf(stderr,"++ %s\n", g_model_ver);
613 }
614
615 /*** make sure there is a reference function to convolve with ***/
616 /* it may be used in reset_stim_aperture_dset */
617
618 if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
619
620 if ( genv_on_grid ) {
621 int pad = (g_iter % 76) - 38;
622 if ( pad < 0 ) pad = -pad;
623 fprintf(stderr,"%*s** not ready for AFNI_MODEL_PRF_ON_GRID **\n",
624 pad,"");
625 for ( ii=0; ii < ts_length; ii++ )
626 ts_array[ii] = 0;
627 return;
628 }
629
630 /* create stim aperture dset */
631 if( g_iter == 0 ) {
632 (void)reset_stim_aperture_dset(ts_length); /* free and reload saset */
633 (void)reset_exp_time_series(); /* pre-compute exp(x) */
634 if( genv_debug ) fprintf(stderr, "== start time %d\n", NI_clock_time());
635 }
636
637 /*** initialize the output before possible early return ***/
638
639 for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
640
641 /* if we had some failure, bail */
642 if( !g_saset ) return;
643
644 if( genv_debug > 1 ) { /* babble */
645 if( genv_on_grid ) iter_step = 100000;
646 else iter_step = 1000;
647 if( (g_iter % iter_step) == 0 ) {
648 if( g_iter % (10*iter_step) ) fputc('\r', stderr);
649 fprintf(stderr, "-- time for %d iter set %5d : %6.1f\n",
650 iter_step, g_iter/iter_step, 0.001*NI_clock_time());
651 }
652 }
653
654 /* to clean up, particularly as a parameter */
655 cur_debug = (g_iter == genv_diter || (g_iter == 0 && genv_debug > 1));
656 if( cur_debug ) disp_floats("+d input params: ", gs, 4);
657
658 /*** initialize the impulse response ***/
659
660 if( nid < ts_length ){ /* make some space for it */
661 if( fid ) free(fid) ;
662 nid = ts_length ;
663 fid = (float *) malloc( sizeof(float) * nid ) ;
664 }
665
666 /* If on grid, get signal directly. Else, convolve signal impulse
667 function with ref time series, and demean. */
668 if( genv_on_grid ) {
669 irfdur = signal_model(gs, ts_length, x_array, ts_array, cur_debug);
670 } else {
671 irfdur = signal_model(gs, ts_length, x_array, fid, cur_debug);
672 convolve_by_ref(ts_array, ts_length, fid, irfdur, 0, 1);
673 if( cur_debug ) disp_floats("+d no_grid conv : ", ts_array, ts_length);
674 }
675
676 return ;
677 }
678
679 /*
680 Convolve reference function (GAM here?) with stim time series.
681
682 globals:
683 int refnum : # pts in refts
684 int refnz : # of nonzero pts
685 float * refts : reference time series
686 int * refin : indexes of nonzero pts
687
688 inputs:
689 result : convolved result (already initialized to 0)
690 rlen : length of input data time series (maybe equal)
691 signal : signal model output
692 siglen : length of signal time series
693 init : flag: init result?
694 demean : flag: demean result?
695 */
convolve_by_ref(float * result,int rlen,float * signal,int siglen,int init,int demean)696 static int convolve_by_ref(float * result, int rlen, float * signal, int siglen,
697 int init, int demean)
698 {
699 int ii, jj, kk, jtop;
700 float val; /* current reference value to apply */
701 double mean; /* for demeaning result */
702
703 if( init ) {
704 for( ii=0 ; ii < rlen ; ii++ ) result[ii] = 0.0;
705 }
706
707 /* TR-locked convolution */
708 for( ii=0 ; ii < refnz ; ii++ ){
709 kk = refin[ii] ; if( kk >= rlen ) break ;
710 val = refts[kk] ;
711
712 /* for each point in the impulse, add its val times irf */
713 /* (top index offset is min(irfdur, length-kk-1)) */
714 jtop = rlen - kk ; if( jtop > siglen ) jtop = siglen ;
715 for( jj=0 ; jj < jtop ; jj++ )
716 result[kk+jj] += val * signal[jj];
717 }
718
719 if( ! demean ) return 0;
720
721 /* demean the result */
722 mean = 0.0;
723 for( ii=0 ; ii < rlen ; ii++ ) mean += result[ii];
724 mean /= rlen;
725 for( ii=0 ; ii < rlen ; ii++ ) result[ii] -= mean;
726
727 return 0;
728 }
729
disp_floats(char * mesg,float * p,int len)730 static int disp_floats(char * mesg, float * p, int len)
731 {
732 int c;
733 if( mesg ) fputs(mesg, stderr);
734 for( c = 0; c < len; c++ ) fprintf(stderr," %f ", p[c]);
735 fprintf(stderr,"\n\n");
736 return 0;
737 }
738
739
740 /*-----------------------------------------------------------------------*/
741
742 /*----------------------------------------------------------------------*
743 * for use in signal model
744 *----------------------------------------------------------------------*/
745 #ifdef PI
746 #undef PI
747 #endif
748 #define PI 3.141592653589793238462643
749
750
751 DEFINE_MODEL_PROTOTYPE
752
initialize_model()753 MODEL_interface * initialize_model ()
754 {
755 MODEL_interface * mi = NULL;
756
757 /*----- first, see if the user wants help -----*/
758 if ( AFNI_yesenv("AFNI_MODEL_HELP_CONV_PRF_6_BAD") ||
759 AFNI_yesenv("AFNI_MODEL_HELP_ALL") ) model_help();
760
761 /*----- allocate memory space for model interface -----*/
762
763 mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));
764
765 /*----- name of this model -----*/
766
767 strcpy (mi->label, "Conv_PRF_6_BAD");
768
769 /*----- this is a signal model -----*/
770
771 mi->model_type = MODEL_SIGNAL_TYPE;
772
773 /*----- number of parameters in the model -----*/
774
775 mi->params = 6;
776
777 /*----- parameter labels -----*/
778
779 strcpy (mi->plabel[0], "Amp");
780 strcpy (mi->plabel[1], "X");
781 strcpy (mi->plabel[2], "Y");
782 strcpy (mi->plabel[3], "Sigma");
783 strcpy (mi->plabel[4], "SigRat");
784 strcpy (mi->plabel[5], "Theta");
785
786 /*----- minimum and maximum parameter constraints -----*/
787
788 /* amplitude, x/y ranges, and sigma range */
789 mi->min_constr[0] = -10.0; mi->max_constr[0] = 10.0;
790
791 mi->min_constr[1] = -1.0; mi->max_constr[1] = 1.0;
792 mi->min_constr[2] = -1.0; mi->max_constr[2] = 1.0;
793
794 mi->min_constr[3] = 0.0; mi->max_constr[3] = 1.0;
795 mi->min_constr[4] = 1.0; mi->max_constr[4] = 10.0;
796 mi->min_constr[5] = -PI/2.0; mi->max_constr[5] = PI/2.0;
797
798 /*----- function which implements the model -----*/
799 mi->call_func = conv_model;
800
801 return (mi);
802 }
803
804
805 /*----------------------------------------------------------------------*/
806 /*
807 Routine to calculate the time series to (hopefully) fit the data.
808
809 Definition of model parameters (gs[2] > 0)
810
811 gs[0] = Amp = amplitude
812 gs[1] = x0 = x-coordinate of gaussian center
813 gs[2] = y0 = y-coordinate of gaussian center
814 gs[3] = sigma = "width" of gaussian curve
815 gs[4] = sigrat = sigma ratio = sigma_y / sigma_x
816 gs[5] = theta = angle from "due north"
817
818 For each TR, integrate g(x,y) over stim aperture dset.
819
820 g(x,y) = e^-[((x-x0)^2+(y-y0)^2)/(2*sigma^2)]
821
822 The resulting returned time series will be convolved in the
823 parent function.
824 */
signal_model(float * gs,int ts_length,float ** x_array,float * ts_array,int debug)825 static int signal_model
826 (
827 float * gs, /* parameters for signal model */
828 int ts_length, /* length of time series data */
829 float ** x_array, /* independent variable matrix */
830 float * ts_array, /* estimated signal model time series */
831 int debug /* make some noise */
832 )
833 {
834 int maxind, tmpmax;/* largest dimension */
835 float A, x, y, sigma;/* model params */
836 float sigrat, theta;
837
838 /* assign parameters */
839 A = gs[0];
840 x = gs[1]; y = gs[2];
841 sigma = gs[3]; sigrat = gs[4]; theta = gs[5];
842
843 if( debug ) fprintf(stderr, "-d model_conv_PRF parameters: "
844 "A = %f, x = %f, y = %f, sigma = %f\n"
845 " sigrat = %f, theta = %f\n"
846 " nz = %d, nvals = %d, ts_len = %d\n",
847 A, x, y, sigma, sigrat, theta,
848 DSET_NZ(g_saset), DSET_NVALS(g_saset), ts_length);
849
850 if( ! ISVALID_3DIM_DATASET(g_saset) ) return 0;
851
852 maxind = ts_length;
853 if( genv_on_grid ) tmpmax = DSET_NX(g_saset);
854 else tmpmax = DSET_NVALS(g_saset);
855
856 if( maxind > tmpmax ) maxind = tmpmax;
857 if( maxind == 0 ) return 0;
858
859 if( debug )
860 fprintf( stderr,"-d NT orig=%d, applied=%d\n", ts_length, maxind);
861
862 /* time array must be ordered according to stim dset */
863 if( genv_on_grid ) /* scale data directly from grid */
864 {
865 fprintf(stderr,"== rcr - need to apply sigrat, theta on grid\n");
866 fill_scaled_farray(ts_array, maxind, g_saset, x, y, sigma, A, debug);
867 }
868 else
869 fill_computed_farray(ts_array, maxind, g_saset, x, y,
870 sigma, sigrat, theta, A, debug);
871
872 if( debug )
873 disp_floats("+d signal model result : ", ts_array, ts_length);
874
875 return maxind;
876 }
877
878 /* get j, k, t from x, y, sigma, then copy data */
fill_scaled_farray(float * fdest,int length,THD_3dim_dataset * dsrc,float x,float y,float sigma,float scale,int debug)879 static int fill_scaled_farray(float *fdest, int length, THD_3dim_dataset *dsrc,
880 float x, float y, float sigma, float scale, int debug)
881 {
882 float * inptr, * outptr;
883 int nx, ny, nz, nt;
884 int i, j, k, t;
885
886 nx = DSET_NX(dsrc); ny = DSET_NY(dsrc); nz = DSET_NZ(dsrc);
887 nt = DSET_NVALS(dsrc);
888
889 if( nx != length ) {
890 fprintf(stderr, "** FSF: nx, len mis-match, %d != %d\n", nx, length);
891 return 1;
892 }
893
894 /* for j,k, map [-1,1] to [0, ny-1] */
895 j = (int)(0.5+ny*(x+1.0)/2.0); /* x-coord -> j index */
896 if (j < 0 ) j = 0;
897 else if(j >= ny ) j = ny-1;
898
899 k = (int)(0.5+nz*(y+1.0)/2.0); /* y-coord -> k index */
900 if (k < 0 ) k = 0;
901 else if(k >= nz ) k = nz-1;
902
903 /* init to round(nsteps * fraction of max) */
904 t = (int)(0.5 + genv_sigma_nsteps * sigma / genv_sigma_max);
905 if ( t < 0 ) t = 0;
906 else if( t >= nt ) t = nt-1;
907
908 if( debug )
909 fprintf(stderr,"-- fill_array from x=%f, y=%f, s=%f\n"
910 " at j=%d, k=%d, t=%d\n",
911 x, y, sigma, j, k, t);
912
913 inptr = ((float *)DSET_ARRAY(dsrc, t)) + j*nx + k*nx*ny;
914 outptr = fdest;
915 for( i = 0; i < length; i++ ){
916 *outptr = scale**inptr;
917 inptr++;
918 outptr++;
919 }
920
921 return 0;
922 }
923
924
925 /* as in fill_scaled_farray, but map x,y,s to i,j,k */
inputs_to_coords(THD_3dim_dataset * dset,float x,float y,float sigma,float sigrat,float theta)926 static int inputs_to_coords(THD_3dim_dataset * dset, float x, float y,
927 float sigma, float sigrat, float theta)
928 {
929 int nx, ny, nz;
930 int i, j, k;
931 float eval;
932 double A, B, C;
933
934 nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset);
935
936 /* for i,j, map [-1,1] to [0, nx-1] */
937 i = (int)(0.5+nx*(x+1.0)/2.0);
938 if (i < 0 ) i = 0;
939 else if(i >= nx ) i = nx-1;
940
941 j = (int)(0.5+ny*(y+1.0)/2.0);
942 if (j < 0 ) j = 0;
943 else if(j >= ny ) j = ny-1;
944
945 /* init to round(nsteps * fraction of max) */
946 k = (int)(0.5 + genv_sigma_nsteps * sigma / genv_sigma_max);
947 if ( k < 0 ) k = 0;
948 else if( k >= nz ) k = nz-1;
949
950 get_ABC(sigma, sigrat, theta, &A, &B, &C); /* get exp coefficients */
951 eval = A*x*x + B*x*y + C*y*y;
952
953 fprintf(stderr,"-- fill_array from x=%f, y=%f, s=%f\n"
954 " at i=%d, j=%d, k=%d\n",
955 x, y, sigma, i, j, k);
956 fprintf(stderr," sigrat=%f, theta=%f, exp=%.3f\n",
957 sigrat, theta, eval);
958
959 return 0;
960 }
961
962
963
964 /* compute the response value directly, given parameters:
965 * - compute at a given index, under the assumption that
966 * such computations will be sequential (at least starting from 0)
967 * - at index 0, a new e^x slice will be computed and
968 * applied across time */
fill_computed_farray(float * ts,int tslen,THD_3dim_dataset * dset,float x0,float y0,float sigma,float sigrat,float theta,float A,int debug)969 static int fill_computed_farray(float * ts, int tslen, THD_3dim_dataset * dset,
970 float x0, float y0, float sigma, float sigrat,
971 float theta, float A, int debug)
972 {
973 static float * sexpgrid = NULL;
974 static int snxy=0;
975 float * eptr;
976 byte * mptr;
977 double sum;
978 int nx, ny, tind, sind, nmask;
979
980 nx = DSET_NX(dset); ny = DSET_NY(dset);
981
982 /* maybe allocate new memory for egrid (hopefully once, "ever") */
983 if( ! sexpgrid || snxy != nx*ny ) {
984 if( genv_debug )
985 fprintf(stderr,"++ alloc egrid, snxy = %d, nxy = %d\n", snxy, nx*ny);
986 snxy = nx*ny;
987 if( sexpgrid ) free(sexpgrid); /* nuke any old copy - maybe never */
988 sexpgrid = (float *)calloc(snxy, sizeof(float *));
989 if( !sexpgrid ) {
990 fprintf(stderr,"** PRF egrid alloc failure, nxy = %d\n", snxy);
991 return 1;
992 }
993 }
994
995 /* get current e^x grid (scaled to unit area) */
996 if( compute_e_x_grid(sexpgrid, nx, ny, x0, y0, sigma, sigrat, theta) ) {
997 fprintf(stderr,"PRF:FCA: e_x_g failure\n");
998 return 1;
999 }
1000
1001 for( tind = 0; tind < tslen; tind++ ) {
1002 mptr = DBLK_ARRAY(dset->dblk, tind);
1003 eptr = sexpgrid;
1004 sum = 0.0;
1005 nmask = 0;
1006 for( sind = 0; sind < nx*ny; sind++ ) {
1007 if( *mptr++ ) {
1008 nmask++;
1009 sum += *eptr;
1010 }
1011 eptr++;
1012 }
1013
1014 if( debug && genv_debug > 2 )
1015 fprintf(stderr,"-- nmask %03d = %d\n",tind,nmask);
1016
1017 ts[tind] = A * sum;
1018 }
1019
1020 /* if requested, write this 2D image (one time only) */
1021 if( genv_gauss_file ) {
1022 char hist[256];
1023 sprintf(hist, "\n == %s\n x = %g, y = %g, "
1024 "sigma = %g, sigrat = %g, theta = %g\n",
1025 g_model_ver, x0, y0, sigma, sigrat, theta);
1026 fprintf(stderr, "++ writing PRF model curve to %s\n%s\n",
1027 genv_gauss_file, hist);
1028 write_gauss_file(genv_gauss_file, sexpgrid, nx, ny, hist);
1029
1030 genv_gauss_file = NULL; /* clear - no further writes */
1031 }
1032
1033 return 0;
1034 }
1035
1036 /* ------------------------------------------------------------ */
1037 /* write_gauss_curve */
write_gauss_file(char * fname,float * curve,int nx,int ny,char * hist)1038 static int write_gauss_file(char * fname, float * curve, int nx, int ny,
1039 char * hist)
1040 {
1041 THD_3dim_dataset * dout;
1042 THD_ivec3 inxyz;
1043 THD_fvec3 origin, delta;
1044 float * mptr, * dptr;
1045 byte * bptr;
1046 int ind;
1047
1048 fprintf(stderr,"++ creating gauss dset (%dx%d)", nx, ny);
1049 dout = EDIT_empty_copy(NULL);
1050 LOAD_IVEC3(inxyz, nx, ny, 1); /* nxyz */
1051 origin.xyz[0] = origin.xyz[1] = -1.0; /* origin */
1052 origin.xyz[2] = 0.0;
1053 delta.xyz[0] = delta.xyz[1] = 2.0/(nx-1); /* delta */
1054 delta.xyz[2] = 1.0;
1055
1056 EDIT_dset_items(dout, ADN_nxyz, inxyz,
1057 ADN_xyzorg, origin,
1058 ADN_xyzdel, delta,
1059 ADN_prefix, fname,
1060 ADN_nvals, 2,
1061 ADN_none);
1062
1063 /* first is gaussian, second is first mask (to compare orientations) */
1064 /* use NULL to create space, then copy results */
1065 EDIT_substitute_brick(dout, 0, MRI_float, NULL);
1066 EDIT_substitute_brick(dout, 1, MRI_float, NULL);
1067
1068 /* first copy 'curve' to slice 0 */
1069 mptr = DBLK_ARRAY(dout->dblk, 0);
1070 dptr = curve;
1071 for(ind = 0; ind < nx*ny; ind++)
1072 *mptr++ = *dptr++;
1073
1074 /* now fill mask slice */
1075 mptr = DBLK_ARRAY(dout->dblk, 1);
1076 bptr = DBLK_ARRAY(g_saset->dblk, 0);
1077 for(ind = 0; ind < nx*ny; ind++, bptr++)
1078 *mptr++ = (float)*bptr;
1079
1080 dout->daxes->xxorient = ORI_L2R_TYPE;
1081 dout->daxes->yyorient = ORI_P2A_TYPE;
1082 dout->daxes->zzorient = ORI_I2S_TYPE;
1083
1084 if( hist ) tross_Append_History(dout, hist);
1085
1086 fprintf(stderr,", writing");
1087 DSET_write(dout);
1088
1089 DSET_delete(dout);
1090 fprintf(stderr,", done\n");
1091
1092 return 0;
1093 }
1094
1095
1096 /* ------------------------------------------------------------ */
1097 /* A = [R^2cos^2(theta) + sin^2(theta)] / [2R^2sigma^2]
1098 * B = -(R^2-1) * sin(2theta) / [4R^2sigma^2]
1099 * C = [R^2sin^2(theta) + cos^2(theta)] / [2R^2sigma^2]
1100 */
get_ABC(float sigma,float sigrat,float theta,double * A,double * B,double * C)1101 static int get_ABC(float sigma, float sigrat, float theta,
1102 double * A, double * B, double * C)
1103 {
1104 double R2, R2S2, C2, S2, So2;
1105
1106 R2 = sigrat*sigrat;
1107 R2S2 = 2*R2*sigma*sigma;
1108 C2 = cos(theta)*cos(theta);
1109 S2 = sin(theta)*sin(theta);
1110 So2 = sin(2*theta);
1111
1112 *A = (R2 * C2 + S2) / R2S2;
1113 *B = -(R2 - 1.0) * So2 / (2.0 * R2S2);
1114 *C = (R2 * S2 + C2) / R2S2;
1115
1116 return 0;
1117 }
1118
1119
1120 /* - compute an elliptic Gaussian curve over the current grid,
1121 * centered at x0, y0, with the given sigma, sigrat and theta
1122 *
1123 * old: fill with e^-[((x-x0)^2 + (y-y0)^2) / (2*sigma^2)]
1124 * new: e^-[A(x-x0)^2 + B(x-x0)(y-y0) + C(y-y0)^2], where
1125 * A = [R^2cos^2(theta) + sin^2(theta)] / [2R^2sigma^2]
1126 * B = -(R^2-1) * sin(2theta) / [4R^2sigma^2]
1127 * C = [R^2sin^2(theta) + cos^2(theta)] / [2R^2sigma^2]
1128 *
1129 * We do not have to be too efficient in computing A,B,C, since those
1130 * are constant across the image. Only x-x0 and y-y0 vary.
1131 *
1132 * Compute:
1133 * R^2, 2R^2sigma^2,
1134 * cos^2(theta), sin^2(theta), sin(2theta)
1135 */
compute_e_x_grid(float * e,int nx,int ny,float x0,float y0,float sigma,float sigrat,float theta)1136 static int compute_e_x_grid(float * e, int nx, int ny, float x0, float y0,
1137 float sigma, float sigrat, float theta)
1138 {
1139 float * eptr, eval;
1140 double wscale, sum;
1141 double xoff, yoff, expval;
1142 int ix, iy, eind;
1143
1144 double A, B, C;
1145
1146 wscale = 2.0/(nx-1.0); /* scale [0,nx-1] to [0,2] */
1147
1148 get_ABC(sigma, sigrat, theta, &A, &B, &C); /* get exp coefficients */
1149
1150 eptr = e; /* base pointer */
1151 sum = 0.0; /* for scaling to unit integral */
1152 for( iy = 0; iy < ny; iy++ ) {
1153 for( ix = 0; ix < nx; ix++ ) {
1154 xoff = ix*wscale - 1.0f - x0; /* map to [-1,1] and dist from x0 */
1155 yoff = iy*wscale - 1.0f - y0;
1156
1157 /* compute (positive) power of e, will scale by g_exp_ipieces */
1158 eval = A*xoff*xoff + B*xoff*yoff + C*yoff*yoff;
1159 if( eval > g_exp_maxval ) {
1160 *eptr++ = 0.0f;
1161 continue;
1162 }
1163
1164 if( eval < 0.0f ) eval = 0.0f;
1165
1166 /* could do expval = exp(-eval); */
1167 eind = (int)(eval*g_exp_ipieces); /* truncate towards 0? */
1168 if ( eind < g_exp_nvals ) {
1169 expval = g_exp_ts[eind];
1170 sum += expval;
1171 } else expval = 0.0f;
1172 *eptr++ = expval;
1173 }
1174 }
1175
1176 /* and rescale by 1/sum */
1177 if( sum != 0.0 ) {
1178 sum = 1.0/sum;
1179 for( ix = 0, eptr = e; ix < nx*ny; ix++, eptr++ )
1180 if( *eptr ) *eptr *= sum;
1181 }
1182
1183 return 0;
1184 }
1185
1186 /*----------------------------------------------------------------------*/
model_help(void)1187 static int model_help(void)
1188 {
1189 printf(
1190 "----------------------------------------------------------------------\n"
1191 "PRF_6_BAD - 6 parameter population receptive field (in visual cortex)\n"
1192 "\n"
1193 " ------------------------------------------------------------\n"
1194 " This is the bad, naughty version of PRF_6, which is missing a\n"
1195 " scalar of 2 on the B term in the main PRF model equation.\n"
1196 " It is left for comparison. Sorry Ed.\n"
1197 "\n"
1198 " This also applies AFNI_MODEL_PRF_GAUSS_FILE, if requested\n"
1199 " ------------------------------------------------------------\n"
1200 "\n"
1201 " Given stimulus images over time s(x,y,t), find x0, y0, sigma, R and\n"
1202 " theta values that produce a best fit of the model to the data. Here\n"
1203 " x0, y0 are taken to be the center of the population receptive field,\n"
1204 " sigma is the minor width of it (sigma_x, below), sigrat R is the ratio\n"
1205 " (sigma_y / sigma_x), and theta is the rotation from the y-direction\n"
1206 " major axis (so zero is in the positive y-direction).\n"
1207 "\n"
1208 " We assume sigma_y >= sigma_x and refer to sigrat >= 1, since that\n"
1209 " sufficiently represents all possibilities. The reciprocol would\n"
1210 " come from the negative complimentary angle, and would therefore be a\n"
1211 " redundant solution.\n"
1212 "\n"
1213 " parameter domains:\n"
1214 " x,y : [-1,1], scaled by the mask, itself\n"
1215 " sigma : (0,1], where 1 means the mask radius\n"
1216 " R (sigrat) : [1,inf), since sigma defines the smaller size\n"
1217 " theta : [-PI/2, PI/2), since rotation by PI has no effect\n"
1218 "\n"
1219 " The model function of x0, y0, sigma, R and theta is constructed as\n"
1220 " follows:\n"
1221 "\n"
1222 " 1. generate a 2-D elliptical Gaussian density function,\n"
1223 " centered at x0, y0, with given sigma, R (=sigma_y/sigma_x),\n"
1224 " and theta (rotation of major direction from positive y):\n"
1225 "\n"
1226 " -> pRF model g(x,y) = generalized 2-D Gaussian\n"
1227 "\n"
1228 " e^-(A(x-x0)^2 + B(x-x0)(y-y0) + C(y-y0)^2), where\n"
1229 "\n"
1230 " cos^2(theta) sin^2(theta)\n"
1231 " A = ------------ + ------------\n"
1232 " 2sigma_x^2 2sigma_y^2\n"
1233 "\n"
1234 " sin(2theta) sin(2theta)\n"
1235 " B = - ----------- + -----------\n"
1236 " 4sigma_x^2 4sigma_y^2\n"
1237 "\n"
1238 " sin^2(theta) cox^2(theta)\n"
1239 " C = ------------ + ------------\n"
1240 " 2sigma_x^2 2sigma_y^2\n"
1241 "\n"
1242 " Substituting sigma_x = sigma, sigma_y = Rsigma_x yields,\n"
1243 " \n"
1244 " R^2cos^2(theta) + sin^2(theta)\n"
1245 " A = ------------------------------\n"
1246 " 2R^2sigma^2\n"
1247 "\n"
1248 " sin(2theta)\n"
1249 " B = -(R^2-1) -----------\n"
1250 " 4R^2sigma^2\n"
1251 "\n"
1252 " R^2sin^2(theta) + cos^2(theta)\n"
1253 " C = ------------------------------\n"
1254 " 2R^2sigma^2\n"
1255 "\n"
1256 "--------------------------------------------------\n"
1257 "To use this model function:\n"
1258 "\n"
1259 " 1. Generate the stimulus time series (currently, images must be square).\n"
1260 "\n"
1261 " This should be a 2D+time dataset of visual stimuli over time. They\n"
1262 " are viewed as binary masks by the model function.\n"
1263 "\n"
1264 " * If results are computed on a restricted grid (which is much faster\n"
1265 " and is the default (see AFNI_MODEL_PRF_ON_GRID)), the resolution of\n"
1266 " those X,Y results will come directly from this stimulus dataset.\n"
1267 " It might be reasonable to have this be 100 or 200 (or 101 or 201)\n"
1268 " voxels on a side.\n"
1269 "\n"
1270 "\n"
1271 " **** Years have gone by, and I still have not implemented the speed-up\n"
1272 " for AFNI_MODEL_PRF_ON_GRID. Please leave it set to NO for now.\n"
1273 "\n"
1274 "\n"
1275 " * The amount of memory used for the precomputation should be the size\n"
1276 " of this dataset (in float format) times AFNI_MODEL_PRF_SIGMA_NSTEPS.\n"
1277 " It is converted to floats because it is blurred internally.\n"
1278 " The default AFNI_MODEL_PRF_SIGMA_NSTEPS is 100.\n"
1279 "\n"
1280 " 2. Scale and demean the input EPI time series data.\n"
1281 "\n"
1282 " Scaling is done to put the amplitude values into a reasonable (i.e.\n"
1283 " expected) range, such as by scaling it to a fraction of the mean\n"
1284 " (or maybe twice that).\n"
1285 "\n"
1286 " Setting the mean to zero is done so that no baseline modeling is\n"
1287 " needed (though it might be good to model drifts in the future).\n"
1288 "\n"
1289 " 3. Generate a convolution reference time series, such as one for GAM.\n"
1290 " This should be on the same TR grid, which is 2 in this example.\n"
1291 "\n"
1292 " 3dDeconvolve -nodata 10 2 -polort -1 \\\n"
1293 " -num_stimts 1 -stim_times 1 '1D:0' GAM \\\n"
1294 " -x1D conv.ref.GAM.1D\n"
1295 "\n"
1296 " 4. Set up environment variables to control execution:\n"
1297 "\n"
1298 " setenv AFNI_CONVMODEL_REF conv.ref.GAM.1D\n"
1299 " setenv AFNI_MODEL_PRF_STIM_DSET stim.144.bmask.resam+tlrc\n"
1300 "\n"
1301 " 5. And execute:\n"
1302 "\n"
1303 " 3dNLfim -input epi.scale.demean+tlrc \\\n"
1304 " -noise Zero \\\n"
1305 " -signal Conv_PRF_6_BAD \\\n"
1306 " -sconstr 0 -10.0 10.0 \\\n"
1307 " -sconstr 1 -1.0 1.0 \\\n"
1308 " -sconstr 2 -1.0 1.0 \\\n"
1309 " -sconstr 3 0.0 1.0 \\\n"
1310 " -sconstr 4 1.0 4.0 \\\n"
1311 " -sconstr 5 -1.571 1.570 \\\n"
1312 " -BOTH \\\n"
1313 " -nrand 10000 \\\n"
1314 " -nbest 5 \\\n"
1315 " -bucket 0 buck.PRF \\\n"
1316 " -snfit snfit.PRF\n"
1317 "\n"
1318 "--------------------------------------------------\n"
1319 "environment variables:\n"
1320 "\n"
1321 " -----------------------------------\n"
1322 " required:\n"
1323 "\n"
1324 " AFNI_CONVMODEL_REF : specify convolution reference file\n"
1325 "\n"
1326 " e.g. setenv AFNI_CONVMODEL_REF conv.ref.GAM.1D\n"
1327 "\n"
1328 " The file this specifies should contain a (short?) impulse\n"
1329 " response function, such as made in step #3, above.\n"
1330 "\n"
1331 " AFNI_MODEL_PRF_STIM_DSET : specify visual stimulus dataset\n"
1332 "\n"
1333 " e.g. setenv AFNI_MODEL_PRF_STIM_DSETstim.144.bmask.resam+tlrc\n"
1334 "\n"
1335 " This should be a 2D+time dataset of stimulus images over TRs.\n"
1336 " It will be converted to a byte mask over the visual field.\n"
1337 "\n"
1338 " -----------------------------------\n"
1339 " optional (for use with pre-computed grid):\n"
1340 "\n"
1341 " AFNI_MODEL_PRF_ON_GRID : Y/N - use pre-computed solutions\n"
1342 "\n"
1343 "\n"
1344 " *** This has not been coded yet. If there is strong interest,\n"
1345 " please let me know, though no promises are being made.\n"
1346 "\n"
1347 "\n"
1348 " AFNI_MODEL_PRF_SIGMA_MAX : specify maximum allowable sigma\n"
1349 "\n"
1350 " e.g. setenv AFNI_MODEL_PRF_SIGMA_MAX 2.0\n"
1351 " e.g. default 1.0\n"
1352 "\n"
1353 " Applies directly to AFNI_MODEL_PRF_ON_GRID.\n"
1354 "\n"
1355 " Use this variable to set the maximum pre-computed sigma.\n"
1356 " This should probably match the sconstr value for sigma.\n"
1357 "\n"
1358 " AFNI_MODEL_PRF_SIGMA_NSTEPS : specify number of pre-computed sigmas\n"
1359 "\n"
1360 " e.g. setenv AFNI_MODEL_PRF_SIGMA_NSTEPS 50\n"
1361 " e.g. default 100\n"
1362 "\n"
1363 " Applies directly to AFNI_MODEL_PRF_ON_GRID.\n"
1364 "\n"
1365 " Use this variable to set the number of pre-computed sigma values.\n"
1366 " Note that the resolution of pre-computed sigma values will be the\n"
1367 " ratio: AFNI_MODEL_PRF_SIGMA_MAX/AFNI_MODEL_PRF_SIGMA_NSTEPS.\n"
1368 "\n"
1369 " -----------------------------------\n"
1370 " helpful:\n"
1371 "\n"
1372 " AFNI_MODEL_HELP_CONV_PRF_6_BAD : Y/N - output this help\n"
1373 "\n"
1374 " e.g. setenv AFNI_MODEL_HELP_CONV_PRF_6_BAD YES\n"
1375 "\n"
1376 " When set, the model initialization function will output this help.\n"
1377 "\n"
1378 " Consider:\n"
1379 "\n"
1380 " 3dNLfim -signal Conv_PRF_6_BAD\n"
1381 "\n"
1382 " or more directly (without setenv):\n"
1383 "\n"
1384 " 3dNLfim -DAFNI_MODEL_HELP_CONV_PRF_6_BAD=Y -signal Conv_PRF_6_BAD \n"
1385 "\n"
1386 " AFNI_MODEL_DEBUG : specify debug/verbosity level\n"
1387 "\n"
1388 " e.g. setenv AFNI_MODEL_DEBUG 2\n"
1389 "\n"
1390 " Be more verbose. Valid levels are from 0 to 3, currently.\n"
1391 "\n"
1392 " AFNI_MODEL_DITER : specify debug iteration\n"
1393 "\n"
1394 " e.g. setenv AFNI_MODEL_DITER 999\n"
1395 "\n"
1396 " Get extra debug info at some iteration.\n"
1397 "\n"
1398 " AFNI_MODEL_PRF_GAUSS_FILE : specify dataset prefix for Gauss curve\n"
1399 "\n"
1400 " e.g. setenv AFNI_MODEL_PRF_GAUSS_FILE guass_curve\n"
1401 "\n"
1402 " Write a 2-D image with the Gaussian curve, which is helpful\n"
1403 " for checking the parameters. This works best when used via\n"
1404 " the get_afni_model_PRF_6 program.\n"
1405 "\n"
1406 "----------------------------------------------------------------------\n"
1407 " Written for E Silson and C Baker.\n"
1408 "\n"
1409 " R. Reynolds 27 June, 2014\n"
1410 "----------------------------------------------------------------------\n"
1411 );
1412
1413 return 0 ;
1414 }
1415
1416