1 /*****************************************************************************
2
3 Common functions for the 4-parameter PRF models.
4
5 These functions can be applied in the basic 4-param version or for the
6 difference of Gaussians version, which has 6 parameters, but is just
7 computed as the difference between 2 4-parameter expressions.
8
9 So at this level, there is little difference between the 2 cases.
10
11 ----------------------------------------------------------------------
12
13 The model is from:
14
15 Population receptive field estimates in human visual cortex
16 NeuroImage 39 (2008) 647-660
17 Serge O. Dumoulin, Brian A. Wandell
18
19 Given stimulus images over time s(x,y,t), find x0, y0, sigma values that
20 produce a best fit of the model to the data, where the model function of
21 x0, y0 and sigma is constructed as follows:
22
23 1. generate a 2-D Gaussian density function, centered at x0, y0,
24 and with given sigma
25 -> pRF model g(x,y) = e^-([(x-x0)^2+(y-y0)^2] / (2*sigma^2))
26 2. integrate (dot product) over binary stimulus image per time point
27 -> pRF response r(t)
28 3. convolve with HRF
29 -> Prediction p(t)
30
31 Then perform a simple regression fit, returning A (amplitude), the input
32 x0, y0, sigma values, and the error (sum squares of the residuals)
33
34 inputs
35
36 s(x,y,t) : 2d+time visual field bitmap timeseries
37 (e.g. rotating stimulus wedge)
38 hrf : HRF convolution curve - assume GAM for now
39 data : time series, per voxel
40
41 outputs
42
43 A : amplitude of fit
44 x, y, sigma : best fit coordinates of center of Gaussian kernel
45
46 R. Reynolds June, 2014
47 ******************************************************************************/
48
49 static int refnum = 0; /* # pts in refts */
50 static int refnz = 0; /* # of nonzero pts */
51 static float * refts = NULL; /* reference time series */
52 static int * refin = NULL; /* indexes of nonzero pts */
53 static int g_iter = -1; /* iteration number */
54
55
56 /* exp variables, parameters */
57 static float g_exp_maxval = 8.0; /* max x in exp(-x) */
58 static int g_exp_ipieces = 1000; /* vals per unit length */
59
60 /* exp variables, implied */
61 static int g_exp_nvals = 0; /* maxval*ipieces + 1 */
62 static float * g_exp_ts = NULL; /* exp(-x) for x at index VAL*pieces */
63
64
65 static THD_3dim_dataset * g_saset=NULL; /* stimulus aperture dataset */
66
67 /* prototypes */
68 static int reset_stim_aperture_dset(int);
69 static int reset_exp_time_series(void);
70 static int show_malloc_stats(char * mesg);
71
72 static int convert_to_blurred_masks(THD_3dim_dataset *, THD_3dim_dataset *);
73 static THD_3dim_dataset * THD_reorg_dset(THD_3dim_dataset*, THD_3dim_dataset*);
74 static THD_3dim_dataset * alloc_reorg_dset(THD_3dim_dataset * din, int dsetnz);
75 static int convolve_dset(THD_3dim_dataset * tset);
76 static float * get_float_volume_copy(THD_3dim_dataset *dset, int index, int nz);
77 static int compute_e_x_grid(float *e, int nx, int ny, float x0, float y0,
78 float sigma);
79 static int get_signal_computed(float * ts, int tslen, THD_3dim_dataset * dset,
80 float x0, float y0, float sigma, float A, int debug);
81 static int get_signal_from_grid(float * fdest, int nt, THD_3dim_dataset * dsrc,
82 float x, float y, float sigma, float scale, int debug);
83 static int inputs_to_coords(THD_3dim_dataset * dset, float x, float y,
84 float sigma);
85 static int insert_reorg_slice(THD_3dim_dataset * rset, float * slice,
86 int din_nxy, int din_zind, int din_tind);
87
88 static int disp_floats(char * mesg, float * p, int len);
89 static int write_gauss_file(char * fname, float * curve,
90 int nx, int ny, char * hist);
91 static int write_dset(THD_3dim_dataset * dset, char * name);
92 static int convolve_by_ref(float *, int, float *, int, int, int);
93
94 static void conv_model( float * gs , int ts_length ,
95 float ** x_array , float * ts_array );
96
97 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
98
99 /* MATH uses exp(x) directly, while ARRAY uses pre-computed values */
100 #define ACCUM_EXP_MATH(val) do { if(val>=0) result += exp(-val); } while (0)
101 #define ACCUM_EXP_ARRAY(val) do { \
102 eind = (int)(val*g_exp_ipieces); /* truncate towards 0? */ \
103 if ( eind >= 0 && eind <= g_exp_nvals ) \
104 result += g_exp_ts[eind]; } while (0)
105
106 /* ---------------------------------------------------------------------- */
107 /* interface to the environment */
108 static char * genv_conv_ref = NULL; /* AFNI_CONVMODEL_REF */
109 static char * genv_prf_stim = NULL; /* AFNI_MODEL_PRF_STIM_DSET */
110 static char * genv_gauss_file = NULL; /* AFNI_MODEL_PRF_GAUSS_FILE */
111
112 static int genv_diter = -1; /* debug iteration */
113 static int genv_debug = 1; /* AFNI_MODEL_DEBUG */
114 static int genv_ram_stats= 0; /* AFNI_MODEL_PRF_RAM_STATS */
115
116 static int genv_on_grid = 0; /* restrict computations to grid */
117 static float genv_sigma_max = 1.0; /* on_grid: maximum blur sigma */
118 static int genv_sigma_nsteps = 100; /* on_grid: number of blur steps */
119
120 static int genv_get_help = 0; /* AFNI_MODEL_HELP_ALL or HELP_CONV_PRF */
121
122 /* possble values for genv_ram_stats, to show malloc_stats() or ps */
123 #define AC_PRF_STAT_NONE 0
124 #define AC_PRF_STAT_MALLOC 1
125 #define AC_PRF_STAT_PS 2
126 #define AC_PRF_STAT_ALL (AC_PRF_STAT_MALLOC|AC_PRF_STAT_PS)
127 #define AC_PRF_STAT_WAIT 4
128
129 /* return one of the AC_PRF_STAT values */
set_env_ram_stats(void)130 static int set_env_ram_stats(void){
131 char * estr = my_getenv("AFNI_MODEL_PRF_RAM_STATS");
132 int rval = AC_PRF_STAT_NONE;
133
134 if( ! estr ) rval = AC_PRF_STAT_NONE;
135 else if( AFNI_noenv("AFNI_MODEL_PRF_RAM_STATS") ) rval = AC_PRF_STAT_NONE;
136 else if( AFNI_yesenv("AFNI_MODEL_PRF_RAM_STATS") ) rval = AC_PRF_STAT_ALL;
137 else if( ! strcmp(estr, "MALLOC") ) rval = AC_PRF_STAT_MALLOC;
138 else if( ! strcmp(estr, "PS") ) rval = AC_PRF_STAT_PS;
139 else if( ! strcmp(estr, "ALL") ) rval = AC_PRF_STAT_ALL;
140 else if( ! strcmp(estr, "WAIT") )
141 rval = AC_PRF_STAT_WAIT | AC_PRF_STAT_ALL;
142 else rval = (int)AFNI_numenv_def("AFNI_MODEL_PRF_RAM_STATS", 0);
143
144 if( ! rval ) { genv_ram_stats = rval; return 0; }
145
146 if( genv_debug ) fprintf(stderr,"-- setting PRF_RAM_STATS to %d\n", rval);
147
148 genv_ram_stats = rval;
149 return 0;
150 }
151
set_env_vars(void)152 static int set_env_vars(void)
153 {
154 genv_conv_ref = my_getenv("AFNI_CONVMODEL_REF"); /* reference file */
155 if( genv_conv_ref ) fprintf(stderr,"-- PRF: have REF %s\n", genv_conv_ref);
156 else fprintf(stderr,"** model PRF: AFNI_CONVMODEL_REF is not set\n");
157
158 genv_prf_stim = my_getenv("AFNI_MODEL_PRF_STIM_DSET"); /* visual stim */
159 if( genv_prf_stim ) fprintf(stderr,"-- PRF: have stim %s\n", genv_prf_stim);
160 else fprintf(stderr,"** model PRF: AFNI_MODEL_PRF_STIM_DSET is not set\n");
161
162 genv_diter = (int)AFNI_numenv_def("AFNI_MODEL_DITER", -1);
163 genv_debug = (int)AFNI_numenv_def("AFNI_MODEL_DEBUG", 1);
164 if( genv_debug )
165 fprintf(stderr,"-- PRF: debug %d, iter %d\n", genv_debug, genv_diter);
166
167 /* on grid - default to yes */
168 genv_on_grid = 1-AFNI_noenv("AFNI_MODEL_PRF_ON_GRID"); /* flag */
169 if( genv_debug )
170 fprintf(stderr,"-- PRF: results on grid: %s\n", genv_on_grid?"yes":"no");
171
172 genv_sigma_max = AFNI_numenv_def("AFNI_MODEL_PRF_SIGMA_MAX", genv_sigma_max);
173 genv_sigma_nsteps = (int)AFNI_numenv_def("AFNI_MODEL_PRF_SIGMA_NSTEPS",
174 genv_sigma_nsteps);
175 if( genv_debug && genv_on_grid )
176 fprintf(stderr,"-- PRF: sigma_max = %f, nsteps = %d\n",
177 genv_sigma_max, genv_sigma_nsteps);
178
179 /* show RAM stats - default to no */
180 set_env_ram_stats();
181
182 /* help */
183 genv_get_help = AFNI_yesenv("AFNI_MODEL_HELP_CONV_PRF")
184 || AFNI_yesenv("AFNI_MODEL_HELP_ALL");
185
186 /* write a Gaussian mask? */
187 genv_gauss_file = my_getenv("AFNI_MODEL_PRF_GAUSS_FILE");
188 if( genv_gauss_file && genv_debug )
189 fprintf(stderr, "-- plan to write gauss file %s\n", genv_gauss_file);
190
191 return 0;
192 }
193 /* ---------------------------------------------------------------------- */
194
195
196 /*----------------------------------------------------------------------
197 Function to set the reference time series, with which the
198 model function is convolved to produce the simulated data.
199 ------------------------------------------------------------------------*/
200
conv_set_ref(int num,float * ref)201 static void conv_set_ref( int num , float * ref )
202 {
203 if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
204 int ii ;
205
206 /* get rid of old data */
207
208 if(refts != NULL){ free(refts); refts = NULL; free(refin); refin = NULL; }
209
210 refnum = num ;
211 refts = (float *) malloc( sizeof(float) * num ) ;
212 refin = (int *) malloc( sizeof(int) * num ) ;
213 memcpy( refts , ref , sizeof(float) * num ) ;
214 for( ii=0,refnz=0 ; ii < num ; ii++ ) /* build list of nonzero */
215 if( refts[ii] != 0.0 ) refin[refnz++] = ii ; /* points in refts */
216 if( refnz == 0 )
217 ERREX("model_conv_PRF: All zero reference timeseries!") ;
218
219 if( genv_debug ) {
220 fprintf(stderr,"+d conv_set_ref: num=%d nonzero=%d\n",num,refnz) ;
221 if( genv_debug > 1 ) {
222 fprintf(stderr," TR locked stimuli :");
223 for( ii = 0; ii < refnz; ii++ )
224 fprintf(stderr," %f", refts[refin[ii]]);
225 fputc('\n',stderr);
226 }
227 }
228
229 return;
230
231 } else { /*** if no inputs, read it from AFNI_CONVMODEL_REF 1D file ***/
232
233 MRI_IMAGE * flim ;
234
235 if( genv_conv_ref == NULL )
236 ERREX("model_conv_PRF: need ref file as AFNI_CONVMODEL_REF") ;
237
238 flim = mri_read_1D(genv_conv_ref);
239 if( flim == NULL ){
240 fprintf(stderr,"** model_conv_PRF: Can't read timeseries file %s",
241 genv_conv_ref);
242 ERREX("failing...");
243 }
244
245 if( genv_debug ) fprintf(stderr,"+d conv_set_ref: refts=%s nx=%d\n",
246 genv_conv_ref, flim->ny);
247
248 conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ; /* recursion! */
249 mri_free(flim) ;
250 }
251 return;
252 }
253
254 /* allocate reorg dset early, to allow free() of all temp dsets
255 * back to OS
256 *
257 * dsetnz : if set, get nz from dset, else get nz from genv_sigma_nsteps
258 */
alloc_reorg_dset(THD_3dim_dataset * din,int dsetnz)259 static THD_3dim_dataset * alloc_reorg_dset(THD_3dim_dataset * din, int dsetnz)
260 {
261 THD_3dim_dataset * dout;
262 float * newvol;
263 int tind;
264
265 THD_ivec3 iv_nxyz;
266 int in_nx, in_ny, in_nz, in_nt;
267 int out_nx, out_ny, out_nz, out_nt;
268
269 dout = EDIT_empty_copy(din);
270
271 in_nx = DSET_NX(din); in_ny = DSET_NY(din); in_nz = DSET_NZ(din);
272 if( !dsetnz ) in_nz = genv_sigma_nsteps;
273 in_nt = DSET_NVALS(din);
274
275 out_nx = in_nt; out_ny = in_nx; out_nz = in_ny; out_nt = in_nz;
276
277 if(genv_debug > 1)
278 fprintf(stderr,"-- alloc_reorg: nxyzt (%d,%d,%d,%d) -> (%d,%d,%d,%d)\n",
279 in_nx, in_ny, in_nz, in_nt, out_nx, out_ny, out_nz, out_nt);
280
281 LOAD_IVEC3(iv_nxyz , out_nx, out_ny, out_nz);
282
283 EDIT_dset_items(dout,
284 ADN_prefix, "alloc.reorg",
285 ADN_nxyz, iv_nxyz, ADN_nvals, out_nt, ADN_ntt, out_nt,
286 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
287 ADN_none);
288
289 /* create and attach each volume */
290 for( tind = 0; tind < DSET_NVALS(dout); tind++ ) {
291 newvol = (float *)malloc(DSET_NVOX(dout) * sizeof(float));
292 mri_fix_data_pointer(newvol, DSET_BRICK(dout, tind));
293 }
294
295 return dout;
296 }
297
298
write_dset(THD_3dim_dataset * dset,char * name)299 static int write_dset(THD_3dim_dataset * dset, char * name)
300 {
301 EDIT_dset_items(dset, ADN_prefix, name, ADN_none);
302 fprintf(stderr, "++ writing dset %s ...\n", name);
303 DSET_write(dset);
304
305 return 0;
306 }
307
308 /* for monitoring the RAM usage */
show_malloc_stats(char * mesg)309 static int show_malloc_stats(char * mesg)
310 {
311 int show_stats, show_ps, get_char;
312
313 show_stats = genv_ram_stats & AC_PRF_STAT_MALLOC;
314 show_ps = genv_ram_stats & AC_PRF_STAT_PS;
315 get_char = genv_ram_stats & AC_PRF_STAT_WAIT;
316
317 if( ! show_stats && ! show_ps ) return 0;
318
319 if( show_stats ) {
320 fprintf(stderr,"\n----- malloc stats: %s\n", mesg);
321 /* apply FreeBSD patch from J Bacon [10 Jun 2019 rickr] */
322 #if defined(__linux__)
323 #include <malloc.h>
324 malloc_stats();
325 #elif defined(__FreeBSD__) && !defined(__DragonFly__)
326 #include <stdlib.h>
327 #include <malloc_np.h>
328 malloc_stats_print(NULL, NULL, "g");
329 #else
330 fprintf(stderr, "No malloc_stats() on this platform.\n");
331 #endif
332 }
333
334 if( show_ps ) {
335 fprintf(stderr,"\n----- ps info: %s\n", mesg);
336 system("ps aux | grep NLfim | grep -v 'bin.time' | grep PRF");
337 }
338
339 if( get_char ) {
340 fprintf(stderr," ... hit <enter> to proceed");
341 getchar();
342 }
343
344 return 0;
345 }
346
347 /* any failure should leave g_saset == NULL
348 *
349 * return 0 on success */
reset_stim_aperture_dset(int needed_length)350 static int reset_stim_aperture_dset(int needed_length)
351 {
352 THD_3dim_dataset * sareorg;
353 int errs=0;
354
355 /* free and reload saset */
356 if( g_saset ) DSET_delete(g_saset);
357 g_saset = THD_open_dataset(genv_prf_stim);
358 if( ! g_saset ) return 1;
359
360 /* check for square dataset and sufficient nt */
361 if( fabs(DSET_DX(g_saset) - DSET_DY(g_saset)) > 0.0001 ) {
362 fprintf(stderr,"** PRF: stimset voxels not square (%f, %f)\n",
363 DSET_DX(g_saset), DSET_DY(g_saset));
364 errs++;
365 }
366
367 if( DSET_NX(g_saset) != DSET_NY(g_saset) ) {
368 fprintf(stderr,"** PRF: stimset not square (%d, %d)\n",
369 DSET_NX(g_saset), DSET_NY(g_saset));
370 errs++;
371 }
372
373 if( DSET_NVALS(g_saset) < needed_length ) {
374 fprintf(stderr,"** PRF: dset nt = %d, tslen = %d\n",
375 DSET_NVALS(g_saset), needed_length);
376 errs++;
377 }
378
379 if( errs ) { DSET_delete(g_saset); g_saset = NULL; return 1; }
380
381
382 if( THD_dset_to_mask(g_saset, 1.0, 0.0) ) return 1;
383
384 /* allocate reorg dset early */
385 sareorg = alloc_reorg_dset(g_saset, 0);
386
387 show_malloc_stats("pre blur");
388
389 if( genv_on_grid ) {
390 if( convert_to_blurred_masks(g_saset, sareorg) ) {
391 fprintf(stderr,"** failed blur/reorg, nuking dset\n");
392 DSET_delete(sareorg); sareorg = NULL;
393 }
394
395 show_malloc_stats("post blur");
396 DSET_delete(g_saset);
397 show_malloc_stats("post blur delete");
398
399 g_saset = sareorg;
400 // write_dset(g_saset, "sa.blur");
401
402 convolve_dset(g_saset);
403 // write_dset(g_saset, "sa.conv");
404
405 if( ! g_saset ) return 1;
406 }
407
408 return 0;
409 }
410
411
412 /* x-axis is convolution time axis, so for each row, convolve */
convolve_dset(THD_3dim_dataset * tset)413 static int convolve_dset(THD_3dim_dataset * tset)
414 {
415 float * result, * volbase, * dbase;
416 int nx, ny, nz, nt, nxy;
417 int xind, yind, zind, tind;
418
419 if( ! tset ) { fprintf(stderr, "** no tset to convolve\n"); return 1; }
420
421 if(genv_debug>1)fprintf(stderr, "++ starting convolution at time %6.1f\n",
422 0.001*NI_clock_time());
423
424 /* allocate room for computation and result */
425 nx = DSET_NX(tset); ny = DSET_NY(tset); nz = DSET_NZ(tset);
426 nt = DSET_NVALS(tset); nxy = nx * ny;
427 result = (float *)malloc(nx * sizeof(float));
428 if( !result ) {
429 fprintf(stderr, "** PRF: conv_dset: failed malloc of %d floats\n", nx);
430 return 1;
431 }
432
433 for( tind = 0; tind < nt; tind++ ) {
434 volbase = (float *)DSET_ARRAY(tset, tind);
435
436 /* for each row (the signal), convolve */
437 for( zind = 0; zind < nz; zind++ ) {
438 for( yind = 0; yind < ny; yind++ ) {
439 /* result = convolved signal (initialize and demean) */
440 dbase = volbase + yind*nx + zind*nxy;
441 convolve_by_ref(result, nx, dbase, nx, 1, 1);
442 /* and put result back in */
443 for( xind = 0; xind < nx; xind++ ) dbase[xind] = result[xind];
444 }
445 }
446 }
447
448 if(genv_debug>1)fprintf(stderr, "-- finished convolution at time %6.1f\n",
449 0.001*NI_clock_time());
450
451 return 0;
452 }
453
454
455 /* Convert 2D+time byte binary mask dataset into 3D+time float
456 * dataset of slices (blur_nsteps) at different blur levels
457 * (multiples of blur_res).
458 *
459 * input: 1 binary mask slice per TR
460 * output: duplicate dset, but where each slice is converted to an array of
461 * slices at different blur levels
462 *
463 * Resample (reorganize) output so that time is the fastest dimension,
464 * then x, y and sigma (sub-brick direction).
465 *
466 * given blur_max (e.g. 0.1)
467 * blur_nsteps (e.g. 50)
468 * compute blur_res (e.g. 0.002)
469 *
470 *------------------------------------------------------------
471 * allocate float slize (NX*NY), SLICE
472 * for each input TR ("volume" is just a slice)
473 * get float copy of slice (mri_to_float()), COPY
474 * for each blur level (genv_sigma_nsteps)
475 * memcpy(SLICE, COPY, NX*NY*sizeof(float));
476 * EDIT_blur_volume(NX, NY, NZ=1, dx, dy, dz, MRI_float, data, SIGMA);
477 * copy slice into reorg dset
478 * mri_free(COPY)
479 * free(SLICE)
480 *------------------------------------------------------------
481 */
convert_to_blurred_masks(THD_3dim_dataset * dset,THD_3dim_dataset * rset)482 static int convert_to_blurred_masks(THD_3dim_dataset * dset,
483 THD_3dim_dataset * rset)
484 {
485 MRI_IMAGE * fim;
486 float * fslice, sigma;
487 int nx, ny, nt;
488 int vind, sind;
489
490 /* test required inputs */
491 if( !dset ) return 1;
492 if( genv_sigma_max <= 0.0 || genv_sigma_max > 1.0 || genv_sigma_nsteps<=1 ){
493 fprintf(stderr,"** PRF on grid: bad sigma max %f or nsteps %d\n",
494 genv_sigma_max, genv_sigma_nsteps);
495 return 1;
496 }
497
498 if( DSET_NZ(dset) != 1 ) {
499 fprintf(stderr,"** invalid stim NZ = %d\n", DSET_NZ(dset));
500 return 1;
501 }
502
503 nx = DSET_NX(dset); ny = DSET_NY(dset); nt = DSET_NVALS(dset);
504
505 /* report RAM and time info */
506 if(genv_debug) { double mem;
507 fprintf(stderr, "++ making blurred time series: %d x %d x %d x nt=%d\n",
508 DSET_NX(dset),DSET_NY(dset), genv_sigma_nsteps, DSET_NVALS(dset));
509 fprintf(stderr, "++ stored as reorg dset: %d x %d x %d x nt=%d\n",
510 DSET_NX(rset),DSET_NY(rset), DSET_NZ(rset), DSET_NVALS(rset));
511 mem = (double)DSET_NVOX(rset)*DSET_NVALS(rset)*sizeof(float) / (1<<30);
512 fprintf(stderr,
513 " --> expected RAM for pre-computed results: %.2f GB\n\n", mem);
514 if(genv_debug)fprintf(stderr, "++ starting blur at time %6.1f\n",
515 0.001*NI_clock_time());
516 }
517
518 fslice = (float *)malloc(nx * ny * sizeof(float));
519 if( !fslice ){fputs("** CTBM N1\n", stderr); return 1;}
520 if( genv_debug ) fprintf(stderr, "-- blur/reorg %d images: 0", nt);
521 for( vind = 0; vind < nt; vind++ ) {
522 if( genv_debug ) fprintf(stderr, "\b\b\b%3d", vind);
523
524 /* get float copy of slice */
525 fim = mri_to_float(DSET_BRICK(dset, vind));
526 if(!fim){fputs("** CTBM N2\n", stderr); return 1;}
527
528 for( sind = 0; sind < genv_sigma_nsteps; sind++ ) {
529 /* start with original slice each time */
530 memcpy(fslice, MRI_FLOAT_PTR(fim), nx*ny*sizeof(float));
531
532 sigma = genv_sigma_max * ((sind + 1.0)/genv_sigma_nsteps);
533 FIR_blur_volume_3d(nx, ny, 1, 2.0/(nx-1), 2.0/(ny-1), 1.0,
534 fslice, sigma, sigma, 0.0);
535 insert_reorg_slice(rset, fslice, nx*ny, sind, vind);
536 }
537
538 mri_free(fim); /* nuke copy of slice */
539 }
540
541 free(fslice);
542
543 if(genv_debug)fprintf(stderr, "\n-- finished blur volume at time %6.1f\n",
544 0.001*NI_clock_time());
545
546 return 0;
547 }
548
insert_reorg_slice(THD_3dim_dataset * rset,float * slice,int din_nxy,int din_zind,int din_tind)549 static int insert_reorg_slice(THD_3dim_dataset * rset, float * slice,
550 int din_nxy, int din_zind, int din_tind)
551 {
552 float * outvol, * inptr;
553 int in_nx, in_ny, in_nz, in_nt;
554 int ix, iy;
555 int out_nx, out_ny, out_nxy, out_tind, out_xind;
556
557 /* Could pass in accum volume, and add tind+1 each time. The result
558 * should be the sum of genv_sigma_nsteps at each voxel (i.e. 5050 when
559 * genv_sigma_nsteps == 100). Or just count by 1 to 100, too lazy? */
560
561 /* get original dimensions, convert output indices */
562 in_nx = DSET_NY(rset); in_ny = DSET_NZ(rset); in_nz = DSET_NVALS(rset);
563 in_nt = DSET_NX(rset);
564 out_nx = DSET_NX(rset);
565 out_ny = DSET_NY(rset);
566 out_nxy = out_nx * out_ny;
567 out_tind = din_zind;
568 out_xind = din_tind;
569
570 if( in_nx * in_ny != din_nxy ) {
571 fprintf(stderr,"** IRS: no nyz match, nx = %d, ny = %d, DINnxy = %d\n",
572 in_nx, in_ny, din_nxy);
573 return 1;
574 }
575
576 if( din_zind >= in_nz || din_tind >= in_nt ) {
577 fprintf(stderr,"** IRS: no nzt match, nz=%d, zind=%d, nt=%d, tind=%d\n",
578 in_nz, din_zind, in_nt, din_tind);
579 return 1;
580 }
581
582 /* okay, copy the slice into new positions */
583 outvol = (float *)DSET_ARRAY(rset, out_tind);
584 inptr = slice;
585 for( iy = 0; iy < in_ny; iy++ )
586 for( ix = 0; ix < in_nx; ix++ )
587 /* index conversion: in_x -> out_y, in_y -> out_z */
588 outvol[out_xind + out_nx*ix + out_nxy * iy] = *inptr++;
589
590 return 0;
591 }
592
593
594 /* rotate list of axes and correspondingly shuffle data:
595 * x,y,z,t -> t,x,y,z
596 *
597 * The main point is to make the time axes the fast direction.
598 *
599 * require MRI_float for now */
THD_reorg_dset(THD_3dim_dataset * din,THD_3dim_dataset * dorg)600 static THD_3dim_dataset * THD_reorg_dset(THD_3dim_dataset * din,
601 THD_3dim_dataset * dorg)
602 {
603 THD_3dim_dataset * dout;
604 float * newvol, * inslice, * outbase, * inbase;
605 int in_nx, in_ny, in_nz, in_nt;
606 int out_nx, out_ny, out_nz, out_nt;
607 int in_nxyz, in_nxy, out_nxyz, out_nxy;
608 int xind, yind, zind, tind;
609
610 if( !din || !dorg ) { fprintf(stderr,"** reorg_dset:NULL\n"); return NULL; }
611
612 dout = dorg;
613
614 in_nx = DSET_NX(din); in_ny = DSET_NY(din); in_nz = DSET_NZ(din);
615 in_nt = DSET_NVALS(din);
616 in_nxy = in_nx * in_ny;
617 in_nxyz = in_nx * in_ny * in_nz;
618
619 out_nx = in_nt; out_ny = in_nx; out_nz = in_ny; out_nt = in_nz;
620 out_nxy = out_nx * out_ny;
621 out_nxyz = out_nx * out_ny * out_nz;
622
623 if(genv_debug > 1)
624 fprintf(stderr,"-- reorg_dset: nxyzt (%d,%d,%d,%d) -> (%d,%d,%d,%d)\n",
625 in_nx, in_ny, in_nz, in_nt, out_nx, out_ny, out_nz, out_nt);
626 if(genv_debug > 1)fprintf(stderr, "\n== reorg starting at %6.1f\n",
627 0.001*NI_clock_time());
628
629 if( out_nx != DSET_NX(dout) || out_ny != DSET_NY(dout) ||
630 out_nz != DSET_NZ(dout) || out_nt != DSET_NVALS(dout) ) {
631 fprintf(stderr,"** reorg mis-match, crash and burn\n");
632 }
633
634 EDIT_dset_items(dout, ADN_prefix, "reorg.boots", ADN_none);
635
636 /* create and attach each volume */
637 for( tind = 0; tind < out_nt; tind++ ) {
638 /* now newvol just points to existing memory */
639 newvol = (float *)DSET_ARRAY(dout, tind);
640
641 for( xind = 0; xind < out_nx; xind++ ) {
642
643 /* copy one slice of input across y and z directions */
644 inslice = ((float *)DSET_ARRAY(din, xind)) + tind*in_nxy;
645
646 for( yind = 0; yind < out_ny; yind++ ) {
647 outbase = newvol + xind + yind*out_nx;
648 inbase = inslice + yind;
649 for( zind = 0; zind < out_nz; zind++ )
650 outbase[zind*out_nxy] = inbase[zind*in_nx];
651 }
652
653 }
654 }
655
656 if(genv_debug > 1)fprintf(stderr, "\n== reorg finished at %6.1f\n",
657 0.001*NI_clock_time());
658
659 /* maybe dump a couple of time series */
660 if( genv_debug > 2 ) {
661 MRI_IMAGE * im;
662 float * fp;
663 int i, j, k;
664 i = in_nx/3; j = in_ny/3; k = in_nz/3;
665 im = THD_extract_series(i+j*in_nx+k*in_nxy, din, 0);
666 disp_floats("== ARY: sig [nxyz/3]: ", MRI_FLOAT_PTR(im), in_nt);
667 fp = ((float *)DSET_ARRAY(dout,k)) + i*out_nx + j*out_nxy;
668 disp_floats("== ARY: reorg : ", fp, out_nx);
669 mri_free(im);
670 }
671
672 return dout;
673 }
674
675
676 /* OLD version: saved for comparison
677 * (new version will input pre-allocated reorg dset)
678 *
679 * rotate list of axes and correspondingly shuffle data:
680 * x,y,z,t -> t,x,y,z
681 *
682 * The main point is to make the time axes the fast direction.
683 *
684 * require MRI_float for now */
THD_reorg_dset_old(THD_3dim_dataset * din)685 static THD_3dim_dataset * THD_reorg_dset_old(THD_3dim_dataset * din)
686 {
687 THD_3dim_dataset * dout;
688 THD_ivec3 iv_nxyz;
689 float * newvol, * inslice, * outbase, * inbase;
690 int in_nx, in_ny, in_nz, in_nt;
691 int out_nx, out_ny, out_nz, out_nt;
692 int in_nxyz, in_nxy, out_nxyz, out_nxy;
693 int xind, yind, zind, tind;
694
695 if( !din ) { fprintf(stderr,"** reorg_dset: missing input\n"); return NULL; }
696
697 dout = EDIT_empty_copy(din);
698
699 in_nx = DSET_NX(din); in_ny = DSET_NY(din); in_nz = DSET_NZ(din);
700 in_nt = DSET_NVALS(din);
701 in_nxy = in_nx * in_ny;
702 in_nxyz = in_nx * in_ny * in_nz;
703
704 out_nx = in_nt; out_ny = in_nx; out_nz = in_ny; out_nt = in_nz;
705 out_nxy = out_nx * out_ny;
706 out_nxyz = out_nx * out_ny * out_nz;
707
708 LOAD_IVEC3(iv_nxyz , out_nx, out_ny, out_nz);
709
710 EDIT_dset_items(dout, ADN_prefix, "reorg.boots",
711 ADN_nxyz, iv_nxyz, ADN_nvals, out_nt, ADN_ntt, out_nt,
712 ADN_none);
713
714 if(genv_debug)
715 fprintf(stderr,"-- reorg_dset: nxyzt (%d,%d,%d,%d) -> (%d,%d,%d,%d)\n",
716 in_nx, in_ny, in_nz, in_nt, out_nx, out_ny, out_nz, out_nt);
717 if(genv_debug>1)fprintf(stderr, "\n== reorg starting at %6.1f\n",
718 0.001*NI_clock_time());
719
720 /* create and attach each volume */
721 for( tind = 0; tind < out_nt; tind++ ) {
722 newvol = (float *)malloc(out_nxyz * sizeof(float));
723 if( !newvol ) {
724 fprintf(stderr,"** PRF reorg: failed to alloc volume %d\n", tind);
725 DSET_delete(dout);
726 return NULL;
727 }
728
729 for( xind = 0; xind < out_nx; xind++ ) {
730
731 /* copy one slice of input across y and z directions */
732 inslice = ((float *)DSET_ARRAY(din, xind)) + tind*in_nxy;
733
734 for( yind = 0; yind < out_ny; yind++ ) {
735 outbase = newvol + xind + yind*out_nx;
736 inbase = inslice + yind;
737 for( zind = 0; zind < out_nz; zind++ )
738 outbase[zind*out_nxy] = inbase[zind*in_nx];
739 }
740
741 }
742
743 mri_fix_data_pointer(newvol, DSET_BRICK(dout, tind));
744 }
745
746 if(genv_debug>1)fprintf(stderr, "\n== reorg finished at %6.1f\n",
747 0.001*NI_clock_time());
748
749 /* maybe dump a couple of time series */
750 if( genv_debug > 2 ) {
751 MRI_IMAGE * im;
752 float * fp;
753 int i, j, k;
754 i = in_nx/3; j = in_ny/3; k = in_nz/3;
755 im = THD_extract_series(i+j*in_nx+k*in_nxy, din, 0);
756 disp_floats("== ARY: sig [nxyz/3]: ", MRI_FLOAT_PTR(im), in_nt);
757 fp = ((float *)DSET_ARRAY(dout,k)) + i*out_nx + j*out_nxy;
758 disp_floats("== ARY: reorg : ", fp, out_nx);
759 mri_free(im);
760 }
761
762 return dout;
763 }
764
765
766 /* 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)767 static float * get_float_volume_copy(THD_3dim_dataset * dset, int index, int nz)
768 {
769 MRI_IMAGE * fim;
770 float * fdata;
771 int nxy, sind;
772
773 /* make space for new data */
774 nxy = DSET_NXY(dset);
775 fdata = (float *)malloc(nxy * nz * sizeof(float));
776 if( ! fdata ) {
777 fprintf(stderr,"** PRF: failed to alloc %d floats at index %d\n",
778 nxy * nz, index);
779 return NULL;
780 }
781
782 /* get slice to duplicate */
783 fim = mri_to_float( DSET_BRICK(dset, index) );
784 if( ! fim ) {
785 fprintf(stderr,"** PRF: failed mri_to_float\n");
786 free(fdata);
787 return NULL;
788 }
789
790 /* dupe over slices */
791 for( sind = 0; sind < nz; sind++ )
792 memcpy(fdata+sind*nxy, MRI_FLOAT_PTR(fim), nxy * sizeof(float));
793
794 mri_free(fim);
795
796 return fdata;
797 }
798
799
800 /* want e^-x computed out to x=7 (c^-7 goes below 0.001)
801 *
802 * g_exp_maxval = 7
803 * g_exp_ipieces = 1000
804 * so step = 1/1000
805 * g_exp_nvals = g_exp_maxval * g_exp_ipieces + 1 (for 0), plus a few
806 */
reset_exp_time_series(void)807 static int reset_exp_time_series(void)
808 {
809 int ind;
810 float resol = 1.0/g_exp_ipieces;
811
812 g_exp_nvals = (int)(g_exp_maxval * g_exp_ipieces) + 1;
813
814 if(genv_debug) fprintf(stderr, "-- exp nvals = %d, max = %f, pieces = %d\n",
815 g_exp_nvals, g_exp_maxval, g_exp_ipieces);
816
817 if( g_exp_ts ) free(g_exp_ts);
818 g_exp_ts = (float *)malloc(g_exp_nvals * sizeof(float));
819 if( ! g_exp_ts ) {
820 fprintf(stderr,"** failed to alloc %d floats, buring...\n", g_exp_nvals);
821
822 if( g_saset ) DSET_delete(g_saset);
823 g_saset = NULL; /* this blocks analysis */
824 return 1;
825 }
826
827 for( ind = 0; ind < g_exp_nvals; ind++ )
828 g_exp_ts[ind] = exp(-resol*ind);
829
830 return 0;
831 }
832
833 /*-----------------------------------------------------------------------
834 Function to compute the simulated time series.
835 -------------------------------------------------------------------------*/
836
conv_model(float * gs,int ts_length,float ** x_array,float * ts_array)837 static void conv_model( float * gs , int ts_length ,
838 float ** x_array , float * ts_array )
839 {
840 int ii, cur_debug = 0, irfdur=0;
841 int iter_step;
842
843 static int nid = 0 ; /* number of pts in impulse */
844 static float * fid = NULL; /* impulse response function */
845
846 /*----- check for env vars -----*/
847
848 g_iter++ ;
849
850 if( g_iter == 0 ) {
851 set_env_vars(); /* process environment variables */
852 if(genv_debug && x_array) fprintf(stderr,"\n+d TR = %f\n",
853 x_array[1][1]-x_array[0][1]);
854 fprintf(stderr,"++ %s\n", g_model_ver);
855 }
856
857 /*** make sure there is a reference function to convolve with ***/
858 /* it may be used in reset_stim_aperture_dset */
859
860 if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
861
862 /* create stim aperture dset */
863 if( g_iter == 0 ) {
864 (void)reset_stim_aperture_dset(ts_length); /* free and reload saset */
865 (void)reset_exp_time_series(); /* pre-compute exp(x) */
866 if( genv_debug ) fprintf(stderr, "== start time %d\n", NI_clock_time());
867 }
868
869 /*** initialize the output before possible early return ***/
870
871 for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
872
873 /* if we had some failure, bail */
874 if( !g_saset ) return;
875
876 if( genv_debug ) { /* babble */
877 if( genv_on_grid ) iter_step = 100000;
878 else iter_step = 100;
879 if( (g_iter % iter_step) == 0 ) {
880 if( genv_debug > 1 )
881 fprintf(stderr, "-- time for %d iter set %5d : %6.1f\n",
882 iter_step, g_iter/iter_step, 0.001*NI_clock_time());
883 else {
884 static int first = 1;
885 if( first ) {
886 fprintf(stderr, "-- iteration size %d .", iter_step);
887 first = 0;
888 } else fputc('.', stderr);
889 }
890 }
891 }
892
893 /* to clean up, particularly as a parameter */
894 cur_debug = (g_iter == genv_diter || (g_iter == 0 && genv_debug > 1));
895 if( cur_debug ) disp_floats("+d input params: ", gs, 4);
896
897 /*** initialize the impulse response ***/
898
899 if( nid < ts_length ){ /* make some space for it */
900 if( fid ) free(fid) ;
901 nid = ts_length ;
902 fid = (float *) malloc( sizeof(float) * nid ) ;
903 }
904
905 /* If on grid, get signal directly. Else, convolve signal impulse
906 function with ref time series, and demean. */
907 if( genv_on_grid ) {
908 irfdur = signal_model(gs, ts_length, x_array, ts_array, cur_debug);
909 } else {
910 irfdur = signal_model(gs, ts_length, x_array, fid, cur_debug);
911 convolve_by_ref(ts_array, ts_length, fid, irfdur, 0, 1);
912 if( cur_debug ) disp_floats("+d no_grid conv : ", ts_array, ts_length);
913 }
914
915 return ;
916 }
917
918 /*
919 Convolve reference function (GAM here?) with stim time series.
920
921 globals:
922 int refnum : # pts in refts
923 int refnz : # of nonzero pts
924 float * refts : reference time series
925 int * refin : indexes of nonzero pts
926
927 inputs:
928 result : convolved result (already initialized to 0)
929 rlen : length of input data time series (maybe equal)
930 signal : signal model output
931 siglen : length of signal time series
932 init : flag: init result?
933 demean : flag: demean result?
934 */
convolve_by_ref(float * result,int rlen,float * signal,int siglen,int init,int demean)935 static int convolve_by_ref(float * result, int rlen, float * signal, int siglen,
936 int init, int demean)
937 {
938 int ii, jj, kk, jtop;
939 float val; /* current reference value to apply */
940 double mean; /* for demeaning result */
941
942 if( init ) {
943 for( ii=0 ; ii < rlen ; ii++ ) result[ii] = 0.0;
944 }
945
946 /* TR-locked convolution */
947 for( ii=0 ; ii < refnz ; ii++ ){
948 kk = refin[ii] ; if( kk >= rlen ) break ;
949 val = refts[kk] ;
950
951 /* for each point in the impulse, add its val times irf */
952 /* (top index offset is min(irfdur, length-kk-1)) */
953 jtop = rlen - kk ; if( jtop > siglen ) jtop = siglen ;
954 for( jj=0 ; jj < jtop ; jj++ )
955 result[kk+jj] += val * signal[jj];
956 }
957
958 if( ! demean ) return 0;
959
960 /* demean the result */
961 mean = 0.0;
962 for( ii=0 ; ii < rlen ; ii++ ) mean += result[ii];
963 mean /= rlen;
964 for( ii=0 ; ii < rlen ; ii++ ) result[ii] -= mean;
965
966 return 0;
967 }
968
disp_floats(char * mesg,float * p,int len)969 static int disp_floats(char * mesg, float * p, int len)
970 {
971 int c;
972 if( mesg ) fputs(mesg, stderr);
973 for( c = 0; c < len; c++ ) fprintf(stderr," %f ", p[c]);
974 fprintf(stderr,"\n\n");
975 return 0;
976 }
977
978
979 /* get j, k, t from x, y, sigma, then copy data */
get_signal_from_grid(float * fdest,int length,THD_3dim_dataset * dsrc,float x,float y,float sigma,float scale,int debug)980 static int get_signal_from_grid(float * fdest, int length,
981 THD_3dim_dataset *dsrc, float x, float y,
982 float sigma, float scale, int debug)
983 {
984 float * inptr, * outptr;
985 int nx, ny, nz, nt;
986 int i, j, k, t;
987
988 nx = DSET_NX(dsrc); ny = DSET_NY(dsrc); nz = DSET_NZ(dsrc);
989 nt = DSET_NVALS(dsrc);
990
991 if( nx != length ) {
992 fprintf(stderr, "** FSF: nx, len mis-match, %d != %d\n", nx, length);
993 return 1;
994 }
995
996 /* for j,k, map [-1,1] to [0, ny-1] */
997 j = (int)(0.5+ny*(x+1.0)/2.0); /* x-coord -> j index */
998 if (j < 0 ) j = 0;
999 else if(j >= ny ) j = ny-1;
1000
1001 k = (int)(0.5+nz*(y+1.0)/2.0); /* y-coord -> k index */
1002 if (k < 0 ) k = 0;
1003 else if(k >= nz ) k = nz-1;
1004
1005 /* init to round(nsteps * fraction of max) */
1006 t = (int)(0.5 + genv_sigma_nsteps * sigma / genv_sigma_max);
1007 if ( t < 0 ) t = 0;
1008 else if( t >= nt ) t = nt-1;
1009
1010 if( debug )
1011 fprintf(stderr,"-- fill_array from x=%f, y=%f, s=%f\n"
1012 " at j=%d, k=%d, t=%d\n",
1013 x, y, sigma, j, k, t);
1014
1015 inptr = ((float *)DSET_ARRAY(dsrc, t)) + j*nx + k*nx*ny;
1016 outptr = fdest;
1017 for( i = 0; i < length; i++ ){
1018 *outptr = scale**inptr;
1019 inptr++;
1020 outptr++;
1021 }
1022
1023 return 0;
1024 }
1025
1026
1027 /* as in get_signal_from_grid, but map x,y,s to i,j,k */
inputs_to_coords(THD_3dim_dataset * dset,float x,float y,float sigma)1028 static int inputs_to_coords(THD_3dim_dataset * dset, float x, float y,
1029 float sigma)
1030 {
1031 int nx, ny, nz;
1032 int i, j, k;
1033
1034 nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset);
1035
1036 /* for i,j, map [-1,1] to [0, nx-1] */
1037 /* note, (nx*(x+1.0)/2.0) is in [0,nx], with p(val=nx) very small, */
1038 /* so it should be enough to just limit floor(result) to nx-1 */
1039 i = (int)(nx*(x+1.0)/2.0);
1040 if (i < 0 ) i = 0;
1041 else if(i >= nx ) i = nx-1;
1042
1043 j = (int)(ny*(y+1.0)/2.0);
1044 if (j < 0 ) j = 0;
1045 else if(j >= ny ) j = ny-1;
1046
1047 /* init to round(nsteps * fraction of max) */
1048 k = (int)(genv_sigma_nsteps * sigma / genv_sigma_max);
1049 if ( k < 0 ) k = 0;
1050 else if( k >= nz ) k = nz-1;
1051
1052 fprintf(stderr,"-- fill_array from x=%f, y=%f, s=%f\n"
1053 " at i=%d, j=%d, k=%d\n",
1054 x, y, sigma, i, j, k);
1055
1056 return 0;
1057 }
1058
1059
1060
1061 /* compute the response value directly, given x,y,sigm:
1062 * - compute at a given index, under the assumption that
1063 * such computations will be sequential (at least
1064 * starting from 0
1065 * - at index 0, a new e^x slice will be computed and
1066 * applied across time */
get_signal_computed(float * ts,int tslen,THD_3dim_dataset * dset,float x0,float y0,float sigma,float A,int debug)1067 static int get_signal_computed(float * ts, int tslen, THD_3dim_dataset * dset,
1068 float x0, float y0, float sigma, float A, int debug)
1069 {
1070 static float * sexpgrid = NULL;
1071 static int snxy=0;
1072 float * eptr;
1073 byte * mptr;
1074 double sum;
1075 int nx, ny, tind, sind, nmask;
1076
1077 nx = DSET_NX(dset); ny = DSET_NY(dset);
1078
1079 /* maybe allocate new memory for egrid (hopefully once, "ever") */
1080 if( ! sexpgrid || snxy != nx*ny ) {
1081 if( genv_debug > 1 )
1082 fprintf(stderr,"++ alloc egrid, snxy = %d, nxy = %d\n", snxy, nx*ny);
1083 snxy = nx*ny;
1084 if( sexpgrid ) free(sexpgrid); /* nuke any old copy - maybe never */
1085 sexpgrid = (float *)calloc(snxy, sizeof(float));
1086 if( !sexpgrid ) {
1087 fprintf(stderr,"** PRF egrid alloc failure, nxy = %d\n", snxy);
1088 return 1;
1089 }
1090 }
1091
1092 /* get current e^x grid (scaled to unit area) */
1093 if( compute_e_x_grid(sexpgrid, nx, ny, x0, y0, sigma) ) {
1094 fprintf(stderr,"PRF:FCA: e_x_g failure\n");
1095 return 1;
1096 }
1097
1098 for( tind = 0; tind < tslen; tind++ ) {
1099 mptr = DBLK_ARRAY(dset->dblk, tind);
1100 eptr = sexpgrid;
1101 sum = 0.0;
1102 nmask = 0;
1103 for( sind = 0; sind < nx*ny; sind++ ) {
1104 if( *mptr++ ) {
1105 nmask++;
1106 sum += *eptr;
1107 }
1108 eptr++;
1109 }
1110
1111 if( debug && genv_debug > 2 )
1112 fprintf(stderr,"-- nmask %03d = %d\n",tind,nmask);
1113
1114 ts[tind] = A * sum;
1115 }
1116
1117 /* if requested, write this 2D image (one time only) */
1118 if( genv_gauss_file ) {
1119 char hist[256];
1120 sprintf(hist, "\n == %s\n x = %g, y = %g, sigma = %g\n",
1121 g_model_ver, x0, y0, sigma);
1122 fprintf(stderr, "++ writing PRF model curve to %s\n%s\n",
1123 genv_gauss_file, hist);
1124 write_gauss_file(genv_gauss_file, sexpgrid, nx, ny, hist);
1125
1126 genv_gauss_file = NULL; /* clear - no further writes */
1127 }
1128
1129 return 0;
1130 }
1131
1132 /* ------------------------------------------------------------ */
1133 /* write_gauss_curve */
write_gauss_file(char * fname,float * curve,int nx,int ny,char * hist)1134 static int write_gauss_file(char * fname, float * curve, int nx, int ny,
1135 char * hist)
1136 {
1137 THD_3dim_dataset * dout;
1138 THD_ivec3 inxyz;
1139 THD_fvec3 origin, delta;
1140 float * mptr, * dptr;
1141 byte * bptr;
1142 int ind;
1143
1144 fprintf(stderr,"++ creating gauss dset (%dx%d)", nx, ny);
1145 dout = EDIT_empty_copy(NULL);
1146 LOAD_IVEC3(inxyz, nx, ny, 1); /* nxyz */
1147 origin.xyz[0] = origin.xyz[1] = -1.0; /* origin */
1148 origin.xyz[2] = 0.0;
1149 delta.xyz[0] = delta.xyz[1] = 2.0/(nx-1); /* delta */
1150 delta.xyz[2] = 1.0;
1151
1152 EDIT_dset_items(dout, ADN_nxyz, inxyz,
1153 ADN_xyzorg, origin,
1154 ADN_xyzdel, delta,
1155 ADN_prefix, fname,
1156 ADN_nvals, 2,
1157 ADN_none);
1158
1159 /* first is gaussian, second is first mask (to compare orientations) */
1160 /* use NULL to create space, then copy results */
1161 EDIT_substitute_brick(dout, 0, MRI_float, NULL);
1162 EDIT_substitute_brick(dout, 1, MRI_float, NULL);
1163
1164 /* first copy 'curve' to slice 0 */
1165 mptr = DBLK_ARRAY(dout->dblk, 0);
1166 dptr = curve;
1167 for(ind = 0; ind < nx*ny; ind++)
1168 *mptr++ = *dptr++;
1169
1170 /* now fill mask slice */
1171 mptr = DBLK_ARRAY(dout->dblk, 1);
1172 bptr = DBLK_ARRAY(g_saset->dblk, 0);
1173 for(ind = 0; ind < nx*ny; ind++, bptr++)
1174 *mptr++ = (float)*bptr;
1175
1176 dout->daxes->xxorient = ORI_L2R_TYPE;
1177 dout->daxes->yyorient = ORI_P2A_TYPE;
1178 dout->daxes->zzorient = ORI_I2S_TYPE;
1179
1180 if( hist ) tross_Append_History(dout, hist);
1181
1182 fprintf(stderr,", writing");
1183 DSET_write(dout);
1184
1185 DSET_delete(dout);
1186 fprintf(stderr,", done\n");
1187
1188 return 0;
1189 }
1190
1191
1192 /* - compute a Gaussian curve over the current grid, centered at x0, y0,
1193 * and with the given sigma (unless x0,y0,sigma are the same as previous)
1194 *
1195 * fill with e^-[((x-x0)^2 + (y-y0)^2) / (2*sigma^2)]
1196 */
compute_e_x_grid(float * e,int nx,int ny,float x0,float y0,float sigma)1197 static int compute_e_x_grid(float * e, int nx, int ny, float x0, float y0,
1198 float sigma)
1199 {
1200 float * eptr;
1201 double wscale, sig2, sum;
1202 float xoff, yoff, maxdiff, eval, expval;
1203 int ix, iy, eind;
1204
1205 wscale = 2.0/(nx-1.0); /* scale [0,nx-1] to [0,2] */
1206 sig2 = 0.5/sigma/sigma; /* get 1/(2*sigma^2) */
1207
1208 /* restrict |x-x0|^2 to 2*sigma^2*g_exp_maxval
1209 * so, diff <= sigma*sqrt(2*g_exp_maxval)
1210 */
1211 maxdiff = sigma * sqrt(2*g_exp_maxval);
1212
1213 eptr = e; /* base pointer */
1214 sum = 0.0; /* for scaling to unit integral */
1215 for( iy = 0; iy < ny; iy++ ) {
1216 for( ix = 0; ix < nx; ix++ ) {
1217 xoff = ix*wscale - 1.0f - x0; /* map to [-1,1] and dist from x0 */
1218 yoff = iy*wscale - 1.0f - y0;
1219
1220 /* if we are out of range, skip */
1221 if( fabs(xoff) > maxdiff || fabs(yoff) > maxdiff ) {
1222 *eptr++ = 0.0f;
1223 continue;
1224 }
1225
1226 eval = (xoff*xoff+yoff*yoff)*sig2;
1227 if( eval < 0.0f ) eval = 0.0f; /* just to be sure */
1228
1229 /* could do expval = exp(-eval); */
1230 eind = (int)(eval*g_exp_ipieces); /* truncate towards 0? */
1231 if ( eind < g_exp_nvals ) {
1232 expval = g_exp_ts[eind];
1233 sum += expval;
1234 } else expval = 0.0f;
1235 *eptr++ = expval;
1236 }
1237 }
1238
1239 /* and rescale by 1/sum */
1240 if( sum != 0.0 ) {
1241 sum = 1.0/sum;
1242 for( ix = 0, eptr = e; ix < nx*ny; ix++, eptr++ )
1243 if( *eptr ) *eptr *= sum;
1244 }
1245
1246 return 0;
1247 }
1248
1249