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