1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*
8 This program estimates the statistical significance levels and statistical
9 power of FMRI datasets by Monte Carlo simulation of random image generation,
10 Gaussian filtering, intensity thresholding, and minimum cluster size
11 thresholding.
12
13 File: AlphaSim.c
14 Author: B. D. Ward
15 Date: 18 June 1997
16
17 Mod: Changed random number generator function from rand to drand48.
18 Date: 26 August 1997
19
20 Mod: Corrected problem: attempt to close a file which was not open.
21 Date: 21 June 1999
22
23 Mod: Corrected problem with count overflow.
24 Date: 30 July 1999
25
26 Mod: Added -mask option. (Adapted from: 3dpc.c)
27 Date: 15 June 2000
28
29 Mod: Added call to AFNI_logger.
30 Date: 15 August 2001
31
32 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
33 Date: 02 December 2002
34
35 Mod: Added -max_clust_size, to override MAX_CLUSTER_SIZE.
36 Also, increased default from 1000 to 10000.
37 Date: 12 Apr 2006 [rickr]
38
39 Mod: Added AFNI_BLUR_* environment variable stuff.
40 Date: 03 Apr 2007 [RWC]
41
42 Mod: Add -fast and -nxyz and -dxyz options.
43 Date: 10 Jan 2008 [RWC]
44
45 Mod: Modify to use OpenMP to parallelize simulations.
46 Date 09 Jun 2009 [RWC]
47
48 Mod: Modify to compute extreme value approximation to Alpha(i).
49 Date: 10 Jun 2009 [RWC]
50
51 */
52
53
54 /*---------------------------------------------------------------------------*/
55
56 #define PROGRAM_NAME "AlphaSim" /* name of this program */
57 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
58 #define PROGRAM_INITIAL "18 Jun 1997" /* date of initial program release */
59 #define PROGRAM_LATEST "09 Jun 2009" /* date of latest program revision */
60
61 /*---------------------------------------------------------------------------*/
62
63 #include "mrilib.h"
64
65 #ifdef USE_OMP
66 #include <omp.h>
67 #endif
68
69 #include <time.h>
70 #include <sys/types.h>
71 #include <unistd.h>
72
73 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
74 #define MAX_CLUSTER_SIZE 10000 /* max. size of cluster for freq. table */
75
76 /*---------------------------------------------------------------------------*/
77 /*
78 Global data
79 */
80
81 static char * mask_filename = NULL; /* file containing the mask */
82 static byte * mask_vol = NULL; /* mask volume */
83 static int mask_ngood = 0; /* number of good voxels in mask volume */
84 /* allow updating via the -max_clust_size option 12 Apr 2006 [rickr] */
85 static int g_max_cluster_size = MAX_CLUSTER_SIZE;
86
87 static int use_zg = 0 ; /* 10 Jan 2008 */
88
89 static unsigned int gseed ; /* global copy of seed */
90 static int gdo_approx = 0 ; /* print out Approx in output table? */
91
92 /*----------------------------------------------------------------------------*/
93 /*! Function to replace use of cdfnor(), which is not thread-safe (not used). */
94
zthresh(double mn,double sd,double pval)95 double zthresh( double mn , double sd , double pval )
96 {
97 double z ;
98
99 if( pval <= 0.0 ) pval = 1.e-15 ;
100 else if( pval >= 1.0 ) pval = 1.0 - 1.e-15 ;
101 z = qginv(pval) ;
102
103 if( sd > 0.0 ) z = sd*z + mn ;
104 return z ;
105 }
106
107 /*---------------------------------------------------------------------------*/
108 /*
109 Routine to display AlphaSim help menu.
110 */
111
display_help_menu()112 void display_help_menu()
113 {
114 printf
115 (
116 "This program performs alpha probability simulations; among other\n"
117 "things, it computes the probability of a random field of noise\n"
118 "producing a cluster of a given size after the noise is thresholded\n"
119 "at a given level ('-pthr').\n"
120 "\n"
121 " *** PLEASE do not use this program any more. Use 3dClustSim! ***\n"
122 "\n"
123 "Usage: \n"
124 "AlphaSim \n"
125 "-nx n1 n1 = number of voxels along x-axis \n"
126 "-ny n2 n2 = number of voxels along y-axis \n"
127 "-nz n3 n3 = number of voxels along z-axis \n"
128 "-dx d1 d1 = voxel size (mm) along x-axis \n"
129 "-dy d2 d2 = voxel size (mm) along y-axis \n"
130 "-dz d3 d3 = voxel size (mm) along z-axis \n"
131 "-nxyz n1 n2 n3 = give all 3 grid dimensions at once \n"
132 "-dxyz d1 d2 d3 = give all 3 voxel sizes at once \n"
133 "[-mask mset] Use the 0 sub-brick of dataset 'mset' as a mask \n"
134 " to indicate which voxels to analyze (a sub-brick \n"
135 " selector is allowed) [default = use all voxels] \n"
136 " Note: The -mask command also REPLACES the \n"
137 " -nx, -ny, -nz, -dx, -dy, and -dz commands, \n"
138 " and takes the volume dimensions from 'mset'. \n"
139 "[-fwhm s] s = Gaussian filter width (FWHM, in mm) \n"
140 "[-fwhmx sx] sx = Gaussian filter width, x-axis (FWHM) \n"
141 "[-fwhmy sy] sy = Gaussian filter width, y-axis (FWHM) \n"
142 "[-fwhmz sz] sz = Gaussian filter width, z-axis (FWHM) \n"
143 "[-sigma s] s = Gaussian filter width (1 sigma, in mm) \n"
144 "[-sigmax sx] sx = Gaussian filter width, x-axis (1 sigma) \n"
145 "[-sigmay sy] sy = Gaussian filter width, y-axis (1 sigma) \n"
146 "[-sigmaz sz] sz = Gaussian filter width, z-axis (1 sigma) \n"
147 "\n"
148 "[-power] perform statistical power calculations \n"
149 "[-ax n1] n1 = extent of active region (in voxels) along x-axis \n"
150 "[-ay n2] n2 = extent of active region (in voxels) along y-axis \n"
151 "[-az n3] n3 = extent of active region (in voxels) along z-axis \n"
152 "[-zsep z] z = z-score separation between signal and noise \n"
153 "\n"
154 "[-rmm r] r = cluster connection radius (mm) \n"
155 " Default is nearest neighbor connection only. \n"
156 "-pthr p p = individual voxel threshold probability \n"
157 "-iter n n = number of Monte Carlo simulations \n"
158 "[-quiet] suppress lengthy per-iteration screen output \n"
159 "[-out file] file = name of output file [default value = screen] \n"
160 "[-max_clust_size size] size = maximum allowed voxels in a cluster \n"
161 "[-seed S] S = random number seed\n"
162 " default seed = 123456789\n"
163 " if seed=0, then program will randomize it\n"
164 "[-fast] Use a faster random number generator:\n"
165 " Can speed program up by about a factor of 2,\n"
166 " but detailed results will differ slightly since\n"
167 " a different sequence of random values will be used.\n"
168 "[-approx] Compute an analytic approximation to the Alpha(i)\n"
169 " result for cluster size i, and print a column of that\n"
170 " value in the output (only if '-power' is NOT used)\n"
171 " ** This analytic approximation is a way to extrapolate\n"
172 " the alpha value for cluster sizes beyond the\n"
173 " reaches of the simulation. The formula for it is\n"
174 " printed above the output table; see the example below.\n"
175 " ** The analytic approximation is only computed if the\n"
176 " table of cluster size vs. alpha is 'large enough'.\n"
177 " ** The approximation formula is of 'extreme value' type,\n"
178 " possibly with an adjustment for smaller i and larger Alpha.\n"
179 "\n"
180 "Unix environment variables you can use:\n"
181 "---------------------------------------\n"
182 " Set AFNI_BLUR_FFT to YES to require blurring be done with FFTs\n"
183 " (the oldest way, and slowest).\n"
184 " Set AFNI_BLUR_FFT to NO and AFNI_BLUR_FIROLD to YES to require\n"
185 " blurring to be done with the old (crude) FIR code (not advised).\n"
186 " If neither of these are set, then blurring is done using the newer\n"
187 " (more accurate) FIR code (recommended).\n"
188 " Results will differ in detail depending on the blurring method\n"
189 " used to generate the simulated noise fields.\n"
190 );
191
192 printf("\n"
193 "SAMPLE OUTPUT:\n"
194 "--------------\n"
195 " AlphaSim -nxyz 64 64 20 -dxyz 3 3 3 -iter 10000 -pthr 0.004 -fwhm 5 \\\n"
196 " -quiet -fast -approx\n"
197 "# Alpha(i) approx 1-exp[-exp(8.720-2.2166*i^0.58-0.05743*posval(12-i)^1.0)]\n"
198 "# Cl Size Frequency CumuProp p/Voxel Max Freq Alpha Approx\n"
199 " 1 1024002 0.584689 0.00414373 0 1.000000 1.000000\n"
200 " 2 358143 0.789183 0.00289373 0 1.000000 1.000000\n"
201 " 3 156346 0.878455 0.00201936 0 1.000000 1.000000\n"
202 " 4 87554 0.928447 0.00144680 0 1.000000 1.000000\n"
203 " 5 48445 0.956108 0.00101929 6 1.000000 1.000000\n"
204 " 6 29126 0.972738 0.00072361 81 0.999400 0.999736\n"
205 " 7 17743 0.982869 0.00051028 407 0.991300 0.992216\n"
206 " 8 11220 0.989276 0.00035867 1082 0.950600 0.948274\n"
207 " 9 6722 0.993114 0.00024910 1453 0.842400 0.844084\n"
208 " 10 4251 0.995541 0.00017525 1564 0.697100 0.697100\n"
209 " 11 2708 0.997087 0.00012336 1426 0.540700 0.543212\n"
210 " 12 1736 0.998079 0.00008700 1132 0.398100 0.407466\n"
211 " 13 1164 0.998743 0.00006157 875 0.284900 0.284900\n"
212 " 14 744 0.999168 0.00004309 615 0.197400 0.195818\n"
213 " 15 485 0.999445 0.00003038 434 0.135900 0.133634\n"
214 " 16 324 0.999630 0.00002150 302 0.092500 0.091099\n"
215 " 17 213 0.999752 0.00001517 196 0.062300 0.062256\n"
216 " 18 140 0.999832 0.00001075 136 0.042700 0.042736\n"
217 " 19 87 0.999881 0.00000767 84 0.029100 0.029499\n"
218 " 20 62 0.999917 0.00000566 61 0.020700 0.020485\n"
219 " 21 49 0.999945 0.00000414 49 0.014600 0.014314\n"
220 " 22 31 0.999962 0.00000289 31 0.009700 0.010064\n"
221 " 23 16 0.999971 0.00000205 16 0.006600 0.007119\n"
222 " 24 10 0.999977 0.00000161 10 0.005000 0.005065\n"
223 " 25 11 0.999983 0.00000131 11 0.004000 0.003624\n"
224 " 26 12 0.999990 0.00000098 12 0.002900 0.002607\n"
225 " 27 3 0.999992 0.00000060 3 0.001700 0.001885\n"
226 " 28 4 0.999994 0.00000050 4 0.001400 0.001370\n"
227 " 29 7 0.999998 0.00000036 7 0.001000 0.001000\n"
228 " 30 1 0.999999 0.00000011 1 0.000300 0.000733\n"
229 " 31 2 1.000000 0.00000008 2 0.000200 0.000540\n"
230 "\n"
231 " That is, thresholded random noise alone (no signal) would produce a cluster\n"
232 " of size 18 or larger about 4.27%% (Alpha) of the time, in a 64x64x20 volume\n"
233 " with cubical 3 mm voxels and a FHWM noise smoothness of 5 mm, and an uncorrected\n"
234 " uncorrected (per voxel) p-value of 0.004 -- this combination of voxel-wise and\n"
235 " cluster-size thresholds would be a logical one to use for a functional map that\n"
236 " had these parameters.\n"
237 "\n"
238 " If you run the exact command above, you will get slightly different results,\n"
239 " due to variations in the random numbers generated in the simulations.\n"
240 "\n"
241 " To plot the approximation on top of the empirical alpha, if the above file\n"
242 " is stored as alp.1D, then the following command can be used:\n"
243 " 1dplot -start 1 -one -ytran 'log(-log(1-a))' alp.1D'[5,6]'\n"
244 " These will plot the log(log) transformed Alpha(i) and the log(log)\n"
245 " transformed approximation together, so you can see how they fit,\n"
246 " especially for the large i and small Alpha cases. Another comparison\n"
247 " technique is to plot the ratio of Approx(i) to Alpha(i):\n"
248 " 1deval -a alp.1D'[5]' -b alp.1D'[6]' -expr 'b/a' | 1dplot -start 1 -stdin\n"
249 " (Since Alpha(i) is always > 0 in the table, there is no division by zero.)\n"
250 "\n"
251 " The analytic approximation formula above uses the function 'posval(x)',\n"
252 " which is defined to be 'max(x,0)' -- this is the correction for small i\n"
253 " (in this example, i < 12). The syntax is compatible with 1deval and 3dcalc.\n"
254 " The breakpoint for the small i/large Alpha correction is set to be at the\n"
255 " cluster size i where Alpha(i) is about 0.3 [in the sample above, 'posval(12-i)'].\n"
256 " For larger i/smaller Alpha, the approximation is of the simple form\n"
257 " Alpha(i) = 1-exp[-exp(a-b*i^p)]\n"
258 " where a, b, p are constants. For a pure extreme value distribution, p=1;\n"
259 " I've found that allowing p < 1 gives slightly better fits in some cases.\n"
260 "\n"
261 "\n"
262 " *** PLEASE do not use this program any more. Use 3dClustSim! ***\n"
263 ) ;
264
265 PRINT_AFNI_OMP_USAGE("AlphaSim",
266 "* OpenMP compilation implies '-fast'\n" ) ;
267
268 exit(0);
269 }
270
271 /*---------------------------------------------------------------------------*/
272 /*
273 Routine to print error message and stop.
274 */
275
276 #if 0
277 void AlphaSim_error (char * message)
278 {
279 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
280 exit(1);
281 }
282 #else
283 # define AlphaSim_error(m) ERROR_exit("AlphaSim: %s",(m))
284 #endif
285
286
287 /*---------------------------------------------------------------------------*/
288 /*
289 Routine to initialize the input options.
290 */
291
initialize_options(int * nx,int * ny,int * nz,float * dx,float * dy,float * dz,int * filter,float * sigmax,float * sigmay,float * sigmaz,int * egfw,int * power,int * ax,int * ay,int * az,float * zsep,float * rmm,float * pthr,int * niter,int * quiet,char ** outfilename,int * seed)292 void initialize_options (
293 int * nx, int * ny, int * nz,
294 float * dx, float * dy, float * dz,
295 int * filter, float * sigmax, float * sigmay, float * sigmaz,
296 int * egfw,
297 int * power, int * ax, int * ay, int * az, float * zsep,
298 float * rmm, float * pthr, int * niter, int * quiet,
299 char ** outfilename, int * seed)
300
301 {
302 *nx = 0; /* number of voxels along x-axis */
303 *ny = 0; /* number of voxels along y-axis */
304 *nz = 0; /* number of voxels along z-axis */
305 *dx = 0.0; /* voxel size along x-axis */
306 *dy = 0.0; /* voxel size along y-axis */
307 *dz = 0.0; /* voxel size along z-axis */
308 *filter = 0; /* flag for Gaussian filtering */
309 *sigmax = 0.0; /* Gaussian filter width, x-axis (1 sigma) */
310 *sigmay = 0.0; /* Gaussian filter width, y-axis (1 sigma) */
311 *sigmaz = 0.0; /* Gaussian filter width, z-axis (1 sigma) */
312 *egfw = 0; /* flag for estimation of filter width */
313 *power = 0; /* flag for power calculations */
314 *ax = 0; /* number of activation voxels along x-axis */
315 *ay = 0; /* number of activation voxels along y-axis */
316 *az = 0; /* number of activation voxels along z-axis */
317 *zsep = 0.0; /* z-score signal and noise separation */
318 *rmm = -1.0 ; /* cluster connection radius (mm) */
319 *pthr = 0.0; /* individual voxel threshold prob. */
320 *niter = 0; /* number of Monte Carlo simulations */
321 *quiet = 0; /* generate screen output (default) */
322 *outfilename = NULL; /* name of output file */
323 *seed = 123456789 ; /* random number seed - user can override */
324 }
325
326
327 /*---------------------------------------------------------------------------*/
328 /*
329 Routine to get user specified input options.
330 */
331
get_options(int argc,char ** argv,int * nx,int * ny,int * nz,float * dx,float * dy,float * dz,int * filter,float * sigmax,float * sigmay,float * sigmaz,int * egfw,int * power,int * ax,int * ay,int * az,float * zsep,float * rmm,float * pthr,int * niter,int * quiet,char ** outfilename,int * seed)332 void get_options (int argc, char ** argv,
333 int * nx, int * ny, int * nz,
334 float * dx, float * dy, float * dz,
335 int * filter, float * sigmax, float * sigmay, float * sigmaz,
336 int * egfw,
337 int * power, int * ax, int * ay, int * az, float * zsep,
338 float * rmm, float * pthr, int * niter, int * quiet,
339 char ** outfilename, int * seed)
340 {
341 int nopt = 1; /* input option argument counter */
342 int ival; /* integer input */
343 float fval; /* float input */
344 int mask_nx=0, mask_ny=0, mask_nz=0, mask_nvox; /* mask dimensions */
345 float mask_dx=0.0, mask_dy=0.0, mask_dz=0.0;
346
347
348 /*----- does user request help menu? -----*/
349 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
350
351 WARNING_message(
352 "*** Please do not use AlphaSim! Use 3dClustSim instead!! ***\n" ) ;
353 ININFO_message("AlphaSim is obsolete, and is no longer kept up-to-date.") ;
354 ININFO_message("In the near future, this program will be disabled.") ;
355
356 /*----- add to program log -----*/
357 AFNI_logger (PROGRAM_NAME,argc,argv);
358
359
360 /*----- initialize the input options -----*/
361 initialize_options (nx, ny, nz, dx, dy, dz, filter, sigmax, sigmay, sigmaz,
362 egfw, power, ax, ay, az, zsep, rmm, pthr, niter, quiet,
363 outfilename, seed);
364
365 /*----- main loop over input options -----*/
366 while (nopt < argc )
367 {
368
369 /*----- -zg [10 Jan 2008] -----*/
370
371 if( strcmp(argv[nopt],"-zg") == 0 || strcmp(argv[nopt],"-fast") == 0 ){
372 use_zg = 1 ; nopt++ ; continue ;
373 }
374 if( strcmp(argv[nopt],"-nozg") == 0 || strcmp(argv[nopt],"-nofast") == 0 ){
375 use_zg = 0 ; nopt++ ; continue ;
376 }
377 if( strcmp(argv[nopt],"-approx") == 0 ){
378 gdo_approx = 1 ; nopt++ ; continue ;
379 }
380
381 /*----- -nxyz n1 n2 n3 [10 Jan 2008: RWC] -----*/
382
383 if( strcmp(argv[nopt],"-nxyz") == 0 ){
384 nopt++ ; if( nopt+2 >= argc ) AlphaSim_error ("need 3 arguments after -nxyz ") ;
385 *nx = (int)strtod(argv[nopt++],NULL); if( *nx <= 0 ) AlphaSim_error("illegal n1 value") ;
386 *ny = (int)strtod(argv[nopt++],NULL); if( *ny <= 0 ) AlphaSim_error("illegal n2 value") ;
387 *nz = (int)strtod(argv[nopt++],NULL); if( *nz <= 0 ) AlphaSim_error("illegal n3 value") ;
388 continue ;
389 }
390
391 /*----- -dxyz d1 d2 d3 [10 Jan 2008: RWC] -----*/
392
393 if( strcmp(argv[nopt],"-dxyz") == 0 ){
394 nopt++ ; if( nopt+2 >= argc ) AlphaSim_error ("need 3 arguments after -dxyz ") ;
395 *dx = strtod(argv[nopt++],NULL); if( *dx <= 0 ) AlphaSim_error("illegal d1 value") ;
396 *dy = strtod(argv[nopt++],NULL); if( *dy <= 0 ) AlphaSim_error("illegal d2 value") ;
397 *dz = strtod(argv[nopt++],NULL); if( *dz <= 0 ) AlphaSim_error("illegal d3 value") ;
398 continue ;
399 }
400
401 /*----- -nx n -----*/
402 if (strncmp(argv[nopt], "-nx", 3) == 0)
403 {
404 nopt++;
405 if (nopt >= argc) AlphaSim_error ("need argument after -nx ");
406 sscanf (argv[nopt], "%d", &ival);
407 if (ival <= 0)
408 AlphaSim_error ("illegal argument after -nx ");
409 *nx = ival;
410 nopt++;
411 continue;
412 }
413
414
415 /*----- -ny n -----*/
416 if (strncmp(argv[nopt], "-ny", 3) == 0)
417 {
418 nopt++;
419 if (nopt >= argc) AlphaSim_error ("need argument after -ny ");
420 sscanf (argv[nopt], "%d", &ival);
421 if (ival <= 0)
422 AlphaSim_error ("illegal argument after -ny ");
423 *ny = ival;
424 nopt++;
425 continue;
426 }
427
428
429 /*----- -nz n -----*/
430 if (strncmp(argv[nopt], "-nz", 3) == 0)
431 {
432 nopt++;
433 if (nopt >= argc) AlphaSim_error ("need argument after -nz ");
434 sscanf (argv[nopt], "%d", &ival);
435 if (ival <= 0)
436 AlphaSim_error ("illegal argument after -nz ");
437 *nz = ival;
438 nopt++;
439 continue;
440 }
441
442
443 /*----- -dx d -----*/
444 if (strncmp(argv[nopt], "-dx", 3) == 0)
445 {
446 nopt++;
447 if (nopt >= argc) AlphaSim_error ("need argument after -dx ");
448 sscanf (argv[nopt], "%f", &fval);
449 if (fval <= 0.0)
450 AlphaSim_error ("illegal argument after -dx ");
451 *dx = fval;
452 nopt++;
453 continue;
454 }
455
456
457 /*----- -dy d -----*/
458 if (strncmp(argv[nopt], "-dy", 3) == 0)
459 {
460 nopt++;
461 if (nopt >= argc) AlphaSim_error ("need argument after -dy ");
462 sscanf (argv[nopt], "%f", &fval);
463 if (fval <= 0.0)
464 AlphaSim_error ("illegal argument after -dy ");
465 *dy = fval;
466 nopt++;
467 continue;
468 }
469
470
471 /*----- -dz d -----*/
472 if (strncmp(argv[nopt], "-dz", 3) == 0)
473 {
474 nopt++;
475 if (nopt >= argc) AlphaSim_error ("need argument after -dz ");
476 sscanf (argv[nopt], "%f", &fval);
477 if (fval <= 0.0)
478 AlphaSim_error ("illegal argument after -dz ");
479 *dz = fval;
480 nopt++;
481 continue;
482 }
483
484
485 /**** -mask mset [14 June 2000] ****/
486
487 if (strcmp(argv[nopt],"-mask") == 0 )
488 {
489 THD_3dim_dataset * mset ; int ii,mc ;
490 nopt++ ;
491 if (nopt >= argc) AlphaSim_error ("need argument after -mask!") ;
492 mask_filename = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
493 if (mask_filename == NULL)
494 AlphaSim_error ("unable to allocate memory");
495 strcpy (mask_filename, argv[nopt]);
496 mset = THD_open_dataset (mask_filename);
497 if (mset == NULL) AlphaSim_error ("can't open -mask dataset!") ;
498
499 mask_vol = THD_makemask( mset , 0 , 1.0,0.0 ) ;
500 if (mask_vol == NULL ) AlphaSim_error ("can't use -mask dataset!") ;
501
502 mask_nvox = DSET_NVOX(mset) ;
503 mask_nx = DSET_NX(mset);
504 mask_ny = DSET_NY(mset);
505 mask_nz = DSET_NZ(mset);
506 mask_dx = DSET_DX(mset);
507 mask_dy = DSET_DY(mset);
508 mask_dz = DSET_DZ(mset);
509
510 DSET_delete(mset) ;
511
512 for( ii=mc=0 ; ii < mask_nvox ; ii++ ) if (mask_vol[ii]) mc++ ;
513 if (mc == 0) AlphaSim_error ("mask is all zeros!") ;
514 if( *quiet < 2 ) printf("++ %d voxels in mask\n",mc) ;
515 mask_ngood = mc;
516 nopt++ ; continue ;
517 }
518
519
520 /*----- -fwhm s -----*/
521 if (strncmp(argv[nopt], "-fwhm", 6) == 0)
522 {
523 nopt++;
524 if (nopt >= argc) AlphaSim_error ("need argument after -fwhm ");
525 sscanf (argv[nopt], "%f", &fval);
526 if (fval < 0.0)
527 AlphaSim_error ("illegal argument after -fwhm ");
528 if (fval > 0.0) *filter = 1;
529 *sigmax = FWHM_TO_SIGMA(fval);
530 *sigmay = FWHM_TO_SIGMA(fval);
531 *sigmaz = FWHM_TO_SIGMA(fval);
532 nopt++;
533 continue;
534 }
535
536
537 /*----- -fwhmx sx -----*/
538 if (strncmp(argv[nopt], "-fwhmx", 6) == 0)
539 {
540 nopt++;
541 if (nopt >= argc) AlphaSim_error ("need argument after -fwhmx ");
542 sscanf (argv[nopt], "%f", &fval);
543 if (fval < 0.0)
544 AlphaSim_error ("illegal argument after -fwhmx ");
545 if (fval > 0.0) *filter = 1;
546 *sigmax = FWHM_TO_SIGMA(fval);
547 nopt++;
548 continue;
549 }
550
551
552 /*----- -fwhmy sy -----*/
553 if (strncmp(argv[nopt], "-fwhmy", 6) == 0)
554 {
555 nopt++;
556 if (nopt >= argc) AlphaSim_error ("need argument after -fwhmy ");
557 sscanf (argv[nopt], "%f", &fval);
558 if (fval < 0.0)
559 AlphaSim_error ("illegal argument after -fwhmy ");
560 if (fval > 0.0) *filter = 1;
561 *sigmay = FWHM_TO_SIGMA(fval);
562 nopt++;
563 continue;
564 }
565
566
567 /*----- -fwhmz sz -----*/
568 if (strncmp(argv[nopt], "-fwhmz", 6) == 0)
569 {
570 nopt++;
571 if (nopt >= argc) AlphaSim_error ("need argument after -fwhmz ");
572 sscanf (argv[nopt], "%f", &fval);
573 if (fval < 0.0)
574 AlphaSim_error ("illegal argument after -fwhmz ");
575 if (fval > 0.0) *filter = 1;
576 *sigmaz = FWHM_TO_SIGMA(fval);
577 nopt++;
578 continue;
579 }
580
581
582 /*----- -sigma s -----*/
583 if (strncmp(argv[nopt], "-sigma", 7) == 0)
584 {
585 nopt++;
586 if (nopt >= argc) AlphaSim_error ("need argument after -sigma ");
587 sscanf (argv[nopt], "%f", &fval);
588 if (fval < 0.0)
589 AlphaSim_error ("illegal argument after -sigma ");
590 if (fval > 0.0) *filter = 1;
591 *sigmax = fval;
592 *sigmay = fval;
593 *sigmaz = fval;
594 nopt++;
595 continue;
596 }
597
598
599 /*----- -sigmax sx -----*/
600 if (strncmp(argv[nopt], "-sigmax", 7) == 0)
601 {
602 nopt++;
603 if (nopt >= argc) AlphaSim_error ("need argument after -sigmax ");
604 sscanf (argv[nopt], "%f", &fval);
605 if (fval < 0.0)
606 AlphaSim_error ("illegal argument after -sigmax ");
607 if (fval > 0.0) *filter = 1;
608 *sigmax = fval;
609 nopt++;
610 continue;
611 }
612
613
614 /*----- -sigmay sy -----*/
615 if (strncmp(argv[nopt], "-sigmay", 7) == 0)
616 {
617 nopt++;
618 if (nopt >= argc) AlphaSim_error ("need argument after -sigmay ");
619 sscanf (argv[nopt], "%f", &fval);
620 if (fval < 0.0)
621 AlphaSim_error ("illegal argument after -sigmay ");
622 if (fval > 0.0) *filter = 1;
623 *sigmay = fval;
624 nopt++;
625 continue;
626 }
627
628
629 /*----- -sigmaz sz -----*/
630 if (strncmp(argv[nopt], "-sigmaz", 7) == 0)
631 {
632 nopt++;
633 if (nopt >= argc) AlphaSim_error ("need argument after -sigmaz ");
634 sscanf (argv[nopt], "%f", &fval);
635 if (fval < 0.0)
636 AlphaSim_error ("illegal argument after -sigmaz ");
637 if (fval > 0.0) *filter = 1;
638 *sigmaz = fval;
639 nopt++;
640 continue;
641 }
642
643
644 /*----- -egfw -----*/
645 if (strncmp(argv[nopt], "-egfw", 5) == 0)
646 {
647 *egfw = 1;
648 nopt++;
649 continue;
650 }
651
652
653 /*----- -power -----*/
654 if (strncmp(argv[nopt], "-power", 6) == 0)
655 {
656 *power = 1;
657 nopt++;
658 continue;
659 }
660
661
662 /*----- -ax n -----*/
663 if (strncmp(argv[nopt], "-ax", 3) == 0)
664 {
665 nopt++;
666 if (nopt >= argc) AlphaSim_error ("need argument after -ax ");
667 sscanf (argv[nopt], "%d", &ival);
668 if (ival <= 0)
669 AlphaSim_error ("illegal argument after -ax ");
670 *ax = ival;
671 nopt++;
672 continue;
673 }
674
675
676 /*----- -ay n -----*/
677 if (strncmp(argv[nopt], "-ay", 3) == 0)
678 {
679 nopt++;
680 if (nopt >= argc) AlphaSim_error ("need argument after -ay ");
681 sscanf (argv[nopt], "%d", &ival);
682 if (ival <= 0)
683 AlphaSim_error ("illegal argument after -ay ");
684 *ay = ival;
685 nopt++;
686 continue;
687 }
688
689
690 /*----- -az n -----*/
691 if (strncmp(argv[nopt], "-az", 3) == 0)
692 {
693 nopt++;
694 if (nopt >= argc) AlphaSim_error ("need argument after -az ");
695 sscanf (argv[nopt], "%d", &ival);
696 if (ival <= 0)
697 AlphaSim_error ("illegal argument after -az ");
698 *az = ival;
699 nopt++;
700 continue;
701 }
702
703
704 /*----- -zsep z -----*/
705 if (strncmp(argv[nopt], "-zsep", 5) == 0)
706 {
707 nopt++;
708 if (nopt >= argc) AlphaSim_error ("need argument after -zsep ");
709 sscanf (argv[nopt], "%f", &fval);
710 if (fval < 0.0)
711 AlphaSim_error ("illegal argument after -zsep ");
712 *zsep = fval;
713 nopt++;
714 continue;
715 }
716
717
718 /*----- -rmm r -----*/
719 if (strncmp(argv[nopt], "-rmm", 4) == 0)
720 {
721 nopt++;
722 if (nopt >= argc) AlphaSim_error ("need argument after -rmm ");
723 sscanf (argv[nopt], "%f", &fval);
724 if (fval <= 0.0) INFO_message("-rmm set to %g ==> NN clustering",fval) ;
725 *rmm = fval;
726 nopt++;
727 continue;
728 }
729
730
731 /*----- -pthr p -----*/
732 if (strncmp(argv[nopt], "-pthr", 5) == 0 || strcmp(argv[nopt],"-pval") == 0 )
733 {
734 nopt++;
735 if (nopt >= argc) AlphaSim_error ("need argument after -pthr ");
736 sscanf (argv[nopt], "%f", &fval);
737 if ((fval <= 0.0) || (fval > 1.0))
738 AlphaSim_error ("illegal argument after -pthr ");
739 *pthr = fval;
740 nopt++;
741 continue;
742 }
743
744
745 /*----- -iter n -----*/
746 if (strncmp(argv[nopt], "-iter", 5) == 0 || strcmp(argv[nopt],"-niter") == 0 )
747 {
748 nopt++;
749 if (nopt >= argc) AlphaSim_error ("need argument after -iter ");
750 sscanf (argv[nopt], "%d", &ival);
751 if (ival <= 0)
752 AlphaSim_error ("illegal argument after -iter ");
753 *niter = ival;
754 nopt++;
755 continue;
756 }
757
758
759 /*----- -out filename -----*/
760 if (strncmp(argv[nopt], "-out", 4) == 0)
761 {
762 nopt++;
763 if (nopt >= argc) AlphaSim_error ("need argument after -out ");
764 *outfilename = AFMALL( char, sizeof(char) * MAX_NAME_LENGTH);
765 strcpy (*outfilename, argv[nopt]);
766 nopt++;
767 continue;
768 }
769
770
771 /*----- -max_clust_size size -----*/
772 if (strncmp(argv[nopt], "-max_clust_size", 10) == 0)
773 {
774 nopt++;
775 if (nopt >= argc)
776 AlphaSim_error ("need argument after -max_clust_size ");
777 g_max_cluster_size = atoi( argv[nopt] );
778 nopt++;
779 continue;
780 }
781
782 /*----- -quiet q -----*/
783 if (strncmp(argv[nopt], "-quiet", 2) == 0)
784 {
785 (*quiet)++;
786 nopt++;
787 continue;
788 }
789
790
791 /*----- -seed S -----*/
792 if (strncmp(argv[nopt], "-seed", 5) == 0)
793 {
794 nopt++;
795 if (nopt >= argc)
796 AlphaSim_error ("need argument after -seed ");
797 *seed = atoi(argv[nopt]);
798 if( *seed == 0 ){
799 *seed = ((int)time(NULL)) + 17*(int)getpid() ;
800 if( *seed < 0 ) *seed = -*seed ;
801 INFO_message("-seed 0 resets to %d",*seed) ;
802 }
803 nopt++;
804 continue;
805 }
806
807
808 /*----- unknown command -----*/
809 ERROR_exit("AlphaSim -- unknown option '%s'",argv[nopt]) ;
810 }
811
812
813 /*----- If mask dataset is used, set dimensions accordingly -----*/
814 if (mask_vol != NULL)
815 {
816 *nx = mask_nx; *ny = mask_ny; *nz = mask_nz;
817 *dx = fabs(mask_dx); *dy = fabs(mask_dy); *dz = fabs(mask_dz);
818 }
819
820 if( *power ) gdo_approx = 0 ;
821
822 }
823
824
825 /*---------------------------------------------------------------------------*/
826 /*
827 Routine to check for valid inputs.
828 */
829
check_for_valid_inputs(int nx,int ny,int nz,float dx,float dy,float dz,int filter,float sigmax,float sigmay,float sigmaz,int power,int ax,int ay,int az,float zsep,float rmm,float pthr,int niter,char * outfilename)830 void check_for_valid_inputs (int nx, int ny, int nz,
831 float dx, float dy, float dz, int filter,
832 float sigmax, float sigmay, float sigmaz,
833 int power, int ax, int ay, int az, float zsep,
834 float rmm, float pthr, int niter,
835 char * outfilename)
836 {
837 FILE * fout = NULL;
838 char message[MAX_NAME_LENGTH]; /* error message */
839
840 if (nx <= 0) AlphaSim_error ("Illegal value for nx ");
841 if (ny <= 0) AlphaSim_error ("Illegal value for ny ");
842 if (nz <= 0) AlphaSim_error ("Illegal value for nz ");
843 if (dx <= 0.0) AlphaSim_error ("Illegal value for dx ");
844 if (dy <= 0.0) AlphaSim_error ("Illegal value for dy ");
845 if (dz <= 0.0) AlphaSim_error ("Illegal value for dz ");
846 if (filter && ((sigmax <= 0.0) || (sigmay <= 0.0) || (sigmaz <= 0.0)))
847 AlphaSim_error ("Illegal Gaussian filter width specification ");
848 if (power && ((ax <= 0) || (ay <= 0) || (az <= 0)))
849 AlphaSim_error ("Illegal dimensions for activation region ");
850 if (power && (zsep <= 0.0)) AlphaSim_error ("Illegal value for zsep ");
851 if ( (rmm < dx) && (rmm < dy) && (rmm < dz) ){
852 rmm = -1.0f ; INFO_message("default NN connectivity being used") ;
853 }
854 if ((pthr <= 0.0) || (pthr >= 1.0))
855 AlphaSim_error ("Illegal value for pthr ");
856 if (niter <= 0) AlphaSim_error ("Illegal value for niter ");
857
858 if (outfilename != NULL && !THD_ok_overwrite())
859 {
860 /*----- see if output file already exists -----*/
861 fout = fopen (outfilename, "r");
862 if (fout != NULL)
863 {
864 fclose(fout) ;
865 sprintf (message, "Output file %s already exists. ", outfilename);
866 AlphaSim_error (message);
867 }
868 }
869
870 /* 12 Apr 2006 [rickr] */
871 /* use a limit of ten million, more than the voxels in a +tlrc dataset */
872 if (g_max_cluster_size < 2 || g_max_cluster_size > 10000000)
873 {
874 sprintf (message, "Invalid -max_clust_size %d. Must be in [%d,%d]. ",
875 g_max_cluster_size, 2, 10000000);
876 AlphaSim_error (message);
877 }
878
879 if( !power ){
880 INFO_message(
881 "For most purposes, program 3dClustSim should be used instead of AlphaSim") ;
882 ININFO_message(
883 "* 3dClustsim lets you run with multiple p-value thresholds at once,") ;
884 ININFO_message(
885 " and will thus be faster than AlphaSim.") ;
886 ININFO_message(
887 "* 3dClustSim will output a table of cluster-size thresholds that can be") ;
888 ININFO_message(
889 " used in AFNI's interactive Clusterize to see cluster alpha levels.") ;
890 }
891
892 return ;
893 }
894
895
896 /*---------------------------------------------------------------------------*/
897 /*
898 Routine to perform program initialization.
899 */
900
initialize(int argc,char ** argv,int * nx,int * ny,int * nz,float * dx,float * dy,float * dz,int * filter,float * sigmax,float * sigmay,float * sigmaz,int * egfw,float * avgsx,float * avgsy,float * avgsz,int * power,int * ax,int * ay,int * az,float * zsep,float * rmm,float * pthr,int * niter,int * quiet,char ** outfilename,long * count,double * sum,double * sumsq,float * power_thr,float ** fim,float ** arfim,long ** freq_table,long ** max_table)901 void initialize (int argc, char ** argv,
902 int * nx, int * ny, int * nz,
903 float * dx, float * dy, float * dz,
904 int * filter, float * sigmax, float * sigmay, float * sigmaz,
905 int * egfw, float * avgsx, float * avgsy, float * avgsz,
906 int * power, int * ax, int * ay, int * az, float * zsep,
907 float * rmm, float * pthr, int * niter, int * quiet,
908 char ** outfilename, long * count,
909 double * sum, double * sumsq, float * power_thr,
910 float ** fim, float ** arfim,
911 long ** freq_table, long ** max_table)
912
913 {
914 int i;
915 int nxyz;
916
917 int which;
918 double p, q, z, mean, sd;
919 int status, seed;
920 double bound;
921
922
923 /*----- get command line inputs -----*/
924 get_options(argc, argv,
925 nx, ny, nz, dx, dy, dz, filter, sigmax, sigmay, sigmaz,
926 egfw, power, ax, ay, az, zsep, rmm, pthr, niter, quiet,
927 outfilename, &seed);
928
929
930 /*----- check for valid inputs -----*/
931 check_for_valid_inputs (*nx, *ny, *nz, *dx, *dy, *dz, *filter,
932 *sigmax, *sigmay, *sigmaz,
933 *power, *ax, *ay, *az, *zsep,
934 *rmm, *pthr, *niter, *outfilename);
935
936
937 /** 09 Jun 2009: for OpenMP, have moved allocation of fim/arfim to main() **/
938
939 /*----- allocate memory space for image data -----*/
940 #if 0
941 nxyz = (*nx) * (*ny) * (*nz);
942 *fim = (float *) malloc(nxyz * sizeof(float));
943 if (*fim == NULL)
944 AlphaSim_error ("memory allocation error");
945 #endif
946
947
948 /*-- if power calculation, allocate memory space for activation region --*/
949 #if 0
950 if (*power)
951 {
952 nxyz = (*ax) * (*ay) * (*az);
953 *arfim = (float *) malloc(nxyz * sizeof(float));
954 if (*arfim == NULL)
955 AlphaSim_error ("memory allocation error");
956 }
957 #endif
958
959 /** 09 Jun 2009: however, the _table variables are common to all threads **/
960
961 /*----- allocate memory space and initialize frequency table -----*/
962 *freq_table = (long *) malloc( g_max_cluster_size * sizeof(long) );
963 if (*freq_table == NULL)
964 AlphaSim_error ("memory allocation error");
965 for (i = 0; i < g_max_cluster_size; i++)
966 (*freq_table)[i] = 0;
967
968
969 /*----- allocate memory space and initialize max cluster size table -----*/
970 *max_table = (long *) malloc( g_max_cluster_size * sizeof(long) );
971 if (*max_table == NULL)
972 AlphaSim_error ("memory allocation error");
973 for (i = 0; i < g_max_cluster_size; i++)
974 (*max_table)[i] = 0;
975
976
977 /*----- initialize voxel intensity sums -----*/
978 *count = 0;
979 *sum = 0.0;
980 *sumsq = 0.0;
981
982
983 /*----- calculate power threshold -----*/
984 if (*power)
985 {
986 which = 2;
987 q = *pthr;
988 p = 1.0 - q;
989 mean = 0.0;
990 sd = 1.0;
991 cdfnor (&which, &p, &q, &z, &mean, &sd, &status, &bound);
992
993 z = *zsep - z;
994 which = 1;
995 cdfnor (&which, &p, &q, &z, &mean, &sd, &status, &bound);
996 *power_thr = p;
997 }
998
999
1000 /*----- initialize ave. est. gaussian filter widths -----*/
1001 if (*egfw)
1002 {
1003 *avgsx = 0.0; *avgsy = 0.0; *avgsz = 0.0;
1004 }
1005
1006 /*----- initialize random number generator -----*/
1007 srand48 (seed); gseed = (unsigned int)seed ;
1008
1009 }
1010
1011 /*---------------------------------------------------------------------------*/
1012 #include "zgaussian.c" /** Ziggurat Gaussian random number generator **/
1013 /*---------------------------------------------------------------------------*/
1014
1015 #ifndef USE_OMP /* these RNGs are not used in the OpenMP code */
1016 /*---------------------------------------------------------------------------*/
1017 /*
1018 Routine to generate a uniform U(0,1) random variate.
1019 */
1020
uniform()1021 float uniform ()
1022 {
1023 return ( (float)drand48() );
1024 }
1025
1026
1027 /*---------------------------------------------------------------------------*/
1028 /*
1029 Routine to generate a normal N(0,1) random variate.
1030 */
1031
normal(float * n1,float * n2)1032 void normal (float * n1, float * n2)
1033 {
1034 float u1, u2;
1035 float r;
1036
1037 /** the -fast way **/
1038 if( use_zg ){ *n1 = zgaussian() ; *n2 = zgaussian() ; return ; }
1039
1040 u1 = 0.0;
1041 while (u1 <= 0.0)
1042 {
1043 u1 = uniform();
1044 }
1045 u2 = uniform();
1046
1047 r = sqrt(-2.0*log(u1));
1048 *n1 = r * cos(2.0*PI*u2);
1049 *n2 = r * sin(2.0*PI*u2);
1050 }
1051 #endif
1052
1053
1054 /*---------------------------------------------------------------------------*/
1055 /*
1056 Routine to calculate the dimensions of the activation region.
1057 */
1058
activation_region(int nx,int ny,int nz,int ax,int ay,int az,int * xbot,int * xtop,int * ybot,int * ytop,int * zbot,int * ztop)1059 void activation_region (int nx, int ny, int nz, int ax, int ay, int az,
1060 int * xbot, int * xtop, int * ybot, int * ytop,
1061 int * zbot, int * ztop)
1062 {
1063 *xbot = nx/2 - ax/2;
1064 *xtop = *xbot + ax;
1065 *ybot = ny/2 - ay/2;
1066 *ytop = *ybot + ay;
1067 *zbot = nz/2 - az/2;
1068 *ztop = *zbot + az;
1069 }
1070
1071
1072 /*---------------------------------------------------------------------------*/
1073 /*
1074 Routine to generate volume of random voxel intensities.
1075 */
1076
generate_image(int nx,int ny,int nz,int power,int ax,int ay,int az,float zsep,float * fim,unsigned short xran[])1077 void generate_image (int nx, int ny, int nz, int power,
1078 int ax, int ay, int az, float zsep, float * fim,
1079 unsigned short xran[] )
1080 {
1081 int nxy, nxyz;
1082 int nxyzdiv2;
1083 int ixyz;
1084 float n1, n2;
1085 int xbot, xtop, ybot, ytop, zbot, ztop;
1086 int ix, jy, kz;
1087
1088
1089 /*----- initialize local variables -----*/
1090 nxy = nx * ny;
1091 nxyz = nxy * nz;
1092 nxyzdiv2 = nxyz / 2;
1093
1094 /*----- generate random image -----*/
1095 #ifndef USE_OMP
1096 for (ixyz = 0; ixyz < nxyzdiv2; ixyz++)
1097 {
1098 normal(&n1, &n2);
1099 fim[ixyz] = n1;
1100 fim[ixyz+nxyzdiv2] = n2;
1101 }
1102 normal(&n1, &n2);
1103 fim[nxyz-1] = n1;
1104 #else
1105 for( ixyz=0 ; ixyz < nxyz ; ixyz++ ) /* OpenMP always uses zgaussian */
1106 fim[ixyz] = zgaussian_sss(xran) ;
1107 #endif
1108
1109 /*----- if power calculation, generate "island" of activation -----*/
1110 if (power)
1111 {
1112 /*--- calculate dimensions of activation region ---*/
1113 activation_region (nx, ny, nz, ax, ay, az,
1114 &xbot, &xtop, &ybot, &ytop, &zbot, &ztop);
1115
1116 /*--- add z-score offset to voxels within activation region ---*/
1117 for (ixyz = 0; ixyz < nxyz; ixyz++)
1118 {
1119 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
1120 if ( (ix >= xbot) && (ix < xtop)
1121 && (jy >= ybot) && (jy < ytop)
1122 && (kz >= zbot) && (kz < ztop) )
1123 fim[ixyz] += zsep;
1124 }
1125 }
1126 }
1127
1128
1129 /*---------------------------------------------------------------------------*/
1130 /*
1131 Routine to apply Gaussian filter to the volume data.
1132 */
1133
gaussian_filter(int nx,int ny,int nz,float dx,float dy,float dz,float rmm,float sigmax,float sigmay,float sigmaz,float * fim)1134 void gaussian_filter (int nx, int ny, int nz, float dx, float dy, float dz,
1135 float rmm, float sigmax, float sigmay, float sigmaz,
1136 float * fim)
1137 {
1138
1139 if( AFNI_yesenv("AFNI_BLUR_FFT") ){
1140 EDIT_blur_allow_fir(0) ; /* 03 Apr 2007 */
1141 EDIT_blur_volume_3d(nx, ny, nz, dx, dy, dz,
1142 MRI_float, fim, sigmax, sigmay, sigmaz);
1143 } else {
1144 FIR_blur_volume_3d( nx,ny,nz , dx,dy,dz , fim , sigmax,sigmay,sigmaz ) ;
1145 }
1146
1147 }
1148
1149
1150 /*---------------------------------------------------------------------------*/
1151 /*
1152 Routine to estimate the Gaussian filter width required to generate the data.
1153 */
1154
estimate_gfw(int nx,int ny,int nz,float dx,float dy,float dz,int niter,int quiet,float * fim,float * avgsx,float * avgsy,float * avgsz)1155 void estimate_gfw (int nx, int ny, int nz, float dx, float dy, float dz,
1156 int niter, int quiet, float * fim,
1157 float * avgsx, float * avgsy, float * avgsz)
1158 {
1159 int nxy, nxyz; /* total number of voxels */
1160 int ixyz; /* voxel index */
1161 int ix, jy, kz, ixyz2;
1162 float fsum, fsq, var;
1163 float dfdx, dfdxsum, dfdxsq, varxx;
1164 float dfdy, dfdysum, dfdysq, varyy;
1165 float dfdz, dfdzsum, dfdzsq, varzz;
1166 int countx, county, countz;
1167 float sx, sy, sz;
1168 float arg;
1169
1170
1171 /*----- initialize local variables -----*/
1172 nxy = nx * ny;
1173 nxyz = nxy * nz;
1174
1175
1176 /*----- estimate the variance of the data -----*/
1177 fsum = 0.0;
1178 fsq = 0.0;
1179 for (ixyz = 0; ixyz < nxyz; ixyz++)
1180 {
1181 fsum += fim[ixyz];
1182 fsq += fim[ixyz] * fim[ixyz];
1183 }
1184 var = (fsq - (fsum * fsum)/nxyz) / (nxyz-1);
1185
1186
1187 /*----- estimate the partial derivatives -----*/
1188 dfdxsum = 0.0; dfdysum = 0.0; dfdzsum = 0.0;
1189 dfdxsq = 0.0; dfdysq = 0.0; dfdzsq = 0.0;
1190 countx = 0; county = 0; countz = 0;
1191 for (ixyz = 0; ixyz < nxyz; ixyz++)
1192 {
1193 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
1194
1195 if (ix+1 < nx)
1196 {
1197 ixyz2 = THREE_TO_IJK (ix+1, jy, kz, nx, nxy);
1198 dfdx = (fim[ixyz2] - fim[ixyz]) / 1.0;
1199 dfdxsum += dfdx;
1200 dfdxsq += dfdx * dfdx;
1201 countx += 1;
1202 }
1203
1204 if (jy+1 < ny)
1205 {
1206 ixyz2 = THREE_TO_IJK (ix, jy+1, kz, nx, nxy);
1207 dfdy = (fim[ixyz2] - fim[ixyz]) / 1.0;
1208 dfdysum += dfdy;
1209 dfdysq += dfdy * dfdy;
1210 county += 1;
1211 }
1212
1213 if (kz+1 < nz)
1214 {
1215 ixyz2 = THREE_TO_IJK (ix, jy, kz+1, nx, nxy);
1216 dfdz = (fim[ixyz2] - fim[ixyz]) / 1.0;
1217 dfdzsum += dfdz;
1218 dfdzsq += dfdz * dfdz;
1219 countz += 1;
1220 }
1221
1222 }
1223
1224 /*----- estimate the variance of the partial derivatives -----*/
1225 if (countx < 2)
1226 varxx = 0.0;
1227 else
1228 varxx = (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1);
1229
1230 if (county < 2)
1231 varyy = 0.0;
1232 else
1233 varyy = (dfdysq - (dfdysum * dfdysum)/county) / (county-1);
1234
1235 if (countz < 2)
1236 varzz = 0.0;
1237 else
1238 varzz = (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1);
1239
1240
1241 /*----- now estimate the equivalent Gaussian filter width -----*/
1242 arg = 1.0 - 0.5*(varxx/var);
1243 if ( (arg <= 0.0) || (varxx == 0.0) )
1244 sx = 0.0;
1245 else
1246 sx = sqrt( -1.0 / (4.0*log(arg)) ) * dx;
1247
1248 arg = 1.0 - 0.5*(varyy/var);
1249 if ( (arg <= 0.0) || (varyy == 0.0) )
1250 sy = 0.0;
1251 else
1252 sy = sqrt( -1.0 / (4.0*log(arg)) ) * dy;
1253
1254 arg = 1.0 - 0.5*(varzz/var);
1255 if ( (arg <= 0.0) || (varzz == 0.0) )
1256 sz = 0.0;
1257 else
1258 sz = sqrt( -1.0 / (4.0*log(arg)) ) * dz;
1259
1260 /*----- save results -----*/
1261 #pragma omp critical (AVG)
1262 {
1263 *avgsx += sx / niter;
1264 *avgsy += sy / niter;
1265 *avgsz += sz / niter;
1266 }
1267
1268 /*----- output results -----*/
1269 if (!quiet){
1270 #pragma omp critical (PRINTF)
1271 {
1272 printf ("var =%f \n", var);
1273 printf ("varxx=%f varyy=%f varzz=%f \n", varxx, varyy, varzz);
1274 printf (" sx=%f sy=%f sz=%f \n", sx, sy, sz);
1275 }
1276 }
1277 }
1278
1279
1280 /*---------------------------------------------------------------------------*/
1281 /*
1282 Routine to copy the activation region into a separate volume.
1283 */
1284
get_activation_region(int nx,int ny,int nz,int ax,int ay,int az,float pthr,float zsep,float * fim,float * arfim)1285 void get_activation_region (int nx, int ny, int nz, int ax, int ay, int az,
1286 float pthr, float zsep, float * fim, float * arfim)
1287 {
1288 int nxy, nxyz;
1289 int axy, axyz;
1290 int ixyz, ixyz2;
1291 int xbot, xtop, ybot, ytop, zbot, ztop;
1292 int ix, jy, kz;
1293
1294
1295 /*----- initialize local variables -----*/
1296 nxy = nx * ny;
1297 nxyz = nxy * nz;
1298 axy = ax * ay;
1299 axyz = axy * az;
1300
1301
1302 /*--- calculate dimensions of activation region ---*/
1303 activation_region (nx, ny, nz, ax, ay, az,
1304 &xbot, &xtop, &ybot, &ytop, &zbot, &ztop);
1305
1306 /*--- copy activation region ---*/
1307 for (ixyz = 0; ixyz < axyz; ixyz++)
1308 {
1309 IJK_TO_THREE (ixyz, ix, jy, kz, ax, axy);
1310 ix += xbot;
1311 jy += ybot;
1312 kz += zbot;
1313 ixyz2 = THREE_TO_IJK (ix, jy, kz, nx, nxy);
1314 arfim[ixyz] = fim[ixyz2];
1315 }
1316
1317
1318 }
1319
1320
1321
1322 /*---------------------------------------------------------------------------*/
1323 /*
1324 Routine to calculate threshold probability.
1325 */
1326
pcalc(int nx,int ny,int nz,float * fim,float zthr)1327 float pcalc (int nx, int ny, int nz, float * fim, float zthr)
1328 {
1329 int nxyz;
1330 int ixyz;
1331 int pcount;
1332 float p;
1333
1334
1335 /*----- initialize local variables -----*/
1336 nxyz = nx * ny * nz;
1337
1338 pcount = 0;
1339 for (ixyz = 0; ixyz < nxyz; ixyz++)
1340 if (fim[ixyz] > zthr) pcount ++;
1341 p = (float)pcount / (float)nxyz;
1342
1343 return (p);
1344 }
1345
1346
1347 /*---------------------------------------------------------------------------*/
1348 /*
1349 Routine to apply threshold to volume data.
1350 */
1351
threshold_data(int nx,int ny,int nz,float * fim,float pthr,long * count,double * sum,double * sumsq,int quiet,int iter)1352 void threshold_data (int nx, int ny, int nz, float * fim,
1353 float pthr, long * count, double * sum, double * sumsq,
1354 int quiet, int iter)
1355 {
1356 int ixyz;
1357 int nxyz;
1358 float zthr;
1359 float pact;
1360
1361 int which;
1362 double p, q, z, mean, sd;
1363 int status;
1364 double bound;
1365 int nzzz=0 ;
1366
1367
1368 /*----- initialize local variables -----*/
1369 nxyz = nx * ny * nz;
1370
1371
1372 /*----- update statistical sums -----*/
1373 if (*count < 1.0e+09)
1374 {
1375 *count += nxyz;
1376 for (ixyz = 0; ixyz < nxyz; ixyz++)
1377 {
1378 *sum += fim[ixyz];
1379 *sumsq += fim[ixyz] * fim[ixyz];
1380 }
1381 }
1382
1383
1384 /*----- calculate z-threshold from statistical sums -----*/
1385 which = 2;
1386 p = 1.0 - pthr;
1387 q = pthr;
1388 mean = (*sum) / (*count);
1389 sd = sqrt(((*sumsq) - ((*sum) * (*sum))/(*count)) / ((*count)-1));
1390
1391 #pragma omp critical (CDFNOR)
1392 cdfnor (&which, &p, &q, &z, &mean, &sd, &status, &bound);
1393 zthr = z;
1394
1395 if (!quiet)
1396 {
1397 #pragma omp critical (PRINTF)
1398 {
1399 pact = pcalc (nx, ny, nz, fim, zthr);
1400 printf ("pthr=%f zthr=%f pact=%f mean=%f sd=%f ", pthr, zthr, pact,mean,sd);
1401 }
1402 }
1403
1404
1405 /*----- apply threshold to image data -----*/
1406 for (ixyz = 0; ixyz < nxyz; ixyz++)
1407 #if 1
1408 if (fim[ixyz] > zthr)
1409 fim[ixyz] = 1.0;
1410 else
1411 fim[ixyz] = 0.0;
1412 #else
1413 if( fim[ixyz] <= zthr ){ fim[ixyz] = 0.0f ; nzzz++; }
1414 #endif
1415
1416 #if 0
1417 #pragma omp critical (PRINTF)
1418 if( nzzz < 0.05f*nxyz ) WARNING_message("nzzz=%d mean=%f sd=%f zthr=%f p=%f q=%f",nzzz,mean,sd,zthr,p,q) ;
1419 #endif
1420
1421 }
1422
1423
1424 /*---------------------------------------------------------------------------*/
1425 /*
1426 Routine to apply mask to volume data.
1427 */
1428
apply_mask(int nx,int ny,int nz,float * fim)1429 void apply_mask (int nx, int ny, int nz, float * fim)
1430
1431 {
1432 int ixyz;
1433 int nxyz;
1434
1435
1436 /*----- initialize local variables -----*/
1437 nxyz = nx * ny * nz;
1438
1439
1440 /*----- apply mask to volume data -----*/
1441 for (ixyz = 0; ixyz < nxyz; ixyz++)
1442 if (! mask_vol[ixyz])
1443 fim[ixyz] = 0.0;
1444
1445
1446 }
1447
1448
1449 /*---------------------------------------------------------------------------*/
1450 /*
1451 Routine to identify clusters.
1452 */
1453
identify_clusters(int nx,int ny,int nz,float dx,float dy,float dz,float rmm,float * fim,int quiet,long * freq_table,long * max_table)1454 void identify_clusters (int nx, int ny, int nz,
1455 float dx, float dy, float dz,
1456 float rmm, float * fim, int quiet,
1457 long * freq_table, long * max_table)
1458 /*
1459 where
1460 rmm = cluster connection radius (mm)
1461 nx = number of voxels along x-axis
1462 ny = number of voxels along y-axis
1463 nz = number of voxels along z-axis
1464 dx = voxel size along x-axis
1465 dy = voxel size along y-axis
1466 dz = voxel size along z-axis
1467 */
1468
1469 {
1470 MCW_cluster_array * clar;
1471 MCW_cluster * cl;
1472 int nxy;
1473 int iclu;
1474 int size, max_size;
1475 int do_save=0 ;
1476
1477
1478 /*----- initialize local variables -----*/
1479 nxy = nx * ny;
1480
1481 /*----- create array of clusters -----*/
1482 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz ,
1483 MRI_float , fim , rmm ) ;
1484
1485 /*----- record cluster sizes -----*/
1486 if ((clar == NULL) || (clar->num_clu == 0))
1487 {
1488 if (!quiet){
1489 #pragma omp critical (PRINTF)
1490 printf ("NumCl=%4d MaxClSz=%4d\n", 0, 0);
1491 }
1492 if (clar != NULL) DESTROY_CLARR(clar);
1493 }
1494 else
1495 {
1496 max_size = 0;
1497 for (iclu = 0; iclu < clar->num_clu; iclu++)
1498 {
1499 cl = clar->clar[iclu] ;
1500 if( cl == NULL ) continue ;
1501
1502 size = cl->num_pt;
1503
1504 if (size < g_max_cluster_size) freq_table[size]++;
1505 else {freq_table[g_max_cluster_size-1]++;
1506 WARNING_message("Cluster size = %d",size); }
1507
1508 if (size > max_size) max_size = size;
1509
1510 }
1511
1512 if (max_size < g_max_cluster_size) max_table[max_size]++;
1513 else max_table[g_max_cluster_size-1]++;
1514
1515 if (!quiet){
1516 #pragma omp critical (PRINTF)
1517 printf ("NumCl=%4d MaxClSz=%4d\n", clar->num_clu, max_size);
1518 }
1519
1520 #if 0
1521 if( do_save ){
1522 for (iclu = 0; iclu < clar->num_clu; iclu++){
1523 cl = clar->clar[iclu] ; if( cl == NULL ) continue ;
1524 MCW_cluster_to_vol( nx,ny,nz , MRI_float , fim , cl ) ;
1525 }
1526 #pragma omp critical (DO_SAVE)
1527 { char fname[32] ; FILE *fp ;
1528 for( iclu=0 ; ; iclu++ ){
1529 sprintf(fname,"Clu%06d",iclu) ;
1530 if( !THD_is_file(fname) ) break ;
1531 }
1532 fp = fopen(fname,"w") ;
1533 fwrite( fim , sizeof(float) , nx*ny*nz , fp ) ;
1534 fclose(fp) ;
1535 WARNING_message("Wrote file %s",fname) ;
1536 }
1537 }
1538 #endif
1539
1540 DESTROY_CLARR(clar);
1541 }
1542
1543 }
1544
1545
1546 /*---------------------------------------------------------------------------*/
1547 /*
1548 Routine to generate requested output.
1549 */
1550
output_results(int nx,int ny,int nz,float dx,float dy,float dz,int filter,float sigmax,float sigmay,float sigmaz,int egfw,float avgsx,float avgsy,float avgsz,int power,int ax,int ay,int az,float zsep,float rmm,float pthr,int niter,char * outfilename,long * freq_table,long * max_table,int quiet)1551 void output_results (int nx, int ny, int nz, float dx, float dy, float dz,
1552 int filter, float sigmax, float sigmay, float sigmaz,
1553 int egfw, float avgsx, float avgsy, float avgsz,
1554 int power, int ax, int ay, int az, float zsep,
1555 float rmm, float pthr, int niter, char * outfilename,
1556 long * freq_table, long * max_table, int quiet)
1557 {
1558 const float EPSILON = 1.0e-6;
1559 int i, j;
1560 float divisor;
1561 float * prob_table;
1562 float * alpha_table;
1563 float * cum_prop_table;
1564 long total_num_clusters;
1565 char message[MAX_NAME_LENGTH]; /* error message */
1566 FILE * fout=NULL;
1567
1568 float afit=0.0f , bfit=0.0f , cfit=0.0f , cpow=-1.0f , ipow=1.0f , val ; /* 10 Jun 2009 */
1569 int ibot , ihigh=0 , ilow , itop , ndim ;
1570
1571 /*----- allocate memory space for probability table -----*/
1572 prob_table = (float *) malloc( g_max_cluster_size * sizeof(float) );
1573 if (prob_table == NULL)
1574 AlphaSim_error ("memory allocation error");
1575 for (i = 1; i < g_max_cluster_size; i++)
1576 prob_table[i] = 0.0;
1577
1578 /*----- allocate memory space for alpha table -----*/
1579 alpha_table = (float *) malloc( g_max_cluster_size * sizeof(float) );
1580 if (alpha_table == NULL)
1581 AlphaSim_error ("memory allocation error");
1582 for (i = 1; i < g_max_cluster_size; i++)
1583 alpha_table[i] = 0.0;
1584
1585 /*----- allocate memory space for cum. prop. of cluster size table -----*/
1586 cum_prop_table = (float *) malloc( g_max_cluster_size * sizeof(float) );
1587 if (cum_prop_table == NULL)
1588 AlphaSim_error ("memory allocation error");
1589 for (i = 1; i < g_max_cluster_size; i++)
1590 cum_prop_table[i] = 0.0;
1591
1592 total_num_clusters = 0;
1593 for (i = 1; i < g_max_cluster_size; i++)
1594 total_num_clusters += freq_table[i];
1595
1596 if (power)
1597 divisor = (float)(niter) * ax * ay * az;
1598 else
1599 if (mask_vol)
1600 divisor = (float)(niter) * mask_ngood;
1601 else
1602 divisor = (float)(niter) * nx * ny * nz;
1603
1604 for (i = 1; i < g_max_cluster_size; i++)
1605 {
1606 prob_table[i] = i * freq_table[i] / divisor;
1607 alpha_table[i] = (float)max_table[i] / (float)niter;
1608 cum_prop_table[i] = (float)freq_table[i] / (float)total_num_clusters;
1609 }
1610
1611 for (i = 1; i < g_max_cluster_size-1; i++)
1612 {
1613 j = g_max_cluster_size - i;
1614 prob_table[j-1] += prob_table[j];
1615 alpha_table[j-1] += alpha_table[j];
1616 cum_prop_table[i+1] += cum_prop_table[i];
1617 }
1618
1619
1620 /*----- if output file has not been specified, use stdout -----*/
1621 if (outfilename == NULL)
1622 fout = stdout;
1623 else
1624 {
1625
1626 if (!THD_ok_overwrite()) {
1627 /*----- see if output file already exists -----*/
1628 fout = fopen (outfilename, "r");
1629 if (fout != NULL)
1630 {
1631 fclose(fout) ;
1632 sprintf (message, "file %s already exists. ", outfilename);
1633 AlphaSim_error (message);
1634 }
1635 }
1636
1637 /*----- open file for output -----*/
1638 fout = fopen (outfilename, "w");
1639 if (fout == NULL)
1640 {
1641 AlphaSim_error ("unable to open output file ");
1642 }
1643 }
1644
1645 /*----- print out the results -----*/
1646 if(quiet<2)fprintf (fout, "\n\n");
1647 if(quiet<2)fprintf (fout, "# Program: %s \n", PROGRAM_NAME);
1648 if(quiet<2)fprintf (fout, "# Author: %s \n", PROGRAM_AUTHOR);
1649 if(quiet<2)fprintf (fout, "# Initial Release: %s \n", PROGRAM_INITIAL);
1650 if(quiet<2)fprintf (fout, "# Latest Revision: %s \n", PROGRAM_LATEST);
1651 if(quiet<2)fprintf (fout, "\n");
1652
1653 if(quiet<2)fprintf (fout, "# Data set dimensions: \n");
1654 if(quiet<2)fprintf (fout, "# nx = %5d ny = %5d nz = %5d (voxels)\n", nx, ny, nz);
1655 if(quiet<2)fprintf (fout, "# dx = %5.2f dy = %5.2f dz = %5.2f (mm)\n", dx, dy, dz);
1656
1657 if (mask_vol)
1658 if(quiet<2)fprintf (fout, "\n# Mask filename = %s \n", mask_filename);
1659 if (mask_vol && !power)
1660 if(quiet<2)fprintf (fout, "# Voxels in mask = %5d \n", mask_ngood);
1661
1662 if(quiet<2)fprintf (fout, "\n# Gaussian filter widths: \n");
1663 if(quiet<2)fprintf (fout, "# sigmax = %5.2f FWHMx = %5.2f \n",
1664 sigmax, sigmax * 2.0*sqrt(2.0*log(2.0)));
1665 if(quiet<2)fprintf (fout, "# sigmay = %5.2f FWHMy = %5.2f \n",
1666 sigmay, sigmay * 2.0*sqrt(2.0*log(2.0)));
1667 if(quiet<2)fprintf (fout, "# sigmaz = %5.2f FWHMz = %5.2f \n\n",
1668 sigmaz, sigmaz * 2.0*sqrt(2.0*log(2.0)));
1669
1670 if (egfw)
1671 {
1672 if(quiet<2)fprintf (fout, "# Estimated Gaussian filter widths: \n");
1673 if(quiet<2)fprintf (fout, "# Ave sx = %f Ave sy = %f Ave sz = %f \n\n",
1674 avgsx, avgsy, avgsz);
1675 }
1676
1677 if (power)
1678 {
1679 if(quiet<2)fprintf (fout, "# Activation Region for Power Calculations: \n");
1680 if(quiet<2)fprintf (fout, "# ax = %5d ay = %5d az = %5d (voxels) \n",
1681 ax, ay, az);
1682 if(quiet<2)fprintf (fout, "# z separation = %f \n\n", zsep);
1683 }
1684
1685 if(quiet<2){
1686 if( rmm > 0.0f )
1687 fprintf (fout, "# Cluster connection radius: rmm = %5.2f \n\n", rmm);
1688 else
1689 fprintf (fout, "# Cluster connection = Nearest Neighbor\n") ;
1690 }
1691 if(quiet<2)fprintf (fout, "# Threshold probability: pthr = %e \n\n", pthr);
1692 if(quiet<2)fprintf (fout, "# Number of Monte Carlo iterations = %5d \n\n", niter);
1693
1694 /** estimate alpha[i] as an extreme value distribution [10 Jun 2009] **/
1695
1696 if( !power ){
1697 #define NIPOW 31
1698 #define DIPOW 0.02f
1699 int ii ;
1700 float *zvec=NULL , yvec[4] , *rvec[3] ;
1701 float *ivv[NIPOW] , ipp[NIPOW] ; int kpp ;
1702 float atop=0.98f, ahigh=0.3f, alow=0.03f, abot=4.0f/niter ;
1703 float cout , pp , cbest ;
1704
1705 for( ibot=1 ; /* find first value of alpha <= atop */
1706 ibot < g_max_cluster_size && alpha_table[ibot] > atop ;
1707 ibot++ ) ; /*nada*/
1708 if( ibot == g_max_cluster_size ) goto EXTREME_DONE ;
1709
1710 for( ihigh=ibot+1 ; /* find last value of alpha >= ahigh */
1711 ihigh < g_max_cluster_size && alpha_table[ihigh] > ahigh ;
1712 ihigh++ ) ; /*nada*/
1713 ihigh-- ;
1714 if( ihigh-ibot < 4 ) goto EXTREME_DONE ;
1715
1716 for( ilow=ihigh+1 ; /* find last value of alpha >= alow */
1717 ilow < g_max_cluster_size && alpha_table[ilow] > alow ;
1718 ilow++ ) ; /*nada*/
1719 ilow-- ;
1720 if( ilow-ihigh < 4 ) goto EXTREME_DONE ;
1721
1722 for( itop=ilow+1 ; /* find last value of alpha >= abot */
1723 itop < g_max_cluster_size && alpha_table[itop] > abot ;
1724 itop++ ) ; /*nada*/
1725 itop-- ;
1726 if( itop-ilow < 4 ) goto EXTREME_DONE ;
1727 ndim = itop-ibot+1 ;
1728
1729 /** setup different powers of i^ipow where ipow is stored in ipp[] */
1730
1731 for( kpp=0 ; kpp < NIPOW ; kpp++ ){
1732 ipp[kpp] = 1.0f - kpp*DIPOW ;
1733 ivv[kpp] = (float *)malloc(sizeof(float)*ndim) ;
1734 for( i=0 ; i < ndim ; i++ ){
1735 val = (float)(i+ibot) ; ivv[kpp][i] = -powf(val,ipp[kpp]) ;
1736 }
1737 }
1738
1739 /** setup other regressors and data **/
1740
1741 rvec[0] = (float *)malloc(sizeof(float)*ndim) ; /* regressor = 1 */
1742 rvec[2] = (float *)malloc(sizeof(float)*ndim) ; /* posfunc(ihigh-i)^cpow */
1743 zvec = (float *)malloc(sizeof(float)*ndim) ; /* 'data' */
1744 for( i=0 ; i < ndim ; i++ ){
1745 ii = i+ibot ;
1746 val = (float)ii ;
1747 rvec[0][i] = 1.0f ;
1748 zvec[i] = logf( -logf( 1.0f - alpha_table[ii] ) ) ;
1749 }
1750
1751 cbest = 1.e+38 ; cpow = -1.0f ;
1752
1753 /** do fits without the posfunc() component, keep best one **/
1754
1755 for( kpp=0 ; kpp < NIPOW ; kpp++ ){
1756 yvec[0] = 0.0f ; yvec[1] = 1.0f ; rvec[1] = ivv[kpp] ;
1757 cout = cl1_solve( ndim , 2 , zvec , rvec , yvec , 1 ) ;
1758 if( cout >= 0.0f && cout < cbest ){
1759 cbest = cout; cpow = 1.0f; afit = yvec[0]; bfit = yvec[1]; cfit = 0.0f; ipow = ipp[kpp];
1760 }
1761 }
1762
1763 cbest *= 0.90f ; /* bias results towards no small i correction */
1764
1765 /** now loop over powers pp for correction for small i values = posfunc(ihigh-i)^pp **/
1766
1767 for( pp=1.0f ; pp <= 2.501f ; pp+=0.1f ){
1768 for( i=0 ; i < ndim ; i++ ){ /* setup rvec[2] regressor */
1769 ii = i+ibot ;
1770 val = (float)(ihigh-ii) ;
1771 if( val <= 0.0f ) val = 0.0f ;
1772 else val = powf(val,pp) ;
1773 rvec[2][i] = val ;
1774 }
1775
1776 for( kpp=0 ; kpp < NIPOW ; kpp++ ){ /* loop over powers for the large i fit */
1777 yvec[0] = 0.0f ; yvec[1] = 1.0f ; yvec[2] = 0.0f ; rvec[1] = ivv[kpp] ;
1778 cout = cl1_solve( ndim , 3 , zvec , rvec , yvec , 1 ) ;
1779 if( cout >= 0.0f && cout < cbest ){
1780 cbest = cout; cpow = pp; afit = yvec[0]; bfit = yvec[1]; cfit = yvec[2]; ipow = ipp[kpp];
1781 }
1782 } /* of loop over ipow power */
1783 } /* end of loop over cpow=pp power */
1784
1785 /** print Approx formula **/
1786
1787 if( cfit != 0.0f ) /* two-part fit was best */
1788 printf(
1789 "# Alpha(i) approx 1-exp[-exp"
1790 "(%.3f-%.4f*i^%.2f%+.4g*posval(%d-i)^%.1f)]\n" ,
1791 afit , bfit , ipow , cfit , ihigh , cpow ) ;
1792 else /* simpler fit was best */
1793 printf(
1794 "# Alpha(i) approx 1-exp[-exp(%.3f-%.4f*i^%.2f)]\n" ,
1795 afit , bfit , ipow ) ;
1796
1797 EXTREME_DONE: /** toss the trash **/
1798 if( zvec != NULL ){
1799 free(zvec); free(rvec[0]); free(rvec[2]);
1800 for( kpp=0 ; kpp < NIPOW ; kpp++ ) free(ivv[kpp]) ;
1801 }
1802 }
1803
1804 /** print the tabular output **/
1805
1806 if( quiet < 2 ){
1807 if (!power){
1808 fprintf (fout, "# Cl Size Frequency CumuProp p/Voxel"
1809 " Max Freq Alpha");
1810 if( gdo_approx && bfit > 0.0f && cpow > 0.0f ) fprintf(fout , " Approx") ;
1811 fprintf(fout,"\n") ;
1812 }
1813 else {
1814 fprintf (fout, "# Cl Size Frequency CumuProp p/Voxel"
1815 " Max Freq Power\n");
1816 }
1817 }
1818
1819 for (i = 1; i < g_max_cluster_size; i++) {
1820 if (alpha_table[i] < EPSILON)
1821 break;
1822 else {
1823 fprintf (fout, "%7d %12ld %10.6f %10.8f %7ld %10.6f",
1824 i, freq_table[i], cum_prop_table[i], prob_table[i],
1825 max_table[i], alpha_table[i]);
1826 if( gdo_approx && bfit > 0.0f && cpow > 0.0f ){
1827 val = (float)(ihigh-i) ;
1828 if( val < 0.0f ) val = 0.0f ;
1829 else val = powf(val,cpow) ;
1830 val = afit - bfit*powf((float)i,ipow) + cfit * val ;
1831 val = 1.0f - expf( -expf(val) ) ;
1832 fprintf(fout,"%10.6f",val) ;
1833 }
1834 fprintf(fout,"\n") ;
1835 }
1836 }
1837
1838 if( fout != stdout ) fclose(fout);
1839
1840 return ;
1841 }
1842
1843
1844 /*---------------------------------------------------------------------------*/
1845 /*
1846 Routine to terminate program.
1847 */
1848
terminate(float ** fim,float ** arfim,long ** freq_table,long ** max_table)1849 void terminate (float ** fim, float ** arfim,
1850 long ** freq_table, long ** max_table)
1851 {
1852 #if 0
1853 if (*fim != NULL)
1854 { free (*fim); *fim = NULL; }
1855 #endif
1856
1857 #if 0
1858 if (*arfim != NULL)
1859 { free (*arfim); *arfim = NULL; }
1860 #endif
1861
1862 if (*freq_table != NULL)
1863 { free (*freq_table); *freq_table = NULL; }
1864
1865 if (*max_table != NULL)
1866 { free (*max_table); *max_table = NULL; }
1867 }
1868
1869
1870 /*---------------------------------------------------------------------------*/
1871 /*
1872 Alpha simulation.
1873 */
1874
main(int argc,char ** argv)1875 int main (int argc, char ** argv)
1876 {
1877 int nx; /* number of voxels along x-axis */
1878 int ny; /* number of voxels along y-axis */
1879 int nz; /* number of voxels along z-axis */
1880 float dx; /* voxel size along x-axis */
1881 float dy; /* voxel size along y-axis */
1882 float dz; /* voxel size along z-axis */
1883 int filter; /* flag for Gaussian filtering */
1884 float sigmax; /* Gaussian filter width, x-axis (1 sigma) */
1885 float sigmay; /* Gaussian filter width, y-axis (1 sigma) */
1886 float sigmaz; /* Gaussian filter width, z-axis (1 sigma) */
1887 int egfw; /* flag for estimation of filter width */
1888 float avgsx; /* est. Gaussian filter width, x-axis (1 sigma) */
1889 float avgsy; /* est. Gaussian filter width, x-axis (1 sigma) */
1890 float avgsz; /* est. Gaussian filter width, x-axis (1 sigma) */
1891 int power; /* flag for perform power calculations */
1892 int ax; /* number of activation voxels along x-axis */
1893 int ay; /* number of activation voxels along y-axis */
1894 int az; /* number of activation voxels along z-axis */
1895 float zsep; /* z-score separation between signal and noise */
1896 float rmm; /* cluster connection radius (mm) */
1897 float pthr; /* individual voxel threshold probability */
1898 int niter; /* number of Monte Carlo simulations */
1899 int quiet; /* set to 1 to suppress screen output */
1900 char *outfilename; /* name of output file */
1901
1902 long count;
1903 double sum, sumsq;
1904 float power_thr;
1905
1906 float *fim = NULL; /* won't be used in OpenMP version */
1907 float *arfim = NULL; /* won't be used in OpenMP version */
1908 long *freq_table = NULL;
1909 long *max_table = NULL;
1910 #ifdef USE_OMP
1911 long **mtab=NULL , **ftab=NULL ; int nthr=1 ; /* arrays of tables */
1912 #endif
1913
1914
1915 /*----- Identify software -----*/
1916 #if 0
1917 printf ("\n\n");
1918 printf ("Program: %s \n", PROGRAM_NAME);
1919 printf ("Author: %s \n", PROGRAM_AUTHOR);
1920 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
1921 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
1922 printf ("\n");
1923 #endif
1924
1925 PRINT_VERSION("AlphaSim") ; AUTHOR(PROGRAM_AUTHOR) ;
1926 mainENTRY("AlphaSim main") ; machdep() ;
1927
1928 /*----- program initialization -----*/
1929 initialize (argc, argv,
1930 &nx, &ny, &nz, &dx, &dy, &dz, &filter, &sigmax, &sigmay, &sigmaz,
1931 &egfw, &avgsx, &avgsy, &avgsz, &power, &ax, &ay, &az, &zsep,
1932 &rmm, &pthr, &niter, &quiet, &outfilename, &count, &sum, &sumsq,
1933 &power_thr, &fim, &arfim, &freq_table, &max_table);
1934
1935
1936 AFNI_OMP_START ;
1937 #pragma omp parallel if( niter > 99 )
1938 {
1939 int iter , qqq=quiet ; float *fim , *arfim=NULL ;
1940 long count=0; double sum=0.0, sumsq=0.0 ;
1941 long *mt , *ft ; int ithr=0 ; unsigned short xran[3] ;
1942
1943 /* create separate tables for each thread, if using OpenMP */
1944 #ifdef USE_OMP
1945 ithr = omp_get_thread_num() ;
1946 #pragma omp master /* only in the master thread */
1947 {
1948 nthr = omp_get_num_threads() ;
1949 mtab = (long **)malloc(sizeof(long *)*nthr) ; /* arrays of tables */
1950 ftab = (long **)malloc(sizeof(long *)*nthr) ;
1951 INFO_message("Using %d OpenMP threads",nthr) ;
1952 }
1953 #pragma omp barrier /* all threads wait until the above is finished */
1954 /* create tables for each thread separately */
1955 mtab[ithr] = mt = (long *) calloc( g_max_cluster_size , sizeof(long) ) ;
1956 ftab[ithr] = ft = (long *) calloc( g_max_cluster_size , sizeof(long) ) ;
1957
1958 /* initialize random seed array for each thread separately */
1959 xran[2] = ( gseed & 0xffff) + (unsigned short)ithr ;
1960 xran[1] = ((gseed >> 16) & 0xffff) - (unsigned short)ithr ;
1961 xran[0] = 0x330e + (unsigned short)ithr ;
1962
1963 #else /* not OpenMP ==> only one set of tables */
1964 mt = max_table ;
1965 ft = freq_table ;
1966 #endif
1967
1968 /** malloc of image space local to each thread [09 Jun 2009] **/
1969
1970 fim = (float *)malloc(sizeof(float)*nx*ny*nz) ;
1971 if( power )
1972 arfim = (float *)malloc(sizeof(float)*ax*ay*az) ;
1973
1974 /*----- Monte Carlo iterations -----*/
1975 #pragma omp for
1976 for (iter = 1; iter <= niter; iter++)
1977 {
1978 if (!qqq){
1979 #pragma omp critical (PRINTF)
1980 printf ("Iter =%5d \n", iter);
1981 }
1982
1983 /*----- generate volume of random voxel intensities -----*/
1984 generate_image (nx, ny, nz, power, ax, ay, az, zsep, fim , xran );
1985
1986
1987 /*----- apply gaussian filter to volume data -----*/
1988 if (filter) gaussian_filter (nx, ny, nz, dx, dy, dz, rmm,
1989 sigmax, sigmay, sigmaz, fim);
1990
1991
1992 /*----- estimate equivalent gaussian filter width -----*/
1993 if (egfw) estimate_gfw (nx, ny, nz, dx, dy, dz,
1994 niter, qqq, fim, &avgsx, &avgsy, &avgsz);
1995
1996
1997 /*----- if power calculation, get volume corresponding to -----*/
1998 /*----- activation region and corresponding power threshold -----*/
1999 if (power) get_activation_region (nx, ny, nz, ax, ay, az, pthr, zsep,
2000 fim, arfim);
2001
2002
2003 /*----- apply threshold to volume data -----*/
2004 if (power) threshold_data (ax, ay, az,
2005 arfim, power_thr, &count, &sum, &sumsq,
2006 qqq, iter);
2007 else
2008 threshold_data (nx, ny, nz,
2009 fim, pthr, &count, &sum, &sumsq,
2010 qqq, iter);
2011
2012
2013 /*----- apply mask to volume data -----*/
2014 if (mask_vol && (!power)) apply_mask (nx, ny, nz, fim);
2015
2016
2017 /*----- identify clusters, add to tables ft[] and mt[] -----*/
2018 if (power)
2019 identify_clusters (ax, ay, az, dx, dy, dz, rmm, arfim, qqq,
2020 ft, mt);
2021 else
2022 identify_clusters (nx, ny, nz, dx, dy, dz, rmm, fim, qqq,
2023 ft, mt);
2024
2025 } /* end of long iteration loop */
2026
2027 if( arfim != NULL ) free(arfim) ; /* toss the local trash */
2028 free(fim) ;
2029
2030 } /* end OpenMP parallelization */
2031 AFNI_OMP_END ;
2032
2033 #ifdef USE_OMP /* sum tables from various threads into one result */
2034 if( nthr == 1 ){
2035 memcpy(freq_table,ftab[0],sizeof(long)*g_max_cluster_size) ;
2036 memcpy(max_table ,mtab[0],sizeof(long)*g_max_cluster_size) ;
2037 } else { /* note this is outside of OpenMP! */
2038 int ithr , ii ; long *ft , *mt ;
2039 for( ithr=0 ; ithr < nthr ; ithr++ ){
2040 ft = ftab[ithr] ; mt = mtab[ithr] ;
2041 for( ii=0 ; ii < g_max_cluster_size ; ii++ ){
2042 freq_table[ii] += ft[ii] ; max_table[ii] += mt[ii] ;
2043 }
2044 }
2045 }
2046 #endif
2047
2048 /*----- generate requested output -----*/
2049 output_results (nx, ny, nz, dx, dy, dz, filter, sigmax, sigmay, sigmaz,
2050 egfw, avgsx, avgsy, avgsz, power, ax, ay, az, zsep,
2051 rmm, pthr, niter, outfilename, freq_table, max_table,quiet);
2052
2053
2054 /*----- terminate program -----*/
2055 terminate (&fim, &arfim, &freq_table, &max_table);
2056
2057 exit(0);
2058 }
2059