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 /*** Adapted from AlphaSim.c ***/
8
9 #include "mrilib.h"
10
11 #ifdef USE_OMP
12 #include <omp.h>
13 #endif
14
15 /**** C source codes to include, for compilation efficency (we hope & pray) ***/
16 /*----------------------------------------------------------------------------*/
17 #include "zgaussian.c" /** Ziggurat Gaussian random number generator **/
18 /*----------------------------------------------------------------------------*/
19 #include "mri_radial_random_field.c" /** 3D FFT-based random field generator **/
20 /*----------------------------------------------------------------------------*/
21 #ifndef IND_T
22 typedef unsigned short ind_t ;
23 # define IND_T unsigned short
24 # undef USE_UBYTE
25 # define MAX_IND 65535u
26 #endif
27 #include "thd_Xdataset.c" /** reading .sdat datasets **/
28 /*----------------------------------------------------------------------------*/
29
30 #include <time.h>
31 #include <sys/types.h>
32 #include <unistd.h>
33
34 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
35 #define MAX_CLUSTER_SIZE 99999 /* max. size of cluster for freq. table */
36
37 #undef ALLOW_LOHI
38 #ifndef ALLOW_LOHI
39 # define do_lohi 0
40 #else
41 static int do_lohi = 0 ;
42 #endif
43
44 /*---------------------------------------------------------------------------*/
45 /*
46 Global data
47 */
48
49 static THD_3dim_dataset *mask_dset = NULL ; /* mask dataset */
50 static byte *mask_vol = NULL; /* mask volume */
51 static int mask_nvox = 0, mask_ngood = 0; /* number of good voxels in mask volume */
52
53 static int const max_cluster_size = MAX_CLUSTER_SIZE ;
54
55 static int nx = 64 ;
56 static int ny = 64 ;
57 static int nz = 32 ;
58 static int nxy ;
59 static int nxyz ;
60 static int nxyz1 ;
61 static float dx = 3.5f ;
62 static float dy = 3.5f ;
63 static float dz = 3.5f ;
64 static float fwhm_x = 0.0f ;
65 static float fwhm_y = 0.0f ;
66 static float fwhm_z = 0.0f ;
67 static int niter = 10000 ;
68
69 /*-- global variables for ACF simulation [30 Nov 2015] --*/
70
71 static int do_acf = 0 ;
72 static float acf_a = 0.0f ;
73 static float acf_b = 0.0f ;
74 static float acf_c = 0.0f ;
75 static float acf_parm[3] = {0.0f,0.0f,0.0f} ;
76 static int acf_nxx=0 , acf_nyy=0 , acf_nzz=0 ;
77 static MRI_IMAGE **acf_aim=NULL , **acf_bim=NULL ;
78 static complex **acf_tar=NULL ; static int acf_ntar=0 ;
79 static MRI_IMAGE *acf_wim=NULL ;
80
81 static int do_classic = 0 ; /* does nothing (yet) */
82
83 static int ex_pad = 0 ; /* 12 May 2015 -- allow for padding the volume */
84 static int ey_pad = 0 ;
85 static int ez_pad = 0 ;
86 static int do_pad = 0 ;
87 static int nx_pad ;
88 static int ny_pad ;
89 static int nz_pad ;
90 static int nxy_pad ;
91 static int nxyz_pad ;
92 static int allow_padding = 1 ; /* on by default */
93
94 static float tdof = 0.0f ; /* 26 May 2015 - secret stuff */
95
96 static float sigmax , sigmay , sigmaz ;
97 static int do_blur = 0 ;
98
99 static int do_ssave = 0 ; /* 24 Apr 2014 */
100 static char *ssave_prefix = NULL ;
101 static int ssave_index = 0 ;
102 static THD_3dim_dataset *ssave_dset = NULL ;
103
104 #define SSAVE_BLURRED 1
105 #define SSAVE_MASKED 2
106 #define SSAVE_ABSOLUTE 3 /* no longer used */
107
108 static int nodec = 0 ;
109 static int do_niml = 0 ;
110 static int do_1D = 1 ;
111 static int do_ball = 0 ;
112
113 static unsigned int gseed = 123456789 ;
114
115 #define PMAX 0.5
116
117 static double pthr_init[9] = { 0.05 , 0.02 , 0.01, 0.005, 0.002, 0.001, 0.0005, 0.0002, 0.0001 } ;
118 static double athr_init[4] = { 0.10 , 0.05 , 0.02 , 0.01 } ;
119
120 static int npthr_lots = 29 ;
121 static double pthr_lots[29] = { 0.10, 0.09, 0.08, 0.07, 0.06,
122 0.05, 0.04, 0.03, 0.02, 0.015, 0.01,
123 0.007, 0.005, 0.003, 0.002, 0.0015, 0.001,
124 0.0007, 0.0005, 0.0003, 0.0002, 0.00015, 0.0001,
125 0.00007, 0.00005, 0.00003, 0.00002, 0.000015, 0.00001 } ;
126
127 static int npthr_mega = 38 ;
128 static double pthr_mega[38] = { 0.100, 0.090, 0.080, 0.070, 0.060,
129 0.050, 0.045, 0.040, 0.035, 0.030,
130 0.025, 0.020, 0.015, 0.010, 0.009,
131 0.008, 0.007, 0.006, 0.005, 0.004,
132 0.003, 0.002, 0.001, 0.0009, 0.0008,
133 0.0007, 0.0006, 0.0005, 0.0004, 0.0003,
134 0.0002, 0.0001, 0.00007, 0.00005, 0.00003,
135 0.00002, 0.000015,0.00001 } ;
136
137 static int nathr_lots = 10 ;
138 static double athr_lots[10] = { 0.10, 0.09, .08, .07, .06, .05, .04, .03, .02, .01 } ;
139
140 static int nathr_mega = 22 ;
141 static double athr_mega[22] = { 0.100,0.095,0.090,0.085,0.080,0.075,0.070,0.065,0.060,0.055,
142 0.050,0.045,0.040,0.035,0.030,0.025,0.020,0.015,0.010,0.005,0.004,0.003 } ;
143
144 static int force_mega = 0 ; /* 15 Dec 2021 */
145 static int did_mega = 0 ;
146
147 #undef SET_MEGA /* 15 Dec 2021 */
148 #define SET_MEGA \
149 do{ npthr = npthr_mega ; \
150 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ; \
151 memcpy( pthr , pthr_mega , sizeof(double)*npthr ) ; \
152 nathr = nathr_mega ; \
153 athr = (double *)realloc(athr,sizeof(double)*nathr) ; \
154 memcpy( athr , athr_mega , sizeof(double)*nathr ) ; \
155 athr_sum_bot = 13 ; athr_sum_top = 31 ; \
156 if( niter < 30000 ) niter = 30000 ; \
157 did_mega = 1 ; \
158 } while(0)
159
160 static int npthr = 9 ;
161 static double *pthr = NULL ;
162
163 static float *zthr_1sid = NULL ;
164 static float *zthr_2sid = NULL ;
165
166 static int nathr = 4 ;
167 static double *athr = NULL ;
168
169 /* the output:
170 2D table of cluster size threshold as function
171 of per-voxel threshold and volumetric alpha
172 (will be interpolated from simulation results) */
173
174 static float **clust_thresh_1sid_NN1 = NULL ;
175 static float **clust_thresh_2sid_NN1 = NULL ;
176 static float **clust_thresh_bsid_NN1 = NULL ;
177
178 static float **clust_thresh_1sid_NN2 = NULL ;
179 static float **clust_thresh_2sid_NN2 = NULL ;
180 static float **clust_thresh_bsid_NN2 = NULL ;
181
182 static float **clust_thresh_1sid_NN3 = NULL ;
183 static float **clust_thresh_2sid_NN3 = NULL ;
184 static float **clust_thresh_bsid_NN3 = NULL ;
185
186 static int do_athr_sum = 0 ; /* 18 Dec 2015 */
187 static int athr_sum_bot=-1 , athr_sum_top=-1 ;
188
189 /* max_table_1sid[nnn][ipthr][cc] is the count
190 of how often the maximum cluster had exact size cc voxels
191 (out of niter trials), at the pthr[ipthr] voxel-wise threshold p-value,
192 for the nnn-th NN level (1-3), for 1-sided voxel-wise thresholding.
193 Mutatis mutandis for max_table_2sid and max_table_bsid (bi-sided threshold) */
194
195 static int **max_table_1sid[4] , **max_table_2sid[4] , **max_table_bsid[4] ;
196
197 /* cmx_table_1sid[nnn][ipthr] = max cluster size found
198 over all simulations for 1-sided clustering with NN level nnn
199 at voxel-wise threshold pthr[ipthr]. That is, it is the largest
200 value of cc such that max_table_1sid[nnn][ipthr][cc] > 0. [23 Dec 2015] */
201
202 static int *cmx_table_1sid[4] , *cmx_table_2sid[4] , *cmx_table_bsid[4] ;
203
204 /* csiz_1D_NN1[ipthr][iter] = max cluster size found
205 at the ipthr-th threshold in the iter-th iteration. [23 Dec 2015] */
206
207 static int **csiz_1sid_NN1 , **csiz_1sid_NN2 , **csiz_1sid_NN3 ;
208 static int **csiz_2sid_NN1 , **csiz_2sid_NN2 , **csiz_2sid_NN3 ;
209 static int **csiz_bsid_NN1 , **csiz_bsid_NN2 , **csiz_bsid_NN3 ;
210
211 #undef USE_SHAVE
212 #ifdef USE_SHAVE
213 #define SHAVE_MALLOC 1
214 #define SHAVE_MMAP 2
215 static int do_shave =0 ;
216 size_t shave_siz=0 ;
217 static int64_t shave_tot=0 ;
218 static short *shave =NULL ;
219 #endif
220
221 static int verb = 1 ;
222 static int nthr = 1 ;
223
224 #undef DECLARE_ithr /* 30 Nov 2015 */
225 #ifdef USE_OMP
226 # define DECLARE_ithr const int ithr=omp_get_thread_num()
227 #else
228 # define DECLARE_ithr const int ithr=0
229 #endif
230
231 static int minmask = 128 ; /* 29 Mar 2011 */
232
233 static char *prefix = NULL ;
234
235 static THD_3dim_dataset **inset = NULL ; /* 02 Feb 2016 */
236 static int num_inset = 0 ;
237 static int *nval_inset = NULL ;
238
239 static Xdataset *xinset = NULL ; /* 29 Aug 2016 */
240
241 #undef PSMALL
242 #define PSMALL 1.e-15
243
244 static char *cmd_fname = "3dClustSim.cmd" ;
245
246 /*----------------------------------------------------------------------------*/
247 /*! Threshold for upper tail probability of N(0,1) */
248
zthresh(double pval)249 double zthresh( double pval )
250 {
251 if( pval <= 0.0 ) pval = PSMALL ;
252 else if( pval >= 1.0 ) pval = 1.0 - PSMALL ;
253 return qginv(pval) ;
254 }
255
256 /*---------------------------------------------------------------------------*/
257
258 static int vsnn=0 ;
259
vstep_reset(void)260 static void vstep_reset(void){ vsnn=0; }
261
vstep_print(void)262 static void vstep_print(void)
263 {
264 static char xx[10] = "0123456789" ;
265 fprintf(stderr , "%c" , xx[vsnn%10] ) ;
266 if( vsnn%10 == 9) fprintf(stderr,".") ;
267 vsnn++ ;
268 }
269
270 /*---------------------------------------------------------------------------*/
271
display_help_menu()272 void display_help_menu()
273 {
274 printf(
275 "Usage: 3dClustSim [options]\n"
276 "\n"
277 "Program to estimate the probability of false positive (noise-only) clusters.\n"
278 "An adaptation of Doug Ward's AlphaSim, streamlined for various purposes.\n"
279 "\n"
280 "-----------------------------------------------------------------------------\n"
281 "This program has several different modes of operation, each one involving\n"
282 "simulating noise-only random volumes, thresholding and clustering them,\n"
283 "and counting statistics of how often data 'survives' these processes at\n"
284 "various threshold combinations (per-voxel and cluster-size).\n"
285 "\n"
286 "OLDEST method = simulate noise volume assuming the spatial auto-correlation\n"
287 " function (ACF) is given by a Gaussian-shaped function, where\n"
288 " this shape is specified using the FWHM parameter. The FWHM\n"
289 " parameter can be estimated by program 3dFWHMx.\n"
290 " ** THIS METHOD IS NO LONGER RECOMMENDED **\n"
291 "\n"
292 "NEWER method = simulate noise volume assuming the ACF is given by a mixed-model\n"
293 " of the form a*exp(-r*r/(2*b*b))+(1-a)*exp(-r/c), where a,b,c\n"
294 " are 3 parameters giving the shape, and can also be estimated\n"
295 " by program 3dFWHMx.\n"
296 " ** THIS METHOD IS ACCEPTABLE **\n"
297 "\n"
298 "NEWEST method = program 3dttest++ simulates the noise volumes by randomizing\n"
299 " and permuting input datasets, and sending those volumes into\n"
300 " 3dClustSim directly. There is no built-in math model for the\n"
301 " spatial ACF.\n"
302 " ** THIS METHOD IS MOST ACCURATE AT CONTROLLING FALSE POSITIVE RATE **\n"
303 " ** You invoke this method with the '-Clustsim' option in 3dttest++ **\n"
304 "\n"
305 "3dClustSim computes a cluster-size threshold for a given voxel-wise p-value\n"
306 "threshold, such that the probability of anything surviving the dual thresholds\n"
307 "is at some given level (specified by the '-athr' option).\n"
308 "\n"
309 "Note that this cluster-size threshold is the same for all brain regions.\n"
310 "There is an implicit assumption that the noise spatial statistics are\n"
311 "the same everywhere.\n"
312 "\n"
313 "Program 3dXClustSim introduces the idea of spatially variable cluster-size\n"
314 "thresholds, which may be more useful in some cases. 3dXClustSim's method is\n"
315 "invoked by using the '-ETAC' option in 3dttest++.\n"
316 "-----------------------------------------------------------------------------\n"
317 "\n"
318 "**** NOTICE ****\n"
319 "You should use the -acf method, NOT the -fwhm method, when determining\n"
320 "cluster-size thresholds for FMRI data. The -acf method will give more\n"
321 "accurate false positive rate (FPR) control.\n"
322 "****************\n"
323 "\n"
324 "In particular, this program lets you run with multiple p-value thresholds\n"
325 "(the '-pthr' option) and only outputs the cluster size threshold at chosen\n"
326 "values of the alpha significance level (the '-athr' option).\n"
327 "\n"
328 "In addition, the program allows the output to be formatted for inclusion\n"
329 "into an AFNI dataset's header, whence it can be used in the AFNI Clusterize\n"
330 "interface to show approximate alpha values for the displayed clusters, where\n"
331 "the per-voxel p-value is taken from the interactive threshold slider in the\n"
332 "AFNI 'Define Overlay' control panel, and then the per-cluster alpha value\n"
333 "is interpolated in this table from 3dClustSim. As you change the threshold\n"
334 "slider, the per-voxel p-value (shown below the slider) changes, and then\n"
335 "the interpolated alpha values are updated.\n"
336 "\n"
337 "************* IMPORTANT NOTE [Dec 2015] ***************************************\n"
338 "A completely new method for estimating and using noise smoothness values is\n"
339 "now available in 3dFWHMx and 3dClustSim. This method is implemented in the\n"
340 "'-acf' options to both programs. 'ACF' stands for (spatial) AutoCorrelation\n"
341 "Function, and it is estimated by calculating moments of differences out to\n"
342 "a larger radius than before.\n"
343 "\n"
344 "Notably, real FMRI data does not actually have a Gaussian-shaped ACF, so the\n"
345 "estimated ACF is then fit (in 3dFWHMx) to a mixed model (Gaussian plus\n"
346 "mono-exponential) of the form\n"
347 " ACF(r) = a * exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)\n"
348 "where 'r' is the radius, and 'a', 'b', 'c' are the fitted parameters.\n"
349 "The apparent FWHM from this model is usually somewhat larger in real data\n"
350 "than the FWHM estimated from just the nearest-neighbor differences used\n"
351 "in the 'classic' analysis.\n"
352 "\n"
353 "The longer tails provided by the mono-exponential are also significant.\n"
354 "3dClustSim has also been modified to use the ACF model given above to generate\n"
355 "noise random fields.\n"
356 "\n"
357 "**----------------------------------------------------------------------------**\n"
358 "** The take-away (TL;DR or summary) message is that the 'classic' 3dFWHMx and **\n"
359 "** 3dClustSim analysis, using a pure Gaussian ACF, is not very correct for **\n"
360 "** FMRI data -- I cannot speak for PET or MEG data. **\n"
361 #if 0
362 "** **\n"
363 "** You should start using the '-acf' options in your own scripts. AFNI **\n"
364 "** scripts from afni_proc.py are moving away from the 'classic' method. At **\n"
365 "** some point in the future, '-acf' will become the default in both programs, **\n"
366 "** and you will have to actively specify '-classic' to use the older method. **\n"
367 #endif
368 "**----------------------------------------------------------------------------**\n"
369 "\n"
370 "** ---------------------------------------------------------------------------**\n"
371 "** IMPORTANT CHANGES -- February 2015 ******************************************\n"
372 "** ---------------------------------------------------------------------------**\n"
373 "** In the past, 3dClustSim did '1-sided' testing; that is, the random dataset\n"
374 "** of Gaussian noise-only values is generated, and then it is thresholded on\n"
375 "** the positive side so that the N(0,1) upper tail probability is pthr.\n"
376 "**\n"
377 "** NOW, 3dClustSim does 3 different types of thresholding:\n"
378 "** 1-sided: as above\n"
379 "** 2-sided: where positive and negative values above the threshold\n"
380 "** are included, and then clustered together\n"
381 "** (in this case, the threshold on the Gaussian values is)\n"
382 "** (fixed so that the 1-sided tail probability is pthr/2.)\n"
383 "** bi-sided: where positive values and negative values above the\n"
384 "** threshold are clustered SEPARATELY (with the 2-sided threshold)\n"
385 "** For high levels of smoothness, the results from bi-sided and 2-sided are\n"
386 "** very similar -- since for smooth data, it is unlikely that large clusters of\n"
387 "** positive and negative values will be next to each other. With high smoothness,\n"
388 "** it is also true that the 2-sided results for 2*pthr will be similar to the\n"
389 "** 1-sided results for pthr, for the same reason. Since 3dClustSim is meant to be\n"
390 "** useful when the noise is NOT very smooth, we provide tables for all 3 cases.\n"
391 "**\n"
392 "** In particular, note that when the AFNI GUI threshold is set to a t-statistic,\n"
393 "** 2-sided testing is what is usually appropriate -- in that case, the cluster\n"
394 "** size thresholds tend to be smaller than the 1-sided case, which means that\n"
395 "** more clusters tend to be significant than in the past.\n"
396 "**\n"
397 "** In addition, the 3 different NN approaches (NN=1, NN=2, NN=3) are ALL\n"
398 "** always computed now. That is, 9 different tables are produced, each\n"
399 "** of which has its proper place when combined with the AFNI Clusterize GUI.\n"
400 "** The 3 different NN methods are:\n"
401 "** 1 = Use first-nearest neighbor clustering\n"
402 "** * above threshold voxels cluster together if faces touch\n"
403 "** 2 = Use second-nearest neighbor clustering\n"
404 "** * voxels cluster together if faces OR edges touch\n"
405 "** 3 = Use third-nearest neighbor clustering\n"
406 "** * voxels cluster together if faces OR edges OR corners touch\n"
407 "** The clustering method only makes a difference at higher (less significant)\n"
408 "** values of pthr. At small values of pthr (more significant), all three\n"
409 "** clustering methods will give very similar results.\n"
410 "**\n"
411 "**** PLEASE NOTE that the NIML outputs from this new version are not named the\n"
412 "**** same as those from the older version. Thus, any script that takes the NIML\n"
413 "**** format tables and inserts them into an AFNI dataset header must be modified\n"
414 "**** to match the new names. The 3drefit command fragment output at the end of\n"
415 "**** this program (and echoed into file '3dClustSim.cmd') shows the new form\n"
416 "**** of the names involved.\n"
417 "**** -------------------------------------------------------------------------**\n"
418 "**** SMOOTHING CHANGE -- May 2015 **********************************************\n"
419 "** ---------------------------------------------------------------------------**\n"
420 "** It was pointed out to me (by Anders Eklund and Tom Nichols) that smoothing\n"
421 "** the simulated data over a finite volume introduces 2 artifacts, which might\n"
422 "** be called 'edge effects'. To minimize these problems, this program now makes\n"
423 "** extra-large (padded) simulated volumes before blurring, and then trims those\n"
424 "** back down to the desired size, before continuing with the thresholding and\n"
425 "** cluster-counting steps. To run 3dClustSim without this padding added, use\n"
426 "** the new '-nopad' option.\n"
427 #if 0
428 "** Also see the manuscript\n"
429 "** A Eklund, T Nichols, M Andersson, and H Knutsson.\n"
430 "** Empirically Investigating the Statistical Validity of SPM, FSL and AFNI\n"
431 "** for Single Subject FMRI Analysis.\n"
432 #endif
433 "**** -------------------------------------------------------------------------**\n"
434 "\n"
435 "-------\n"
436 "OPTIONS [at least 1 option is required, or you'll get this help message!]\n"
437 "-------\n"
438 " ******* Specify the volume over which the simulation will occur *******\n"
439 "\n"
440 " -----** (a) Directly give the spatial domain that will be used **-----\n"
441 "\n"
442 "-nxyz n1 n2 n3 = Size of 3D grid to use for simulation\n"
443 " [default values = 64 64 32]\n"
444 "-dxyz d1 d2 d3 = give all 3 voxel sizes at once\n"
445 " [default values = 3.5 3.5 3.5]\n"
446 "-BALL = inside the 3D grid, mask off points outside a ball\n"
447 " at the center of the grid and touching the edges;\n"
448 " this will keep about 1/2 the points in the 3D grid.\n"
449 " [default = use all voxels in the 3D grid]\n"
450 "\n"
451 " -----** OR: (b) Specify the spatial domain using a dataset mask **-----\n"
452 "\n"
453 "-mask mset = Use the 0 sub-brick of dataset 'mset' as a mask\n"
454 " to indicate which voxels to analyze (a sub-brick\n"
455 " selector '[]' is allowed) \n"
456 "\n"
457 "-OKsmallmask = Allow small masks. Normally, a mask volume must have\n"
458 " 128 or more nonzero voxels. However, IF you know what\n"
459 " you are doing, and IF you are willing to live life on\n"
460 " the edge of statistical catastrophe, then you can use\n"
461 " this option to allow smaller masks -- in a sense, this\n"
462 " is the 'consent form' for such strange shenanigans.\n"
463 " * If you use this option, it must come BEFORE '-mask'.\n"
464 " * Also read the 'CAUTION and CAVEAT' section, far below.\n"
465 " -->>** This option is really only recommended for users who\n"
466 " understand what they are doing. Misuse of this option\n"
467 " could easily be construed as 'p-hacking'; for example,\n"
468 " finding results, but your favorite cluster is too small\n"
469 " to survive thresholding, so you post-hoc put a small mask\n"
470 " down in that region. DON'T DO THIS!\n"
471 "\n"
472 " ** '-mask' means that '-nxyz' & '-dxyz' & '-BALL' will be ignored. **\n"
473 "\n"
474 " -----** OR: (c) Specify the spatial domain by directly giving simulated volumes **-----\n"
475 "\n"
476 "-inset iset [iset ...] = Read the 'iset' dataset(s) and use THESE volumes\n"
477 " as the simulations to threshold and clusterize,\n"
478 " [Feb 2016] rather than create the simulations internally.\n"
479 " * For example, these datasets could come from\n"
480 " 3dttest++ -toz -randomsign 1000 -setA ...\n"
481 " * This can be combined with '-mask'.\n"
482 " * Using '-inset' means that '-fwhm', '-acf', '-nopad',\n"
483 " '-niter', and '-ssave' are ignored as meaningless.\n"
484 "\n"
485 " ---** the remaining options control how the simulation is done **---\n"
486 "\n"
487 "-fwhm s = Gaussian filter width (all 3 dimensions) in mm (non-negative)\n"
488 " [default = 0.0 = no smoothing]\n"
489 " * If you wish to set different smoothing amounts for each\n"
490 " axis, you can instead use option\n"
491 " -fwhmxyz sx sy sz\n"
492 " to specify the three values separately.\n"
493 " **** This option is no longer recommended, since FMRI data ****\n"
494 " **** does not have a Gaussian-shaped spatial autocorrelation. ****\n"
495 " **** Consider using '-acf' or '3dttest++ -Clustsim' instead. ****\n"
496 "\n"
497 "-acf a b c = Alternative to Gaussian filtering: use the spherical\n"
498 " autocorrelation function parameters output by 3dFWHMx\n"
499 " to do non-Gaussian (long-tailed) filtering.\n"
500 " * Using '-acf' will make '-fwhm' pointless!\n"
501 " * The 'a' parameter must be between 0 and 1.\n"
502 " * The 'b' and 'c' parameters (scale radii) must be positive.\n"
503 " * The spatial autocorrelation function is given by\n"
504 " ACF(r) = a * exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)\n"
505 " >>---------->>*** Combined with 3dFWHMx, the '-acf' method is now a\n"
506 " recommended way to generate clustering statistics in AFNI!\n"
507 " *** Alternative methods we also recommend:\n"
508 " 3dttest++ with the -Clustsim and/or -ETAC options.\n"
509 "\n"
510 "-nopad = The program now [12 May 2015] adds 'padding' slices along\n"
511 " each face to allow for edge effects of the smoothing process.\n"
512 " If you want to turn this feature off, use the '-nopad' option.\n"
513 " * For example, if you want to compare the 'old' (un-padded)\n"
514 " results with the 'new' (padded) results.\n"
515 " * '-nopad' has no effect when '-acf' is used, since that option\n"
516 " automatically pads the volume when creating it (via FFTs) and\n"
517 " then truncates it back to the desired size for clustering.\n"
518 "\n"
519 "-pthr p1 .. pn = list of uncorrected (per voxel) p-values at which to\n"
520 " threshold the simulated images prior to clustering.\n"
521 " [default = 0.05 0.02 0.01 0.005 0.002 0.001 0.0005 0.0002 0.0001]\n"
522 "\n"
523 "-athr a1 .. an = list of corrected (whole volume) alpha-values at which\n"
524 " the simulation will print out the cluster size\n"
525 " thresholds. For each 'p' and 'a', the smallest cluster\n"
526 " size C(p,a) for which the probability of the 'p'-thresholded\n"
527 " image having a noise-only cluster of size C is less than 'a'\n"
528 " is the output (cf. the sample output, below)\n"
529 " [default = 0.10 0.05 0.02 0.01]\n"
530 "\n"
531 " ** It is possible to use only ONE value in each of '-pthr' and **\n"
532 " ** '-athr', and then you will get exactly one line of output **\n"
533 " ** for each sided-ness and NN case. For example: **\n"
534 " ** -pthr 0.001 -athr 0.05 **\n"
535 "\n"
536 " ** Both lists '-pthr' and '-athr' (of values between 0 and 0.2) **\n"
537 " ** should be given in DESCENDING order. They will be sorted to be **\n"
538 " ** that way in any case, and such is how the output will be given. **\n"
539 "\n"
540 " ** The list of values following '-pthr' or '-athr' can be replaced **\n"
541 " ** with the single word 'LOTS', which will tell the program to use **\n"
542 " ** a longer list of values for these probabilities [try it & see!] **\n"
543 " ** (i.e., '-pthr LOTS' and/or '-athr LOTS' are legal options) **\n"
544 "\n"
545 "-LOTS = the same as using '-pthr LOTS -athr LOTS'\n"
546 "-MEGA = adds even MORE values to the '-pthr' and '-athr' grids.\n"
547 " * NOTE: you can also invoke '-MEGA' by setting environment\n"
548 " variable AFNI_CLUSTSIM_MEGA to YES.\n"
549 " * Doing this will over-ride any use of other options to set\n"
550 " the '-pthr' and '-athr' lists!\n"
551 "\n"
552 "-iter n = number of Monte Carlo simulations [default = 10000]\n"
553 "\n"
554 "-nodec = normally, the program prints the cluster size threshold to\n"
555 " 1 decimal place (e.g., 27.2). Of course, clusters only come\n"
556 " with an integer number of voxels -- this fractional value\n"
557 " is interpolated to give the desired alpha level. If you\n"
558 " want no decimal places (so that 27.2 becomes 28), use '-nodec'.\n"
559 "\n"
560 "-seed S = random number seed [default seed = 123456789]\n"
561 " * if seed=0, then program will quasi-randomize it\n"
562 "\n"
563 "-niml = Output the table in an XML/NIML format, rather than a .1D format.\n"
564 " * This option is for use with other software programs;\n"
565 " see the NOTES section below for details.\n"
566 " * '-niml' also implicitly means '-LOTS'.\n"
567 "\n"
568 "-both = Output the table in XML/NIML format AND in .1D format.\n"
569 " * You probably want to use '-prefix' with this option!\n"
570 " Otherwise, everything is mixed together on stdout.\n"
571 " * '-both' implies 'niml' which implies '-LOTS' (unless '-MEGA').\n"
572 " So '-pthr' (if desired) should follow '-both'/'-niml'\n"
573 "\n"
574 "-prefix ppp = Write output for NN method #k to file 'ppp.NNk_Xsided.1D',\n"
575 " for k=1, 2, 3, and for X=1sided, 2sided, bisided.\n"
576 " * If '-prefix is not used, all results go to standard output.\n"
577 " You will probably find this confusing.\n"
578 " * If '-niml' is used, the filename is 'ppp.NNk_Xsided.niml'.\n"
579 " To be clear, the 9 files that will be named\n"
580 " ppp.NN1_1sided.niml ppp.NN1_2sided.niml ppp.NN1_bisided.niml\n"
581 " ppp.NN2_1sided.niml ppp.NN2_2sided.niml ppp.NN2_bisided.niml\n"
582 " ppp.NN3_1sided.niml ppp.NN3_2sided.niml ppp.NN3_bisided.niml\n"
583 " * If '-niml' AND '-mask' are both used, then a compressed ASCII\n"
584 " encoding of the mask volume is stored into file 'ppp.mask'.\n"
585 " This string can be stored into a dataset header as an attribute\n"
586 " with name AFNI_CLUSTSIM_MASK, and will be used in the AFNI\n"
587 " Clusterize GUI, if present, to mask out above-threshold voxels\n"
588 " before the clusterizing is done (which is how the mask is used\n"
589 " here in 3dClustSim).\n"
590 " * If the ASCII mask string is NOT stored into the statistics dataset\n"
591 " header, then the Clusterize GUI will try to find the original\n"
592 " mask dataset and use that instead. If that fails, then masking\n"
593 " won't be done in the Clusterize process.\n"
594 "\n"
595 " -cmd ccc = Write command for putting results into a file's header to a file\n"
596 " named 'ccc' instead of '3dClustSim.cmd'. This option is mostly\n"
597 " to help with scripting, as in\n"
598 " 3dClustSim -cmd XXX.cmd -prefix XXX.nii ...\n"
599 " `cat XXX.cmd` XXX.nii\n"
600 "\n"
601 " -quiet = Don't print out the progress reports, etc.\n"
602 " * Put this option first to silence most informational messages.\n"
603 "\n"
604 " -ssave:TYPE ssprefix = Save the un-thresholded generated random volumes into\n"
605 " datasets ('-iter' of them). Here, 'TYPE' is one of these:\n"
606 " * blurred == save the blurred 3D volume before masking\n"
607 " * masked == save the blurred volume after masking\n"
608 " The output datasets will actually get prefixes generated\n"
609 " with the string 'ssprefix' being appended by a 6 digit\n"
610 " integer (the iteration index), starting at 000000.\n"
611 " (You can use SOMETHING.nii as a prefix; it will work OK.)\n"
612 " N.B.: This option will slow the program down a lot,\n"
613 " and was intended to help just one specific user.\n"
614 "\n"
615 "------\n"
616 "NOTES:\n"
617 "------\n"
618 "* This program is like running AlphaSim once for each '-pthr' value and then\n"
619 " extracting the relevant information from its 'Alpha' output column.\n"
620 " ++ One reason for 3dClustSim to be used in place of AlphaSim is that it will\n"
621 " be much faster than running AlphaSim multiple times.\n"
622 " ++ Another reason is that the resulting table can be stored in an AFNI\n"
623 " dataset's header, and used in the AFNI Clusterize GUI to see estimated\n"
624 " cluster significance (alpha) levels.\n"
625 "\n"
626 "* To be clear, the C(p,alpha) thresholds that are calculated are for\n"
627 " alpha = probability of a noise-only smooth random field, after masking\n"
628 " and then thresholding at the given per-voxel p value, producing a cluster\n"
629 " of voxels at least this big.\n"
630 " ++ So if your cluster is larger than the C(p,0.01) threshold in size (say),\n"
631 " then it is very unlikely that noise BY ITSELF produced this result.\n"
632 " ++ This statement does not mean that ALL the voxels in the cluster are\n"
633 " 'truly' active -- it means that at least SOME of them are (very probably)\n"
634 " active. The statement of low probability (0.01 in this example) of a\n"
635 " false positive result applies to the cluster as a whole, not to each\n"
636 " voxel within the cluster.\n"
637 "\n"
638 "* To add the cluster simulation C(p,alpha) table to the header of an AFNI\n"
639 " dataset, something like the following can be done [tcsh syntax]:\n"
640 " set fx = ( `3dFWHMx -detrend time_series_dataset+orig` )\n"
641 " 3dClustSim -mask mask+orig -acf $fx[5] $fx[6] $fx[7] -niml -prefix CStemp\n"
642 " 3drefit -atrstring AFNI_CLUSTSIM_NN1_1sided file:CStemp.NN1_1sided.niml \\\n"
643 " -atrstring AFNI_CLUSTSIM_MASK file:CStemp.mask \\\n"
644 " statistics_dataset+orig\n"
645 " rm -f CStemp.*\n"
646 " AFNI's Clusterize GUI makes use of these attributes, if stored in a\n"
647 " statistics dataset (e.g., something from 3dDeconvolve, 3dREMLfit, etc.).\n"
648 "\n"
649 " ** Nota Bene: afni_proc.py will automatically run 3dClustSim, and **\n"
650 " *** put the results into the statistical results dataset for you. ***\n"
651 " **** Another reason to use afni_proc.py for single-subject analyses! ****\n"
652 "\n"
653 "* 3dClustSim will print (to stderr) a 3drefit command fragment, similar\n"
654 " to the one above, that you can use to add cluster tables to any\n"
655 " relevant statistical datasets you have lolling about.\n"
656 "\n"
657 "* The C(p,alpha) table will be used in Clusterize to provide the cluster\n"
658 " level alpha value when the AFNI GUI is set so that the Overlay threshold\n"
659 " sub-brick is a statistical parameter (e.g., a t- or F-statistic), from which\n"
660 " a per-voxel p-value can be calculated, so that Clusterize can interpolate\n"
661 " in the C(p,alpha) table.\n"
662 " ++ To be clear, the per-voxel p-value is taken from the AFNI GUI threshold\n"
663 " slider (the p-value is shown beneath the slider), and then the C(p,alpha)\n"
664 " table is inverse-interpolated to find the per-cluster alpha value for\n"
665 " each different cluster size.\n"
666 " ++ As you move the AFNI threshold slider, the per-voxel (uncorrected for\n"
667 " multiple comparisons) p-value changes, the cluster sizes change (as fewer\n"
668 " or more voxels are included), and so the reported per-cluster alpha\n"
669 " values change for both reasons -- different p and different cluster size.\n"
670 " ++ The alpha values reported are 'per-cluster', and are not themselves\n"
671 " corrected for multiple comparisons ACROSS clusters. These alpha values\n"
672 " are corrected for multiple comparisons WITHIN a cluster.\n"
673 "\n"
674 "* AFNI will use the NN1, NN2, NN3 tables as needed in its Clusterize\n"
675 " interface if they are all stored in the statistics dataset header,\n"
676 " depending on the NN level chosen in the Clusterize controller.\n"
677 "\n"
678 "* The blur estimates (provided to 3dClustSim via -acf) comes from using\n"
679 " program 3dFWHMx.\n"
680 "\n"
681 "-------------------\n"
682 "CAUTION and CAVEAT: [January 2011]\n"
683 "-------------------\n"
684 "* If you use a small ROI mask and also have a large blur, then it might happen\n"
685 " that it is impossible to find a cluster size threshold C that works for a\n"
686 " given (p,alpha) combination.\n"
687 "\n"
688 "* Generally speaking, C(p,alpha) gets smaller as p gets smaller and C(p,alpha)\n"
689 " gets smaller as alpha gets larger. As a result, in a small mask with small p\n"
690 " and large alpha, C(p,alpha) might shrink below 1. But clusters of size C\n"
691 " less than 1 don't make any sense!\n"
692 "\n"
693 "* For example, suppose that for p=0.0005 that only 6%% of the simulations\n"
694 " have ANY above-threshold voxels inside the ROI mask. In that case,\n"
695 " C(p=0.0005,alpha=0.06) = 1. There is no smaller value of C where 10%%\n"
696 " of the simulations have a cluster of size C or larger. Thus, it is\n"
697 " impossible to find the cluster size threshold for the combination of\n"
698 " p=0.0005 and alpha=0.10 in this case.\n"
699 "\n"
700 "* 3dClustSim will report a cluster size threshold of C=1 for such cases.\n"
701 " It will also print (to stderr) a warning message for all the (p,alpha)\n"
702 " combinations that had this problem.\n"
703 #if 0
704 "\n"
705 "* This issue arises because 3dClustSim reports C for a given alpha.\n"
706 " In contrast, AlphaSim reports alpha for each given C, and leaves\n"
707 " you to interpret the resulting table; it doesn't try to find C(p,alpha)\n"
708 " for a given alpha, but finds alpha for various values of C.\n"
709 "\n"
710 "* If you wish to see this effect in action, the following commands\n"
711 " can be used as a starting point:\n"
712 "\n"
713 "3dClustSim -nxyz 8 8 8 -dxyz 2 2 2 -fwhm 8 -niter 10000\n"
714 "AlphaSim -nxyz 8 8 8 -dxyz 2 2 2 -fwhm 8 -niter 10000 -quiet -fast -pthr 0.0005\n"
715 "\n"
716 " From the 3dClustSim command above, you should get a warning message\n"
717 " similar to this, right after the table (only the first 2 lines are shown):\n"
718 "\n"
719 "*+ WARNING: Simulation not effective for these cases:\n"
720 " NN=1 pthr= 0.001000 alpha= 0.100 [max simulated alpha= 0.087]\n"
721 #endif
722 "\n"
723 "-----------------------------\n"
724 "---- RW Cox -- July 2010 ----\n"
725 ) ;
726
727 printf(
728 "\n"
729 "-------------\n"
730 "SAMPLE OUTPUT from the command '3dClustSim -fwhm 7' [only the NN=1 1-sided results]\n"
731 "-------------\n"
732 "# 3dClustSim -fwhm 7\n"
733 "# 1-sided thresholding\n"
734 "# Grid: 64x64x32 3.50x3.50x3.50 mm^3 (131072 voxels)\n"
735 "#\n"
736 "# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels\n"
737 "# -NN 1 | alpha = Prob(Cluster >= given size)\n"
738 "# pthr | 0.100 0.050 0.020 0.010\n"
739 "# ------ | ------ ------ ------ ------\n"
740 " 0.050000 162.5 182.2 207.8 225.7\n"
741 " 0.020000 64.3 71.0 80.5 88.5\n"
742 " 0.010000 40.3 44.7 50.7 55.1\n"
743 " 0.005000 28.0 31.2 34.9 38.1\n"
744 " 0.002000 19.0 21.2 24.2 26.1\n"
745 " 0.001000 14.6 16.3 18.9 20.5\n"
746 " 0.000500 11.5 13.0 15.1 16.7\n"
747 " 0.000200 8.7 10.0 11.6 12.8\n"
748 " 0.000100 7.1 8.3 9.7 10.9\n"
749 "\n"
750 "e.g., for this sample volume, if the per-voxel p-value threshold is set\n"
751 "at 0.005, then to keep the probability of getting a single noise-only\n"
752 "cluster at 0.05 or less, the cluster size threshold should be 32 voxels\n"
753 "(the next integer above 31.2).\n"
754 "\n"
755 "If you ran the same simulation with the '-nodec' option, then the last\n"
756 "line above would be\n"
757 " 0.000100 8 9 10 11\n"
758 "If you set the per voxel p-value to 0.0001 (1e-4), and want the chance\n"
759 "of a noise-only false-positive cluster to be 5%% or less, then the cluster\n"
760 "size threshold would be 9 -- that is, you would keep all NN clusters with\n"
761 "9 or more voxels.\n"
762 "\n"
763 "The header lines start with the '#' (commenting) character so that the result\n"
764 "is a correctly formatted AFNI .1D file -- it can be used in 1dplot, etc.\n"
765 ) ;
766
767 PRINT_AFNI_OMP_USAGE("3dClustSim",NULL) ;
768 PRINT_COMPILE_DATE ;
769 exit(0);
770 }
771
772 /*---------------------------------------------------------------------------*/
773 /* Routine to initialize the input options (values are in global variables). */
774
get_options(int argc,char ** argv)775 void get_options( int argc , char **argv )
776 {
777 char * ep;
778 int nopt=1 , ii , have_pthr=0;
779
780 ENTRY("get_options") ;
781
782 pthr = (double *)malloc(sizeof(double)*npthr) ;
783 memcpy( pthr , pthr_init , sizeof(double)*npthr ) ;
784
785 athr = (double *)malloc(sizeof(double)*nathr) ;
786 memcpy( athr , athr_init , sizeof(double)*nathr ) ;
787
788 if( AFNI_yesenv("AFNI_CLUSTSIM_MEGA") ){ /* 15 Dec 2021 */
789 force_mega = 1 ; /* MEGA it is and MEGA it will be! */
790 do_niml = 1 ;
791 INFO_message("NOTE: AFNI_CLUSTSIM_MEGA forces use of -MEGA") ;
792 SET_MEGA ;
793 }
794
795 if( getenv("AFNI_RANDOM_SEEDVAL") != NULL ){
796 gseed = (unsigned int)AFNI_numenv("AFNI_RANDOM_SEEDVAL") ;
797 if( gseed == 0 )
798 gseed = ((unsigned int)time(NULL)) + 17*(unsigned int)getpid() ;
799 INFO_message("random seed set to %u from AFNI_RANDOM_SEEDVAL",gseed) ;
800 }
801
802 while( nopt < argc ){
803
804 /*----- -inset iii -----*/
805
806 if( strcmp(argv[nopt],"-inset") == 0 ){ /* 02 Feb 2016 */
807 int ii,nbad=0 ; THD_3dim_dataset *qset ;
808 if( num_inset > 0 )
809 ERROR_exit("You can't use '-inset' more than once!") ;
810 if( xinset != NULL )
811 ERROR_exit("You can't use '-inset' after '-insdat'!") ;
812 if( ++nopt >= argc )
813 ERROR_exit("You need at least 1 argument after option '-inset'") ;
814 for( ; nopt < argc && argv[nopt][0] != '-' ; nopt++ ){
815 /* ININFO_message("opening -inset '%s'",argv[nopt]) ; */
816 qset = THD_open_dataset(argv[nopt]) ;
817 if( qset == NULL ){
818 ERROR_message("-inset '%s': failure to open dataset",argv[nopt]) ;
819 nbad++ ; continue ;
820 }
821 for( ii=0 ; ii < DSET_NVALS(qset) ; ii++ ){
822 if( DSET_BRICK_TYPE(qset,ii) != MRI_float ){
823 ERROR_message("-inset '%s': all sub-bricks must be float :-(",argv[nopt]) ;
824 nbad++ ; break ;
825 }
826 }
827 if( num_inset > 0 && DSET_NVOX(qset) != DSET_NVOX(inset[0]) ){
828 ERROR_message("-inset '%s': grid size doesn't match other datasets",argv[nopt]) ;
829 nbad++ ;
830 }
831 inset = (THD_3dim_dataset **)realloc( inset , sizeof(THD_3dim_dataset *)*(num_inset+1)) ;
832 nval_inset = (int *) realloc( nval_inset, sizeof(int) *(num_inset+1)) ;
833 inset[num_inset] = qset ; nval_inset[num_inset] = DSET_NVALS(qset) ; num_inset++ ;
834 }
835 if( num_inset == 0 ) ERROR_exit("no valid datasets opened after -inset :-(") ;
836 if( nbad > 0 ) ERROR_exit("can't continue after above -inset problems") ;
837 continue ;
838 }
839
840 /*----- -insdat mask sdata -----*/
841
842 if( strcasecmp(argv[nopt],"-insdat") == 0 ){
843 int nfile , ii ;
844 if( xinset != NULL )
845 ERROR_exit("You can't use option '%s' more than once!",argv[nopt]) ;
846 if( num_inset > 0 )
847 ERROR_exit("You can't use '-insdat' after '-inset'!") ;
848 if( ++nopt >= argc-1 )
849 ERROR_exit("You need at least 2 arguments after option '%s'",argv[nopt-1]) ;
850
851 for( ii=nopt ; ii < argc && argv[ii][0] != '-' ; ii++ ) ; /*nada*/
852 nfile = ii-nopt ;
853
854 if( verb )
855 INFO_message("Loading %s datasets",argv[nopt-1]) ;
856
857 xinset = open_Xdataset( argv[nopt], nfile-1,argv+(nopt+1) ) ;
858
859 nopt += nfile ; continue ;
860 }
861
862 /*----- -nxyz n1 n2 n3 -----*/
863
864 if( strcmp(argv[nopt],"-nxyz") == 0 ){
865 nopt++ ; if( nopt+2 >= argc ) ERROR_exit("need 3 arguments after -nxyz") ;
866 nx = (int)strtod(argv[nopt++],NULL); if( nx <= 0 ) ERROR_exit("illegal nx value") ;
867 ny = (int)strtod(argv[nopt++],NULL); if( ny <= 0 ) ERROR_exit("illegal ny value") ;
868 nz = (int)strtod(argv[nopt++],NULL); if( nz <= 0 ) ERROR_exit("illegal nz value") ;
869 continue ;
870 }
871
872 /*----- -dxyz d1 d2 d3 -----*/
873
874 if( strcmp(argv[nopt],"-dxyz") == 0 ){
875 nopt++ ; if( nopt+2 >= argc ) ERROR_exit("need 3 arguments after -dxyz") ;
876 dx = strtod(argv[nopt++],NULL); if( dx <= 0.0f ) ERROR_exit("illegal dx value") ;
877 dy = strtod(argv[nopt++],NULL); if( dy <= 0.0f ) ERROR_exit("illegal dy value") ;
878 dz = strtod(argv[nopt++],NULL); if( dz <= 0.0f ) ERROR_exit("illegal dz value") ;
879 continue ;
880 }
881
882 /**** -OKsmallmask ****/
883
884 if( strcmp(argv[nopt],"-OKsmallmask") == 0 ){ /* 29 Mar 2011 */
885 minmask = 2 ; nopt++ ; continue ;
886 }
887
888 /**** -tdof dof ****/
889
890 if( strcmp(argv[nopt],"-tdof") == 0 ){ /* 26 May 2015 -- hidden option */
891 nopt++ ; if( nopt >= argc ) ERROR_exit("need argument after -tdof!") ;
892 tdof = (float)strtod(argv[nopt],NULL) ;
893 if( tdof > 0.0f && tdof < 4.0f ) ERROR_exit("illegal value after -tdof") ;
894 nopt++ ; continue ;
895 }
896
897 /**** -mask mset ****/
898
899 if( strcmp(argv[nopt],"-mask") == 0 ){
900 if( mask_dset != NULL ) ERROR_exit("Can't use -mask twice!") ;
901 nopt++ ; if( nopt >= argc ) ERROR_exit("need argument after -mask!") ;
902 mask_dset = THD_open_dataset(argv[nopt]);
903 if( mask_dset == NULL ) ERROR_exit("can't open -mask dataset!") ;
904 mask_vol = THD_makemask( mask_dset , 0 , 1.0,0.0 ) ;
905 if( mask_vol == NULL ) ERROR_exit("can't use -mask dataset!") ;
906 mask_nvox = DSET_NVOX(mask_dset) ;
907 DSET_unload(mask_dset) ;
908 mask_ngood = THD_countmask( mask_nvox , mask_vol ) ;
909 if( mask_ngood < minmask ){
910 if( minmask > 2 && mask_ngood > 2 ){
911 ERROR_message("-mask has only %d nonzero voxels; minimum allowed is %d.",
912 mask_ngood , minmask ) ;
913 ERROR_message("To run this simulation, please read the CAUTION and CAVEAT in -help,") ;
914 ERROR_message("and then you can use the '-OKsmallmask' option if you so desire.") ;
915 ERROR_exit("Cannot continue -- may we meet under happier circumstances!") ;
916 } else if( mask_ngood == 0 ){
917 ERROR_exit("-mask has no nonzero voxels -- cannot use this at all :-(") ;
918 } else {
919 ERROR_exit("-mask has only %d nonzero voxel%s -- cannot use this :-(",
920 mask_ngood , (mask_ngood > 1) ? "s" : "\0" ) ;
921 }
922 }
923 if( verb ) INFO_message("%d voxels in mask (%.2f%% of total)",
924 mask_ngood,100.0*mask_ngood/(double)mask_nvox) ;
925 nopt++ ; continue ;
926 }
927
928 /*----- -fwhm s -----*/
929
930 if( strcasecmp(argv[nopt],"-fwhm") == 0 ){
931 nopt++; if( nopt >= argc ) ERROR_exit("need argument after -fwhm");
932 fwhm_x = (float)strtod(argv[nopt],&ep) ;
933 if( fwhm_x < 0.0f || ep == argv[nopt] )
934 ERROR_exit("illegal value after -fwhm") ;
935 WARNING_message(
936 "Use of -fwhm is not advised; consider using -acf instead\n"
937 " FMRI data does not have Gaussian-shaped smoothness!" ) ;
938 fwhm_y = fwhm_z = fwhm_x ; nopt++; continue;
939 }
940
941 if( strncasecmp(argv[nopt],"-classic",6) == 0 ){ /* 01 Dec 2015 */
942 do_classic = 1 ; nopt++ ; continue ;
943 }
944
945 /*----- -acf a b c [30 Nov 2015] -----*/
946
947 if( strcasecmp(argv[nopt],"-acf") == 0 ){
948 nopt++; if( nopt+2 >= argc ) ERROR_exit("need 3 arguments after -acf");
949 acf_a = (float)strtod(argv[nopt++],&ep) ;
950 if( acf_a < 0.0f || acf_a > 1.0f || ep == argv[nopt-1] )
951 ERROR_exit("-acf: 'a' value should be between 0 and 1 :-(") ;
952 acf_b = (float)strtod(argv[nopt++],&ep) ;
953 if( acf_b <= 0.0f || ep == argv[nopt-1] )
954 ERROR_exit("-acf: 'b' value should be positive :-(") ;
955 acf_c = (float)strtod(argv[nopt++],&ep) ;
956 if( acf_c <= 0.0f || ep == argv[nopt-1] )
957 ERROR_exit("-acf: 'c' value should be positive :-(") ;
958 do_acf = 1 ; continue ;
959 }
960
961 /*----- -nopad -----*/
962
963 if( strcasecmp(argv[nopt],"-nopad") == 0 ){ /* 12 May 2015 */
964 allow_padding = 0 ; nopt++ ; continue ;
965 }
966
967 /*----- -fwhmxyz s s s -----*/
968
969 if( strcasecmp(argv[nopt],"-fwhmxyz") == 0 ){
970 nopt++; if(nopt+2 >= argc) ERROR_exit("need 3 arguments after -fwhmxyz");
971 fwhm_x = (float)strtod(argv[nopt++],&ep) ;
972 if( fwhm_x < 0.0f || ep == argv[nopt-1] )
973 ERROR_exit("illegal x value after -fwhmxyz") ;
974 fwhm_y = (float)strtod(argv[nopt++],&ep) ;
975 if( fwhm_y < 0.0f || ep == argv[nopt-1] )
976 ERROR_exit("illegal y value after -fwhmxyz") ;
977 fwhm_z = (float)strtod(argv[nopt++],&ep) ;
978 if( fwhm_z < 0.0f || ep == argv[nopt-1] )
979 ERROR_exit("illegal z value after -fwhmxyz") ;
980 WARNING_message(
981 "Use of -fwhmxyz is not advised; consider using -acf instead\n"
982 " FMRI data does not have Gaussian-shaped smoothness!" ) ;
983 continue;
984 }
985
986 /*----- -iter n -----*/
987
988 if( strcmp(argv[nopt],"-iter") == 0 || strncmp(argv[nopt],"-nite",5) == 0 ){
989 nopt++;
990 if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
991 niter = (int)strtod(argv[nopt],NULL) ;
992 if( niter < 2000 && !do_ssave ){
993 WARNING_message("-iter %d replaced by 2000",niter) ; niter = 2000 ;
994 } else if( niter < 1 ){
995 ERROR_exit("-iter %s is illegal :-(",argv[nopt]) ;
996 }
997 nopt++; continue;
998 }
999
1000 /*----- -seed S -----*/
1001
1002 if( strcasecmp(argv[nopt],"-seed") == 0 ){
1003 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
1004 gseed = (unsigned int)strtol(argv[nopt],NULL,10) ;
1005 if( gseed == 0 ){
1006 gseed = ((unsigned int)time(NULL)) + 17*(unsigned int)getpid() ;
1007 if( verb ) INFO_message("-seed 0 resets to %u",gseed) ;
1008 }
1009 nopt++; continue;
1010 }
1011
1012 /*----- -LOTS -----*/
1013
1014 if( strcasecmp(argv[nopt],"-LOTS") == 0 ){
1015 if( !did_mega ){
1016 npthr = npthr_lots ;
1017 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ;
1018 memcpy( pthr , pthr_lots , sizeof(double)*npthr ) ;
1019 nathr = nathr_lots ;
1020 athr = (double *)realloc(athr,sizeof(double)*nathr) ;
1021 memcpy( athr , athr_lots , sizeof(double)*nathr ) ;
1022 athr_sum_bot = 10 ; athr_sum_top = 22 ;
1023 }
1024 nopt++ ; continue ;
1025 }
1026
1027 /*----- -MEGA -----*/
1028
1029 if( strcasecmp(argv[nopt],"-MEGA") == 0 ){ /* 18 Dec 2015 */
1030 if( !did_mega ){
1031 npthr = npthr_mega ;
1032 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ;
1033 memcpy( pthr , pthr_mega , sizeof(double)*npthr ) ;
1034 nathr = nathr_mega ;
1035 athr = (double *)realloc(athr,sizeof(double)*nathr) ;
1036 memcpy( athr , athr_mega , sizeof(double)*nathr ) ;
1037 athr_sum_bot = 13 ; athr_sum_top = 31 ;
1038 if( niter < 30000 ) niter = 30000 ;
1039 }
1040 nopt++ ; continue ;
1041 }
1042
1043 /*----- -sumup -----*/
1044
1045 if( strcasecmp(argv[nopt],"-sumup") == 0 ){ /* 18 Dec 2015 */
1046 do_athr_sum++ ; nopt++ ; continue ;
1047 }
1048
1049 #if 0
1050 /*----- -2sided -----*/
1051
1052 if( strcasecmp(argv[nopt],"-2sided") == 0 ){
1053 do_2sid = 1 ; nopt++ ; continue ;
1054 }
1055 if( strcasecmp(argv[nopt],"-1sided") == 0 ){
1056 do_2sid = 0 ; nopt++ ; continue ;
1057 }
1058 #endif
1059
1060 /*----- -pthr p -----*/
1061
1062 if( strcmp(argv[nopt],"-pthr") == 0 || strcmp(argv[nopt],"-pval") == 0 ){
1063 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
1064 if( force_mega ){
1065 WARNING_message("AFNI_CLUSTSIM_MEGA is YES -- skipping '-pthr' option") ;
1066 for( ii=nopt ; ii < argc && argv[ii][0] != '-' ; ii++ ) ; /*nada*/
1067 continue ;
1068 }
1069 if( did_mega )
1070 WARNING_message("-pthr option is over-riding earlier -MEGA option!") ;
1071 if( strcmp(argv[nopt],"LOTS") == 0 ){ /* built in big table of values */
1072 npthr = npthr_lots ;
1073 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ;
1074 memcpy( pthr , pthr_lots , sizeof(double)*npthr ) ;
1075 nopt++ ;
1076 } else { /* a list of values */
1077 for( ii=nopt ; ii < argc && argv[ii][0] != '-' ; ii++ ) ; /*nada*/
1078 npthr = ii-nopt ;
1079 if( npthr <= 0 ) ERROR_exit("No positive values found after %s",argv[nopt-1]) ;
1080 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ;
1081 for( ii=0 ; ii < npthr ; ii++ ){
1082 pthr[ii] = strtod(argv[nopt+ii],NULL) ;
1083 if( pthr[ii] <= 0.0 || pthr[ii] > PMAX )
1084 ERROR_exit("value '%s' after '%s' is illegal!",argv[nopt+ii],argv[nopt-1]) ;
1085 }
1086 if( npthr > 1 ){
1087 for( ii=0 ; ii < npthr ; ii++ ) pthr[ii] = -pthr[ii] ; /* sort into */
1088 qsort_double( npthr , pthr ) ; /* descending */
1089 for( ii=0 ; ii < npthr ; ii++ ) pthr[ii] = -pthr[ii] ; /* order */
1090 for( ii=1 ; ii < npthr ; ii++ ){
1091 if( pthr[ii] == pthr[ii-1] )
1092 WARNING_message("duplicate value %g after '%s'",pthr[ii],argv[nopt-1]) ;
1093 }
1094 }
1095 nopt += npthr ;
1096 }
1097 have_pthr = 1 ; /* warn if -pthr is before -both/-niml */
1098 continue ;
1099 }
1100
1101 /*----- -athr p -----*/
1102
1103 if( strcmp(argv[nopt],"-athr") == 0 || strcmp(argv[nopt],"-aval") == 0 ){
1104 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
1105 if( force_mega ){
1106 WARNING_message("AFNI_CLUSTSIM_MEGA is YES -- skipping '-athr' option") ;
1107 for( ii=nopt ; ii < argc && argv[ii][0] != '-' ; ii++ ) ; /*nada*/
1108 continue ;
1109 }
1110 if( did_mega )
1111 WARNING_message("-athr option is over-riding earlier -MEGA option!") ;
1112 if( strcmp(argv[nopt],"LOTS") == 0 ){ /* built in big table of values */
1113 nathr = nathr_lots ;
1114 athr = (double *)realloc(athr,sizeof(double)*nathr) ;
1115 memcpy( athr , athr_lots , sizeof(double)*nathr ) ;
1116 nopt++ ;
1117 } else { /* a list of values */
1118 for( ii=nopt ; ii < argc && argv[ii][0] != '-' ; ii++ ) ; /*nada*/
1119 nathr = ii-nopt ;
1120 if( nathr <= 0 ) ERROR_exit("No positive values found after %s",argv[nopt-1]) ;
1121 athr = (double *)realloc(athr,sizeof(double)*nathr) ;
1122 for( ii=0 ; ii < nathr ; ii++ ){
1123 athr[ii] = strtod(argv[nopt+ii],NULL) ;
1124 if( athr[ii] <= 0.0 || athr[ii] > PMAX )
1125 ERROR_exit("value '%s' after '%s' is illegal!",argv[nopt+ii],argv[nopt-1]) ;
1126 }
1127 if( nathr > 1 ){
1128 for( ii=0 ; ii < nathr ; ii++ ) athr[ii] = -athr[ii] ;
1129 qsort_double( nathr , athr ) ;
1130 for( ii=0 ; ii < nathr ; ii++ ) athr[ii] = -athr[ii] ;
1131 for( ii=1 ; ii < nathr ; ii++ ){
1132 if( athr[ii] == athr[ii-1] )
1133 WARNING_message("duplicate value %g after '%s'",athr[ii],argv[nopt-1]) ;
1134 }
1135 }
1136 nopt += nathr ;
1137 }
1138 continue ;
1139 }
1140
1141 /*---- -nodec ----*/
1142
1143 if( strcasecmp(argv[nopt],"-nodec") == 0 ){
1144 nodec = 1 ; nopt++ ; continue ;
1145 }
1146
1147 #ifdef ALLOW_LOHI
1148 if( strcasecmp(argv[nopt],"-alo") == 0 ){ /* 20 Jan 2011: debug stuff */
1149 do_lohi = 1 ; nopt++ ; continue ;
1150 }
1151 #endif
1152
1153 /*---- -niml ----*/
1154
1155 if( strcasecmp(argv[nopt],"-niml") == 0 ||
1156 strcasecmp(argv[nopt],"-both") == 0 ){
1157 /* let users know what to do with -pthr 19 May 2015 [rickr] */
1158 if( have_pthr )
1159 ERROR_exit("-both or -niml cannot follow -pthr, as they override it\n"
1160 " i.e. -pthr should be after -both or -niml") ;
1161 if( !did_mega ){
1162 npthr = npthr_lots ;
1163 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ;
1164 memcpy( pthr , pthr_lots , sizeof(double)*npthr ) ;
1165 nathr = nathr_lots ;
1166 athr = (double *)realloc(athr,sizeof(double)*nathr) ;
1167 memcpy( athr , athr_lots , sizeof(double)*nathr ) ;
1168 }
1169 do_niml = 1 ;
1170 do_1D = ( strcasecmp(argv[nopt],"-both") == 0 ) ;
1171 nopt++ ; continue ;
1172 }
1173
1174 /*---- -BALL ----*/
1175
1176 if( strcasecmp(argv[nopt],"-BALL") == 0 ){
1177 do_ball = 1 ; nopt++ ; continue ;
1178 }
1179
1180 /*---- -NN ----*/
1181
1182 if( strcasecmp(argv[nopt],"-NN") == 0 ){
1183 WARNING_message("-NN option is no longer supported! All NN cases are computed now :)") ;
1184 nopt++ ; if( isdigit(argv[nopt][0]) ) nopt++ ;
1185 continue ;
1186 }
1187
1188 /*----- -prefix -----*/
1189
1190 if( strcmp(argv[nopt],"-prefix") == 0 ){
1191 nopt++ ; if( nopt >= argc ) ERROR_exit("need argument after -prefix!") ;
1192 prefix = strdup(argv[nopt]) ;
1193 if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix option!") ;
1194 nopt++ ; continue ;
1195 }
1196
1197 /*----- -prefix -----*/
1198
1199 if( strcmp(argv[nopt],"-cmd") == 0 ){
1200 nopt++ ; if( nopt >= argc ) ERROR_exit("need argument after -cmd!") ;
1201 cmd_fname = strdup(argv[nopt]) ;
1202 if( !THD_filename_ok(cmd_fname) ) ERROR_exit("bad -cmd option!") ;
1203 nopt++ ; continue ;
1204 }
1205
1206 /*---- -quiet ----*/
1207
1208 if( strcasecmp(argv[nopt],"-quiet") == 0 ){
1209 verb = 0 ; nopt++ ; continue ;
1210 }
1211
1212 /*----- -ssave -----*/
1213
1214 if( strncasecmp(argv[nopt],"-ssave:",7) == 0 ){
1215 if( nopt+1 >= argc ) ERROR_exit("need argument after %s",argv[nopt]) ;
1216 if( strcasecmp(argv[nopt]+7,"blurred") == 0 ){
1217 do_ssave = SSAVE_BLURRED ;
1218 } else if( strcasecmp(argv[nopt]+7,"masked") == 0 ){
1219 do_ssave = SSAVE_MASKED ;
1220 } else {
1221 ERROR_exit("this form of '-ssave' is unknown: '%s'",argv[nopt]) ;
1222 }
1223 ssave_prefix = strdup(argv[++nopt]) ;
1224 if( !THD_filename_ok(ssave_prefix) )
1225 ERROR_exit("bad prefix after option '%s'",argv[nopt-1]) ;
1226 nopt++ ; continue ;
1227 }
1228
1229 /*----- unknown option -----*/
1230
1231 ERROR_exit("3dClustSim -- unknown option '%s'",argv[nopt]) ;
1232 }
1233
1234 /*------- finalize some simple setup stuff --------*/
1235
1236 if( fwhm_x > 0.0f ){
1237 WARNING_message("**************************************************") ;
1238 WARNING_message(
1239 "I repeat:\n"
1240 " -fwhm uses a Gaussian-shaped autocorrelation function,\n"
1241 " which is not accurate for most FMRI data :(\n"
1242 " Consider using -acf instead, or 3dttest++ -Clustsim :)"
1243 ) ;
1244 WARNING_message("**************************************************") ;
1245 }
1246
1247 if( do_athr_sum && (athr_sum_bot < 0 || athr_sum_top < 0) ){ /* 18 Dec 2015 */
1248 do_athr_sum = 0 ;
1249 WARNING_message("-sumup canceled") ;
1250 }
1251
1252 if( xinset != NULL && num_inset > 0 ) /* should not be possible */
1253 ERROR_exit("-insdat and -inset are incompatible :-(") ;
1254
1255 #define INSET_PRELOAD
1256
1257 if( num_inset > 0 ){ /* 02 Feb 2016 */
1258 int qq,nbad=0 ;
1259 nx = DSET_NX(inset[0]) ;
1260 ny = DSET_NY(inset[0]) ;
1261 nz = DSET_NZ(inset[0]) ;
1262 dx = fabsf(DSET_DX(inset[0])) ;
1263 dy = fabsf(DSET_DY(inset[0])) ;
1264 dz = fabsf(DSET_DZ(inset[0])) ;
1265 for( niter=qq=0 ; qq < num_inset ; qq++ ){
1266 niter += nval_inset[qq] ;
1267 #ifdef INSET_PRELOAD
1268 ININFO_message("loading -inset '%s' with %d volumes",DSET_HEADNAME(inset[qq]),DSET_NVALS(inset[qq])) ;
1269 DSET_load(inset[qq]) ;
1270 if( !DSET_LOADED(inset[qq]) ){
1271 ERROR_message("Can't load dataset -inset '%s'",DSET_HEADNAME(inset[qq])) ;
1272 nbad ++ ;
1273 }
1274 #endif
1275 }
1276 if( nbad > 0 ) ERROR_exit("Can't continue without all -inset datasets :-(") ;
1277 if( niter < 100 )
1278 WARNING_message("-inset has only %d volumes total (= new '-niter' value)",niter) ;
1279 else if( verb )
1280 INFO_message("-inset had %d volumes = new '-niter' value",niter) ;
1281 if( mask_dset != NULL ){
1282 if( nx != DSET_NX(mask_dset) ||
1283 ny != DSET_NY(mask_dset) ||
1284 nz != DSET_NZ(mask_dset) )
1285 ERROR_exit("-mask and -inset don't match in grid dimensions :-(") ;
1286 }
1287 if( do_ssave ){ WARNING_message("-inset turns off -ssave"); do_ssave = 0; }
1288 if( do_acf ){ WARNING_message("-inset turns off -acf") ; do_acf = 0 ; }
1289
1290 } else if( xinset != NULL ){ /*--- 29 Aug 2016 ---*/
1291 mask_dset = xinset->mask_dset ;
1292 mask_ngood = xinset->ngood ;
1293 mask_vol = xinset->mask_vol ;
1294 nx = DSET_NX(mask_dset) ;
1295 ny = DSET_NY(mask_dset) ;
1296 nz = DSET_NZ(mask_dset) ;
1297 dx = fabsf(DSET_DX(mask_dset)) ;
1298 dy = fabsf(DSET_DY(mask_dset)) ;
1299 dz = fabsf(DSET_DZ(mask_dset)) ;
1300 niter = xinset->nvtot ;
1301 if( niter < 100 )
1302 WARNING_message("-insdat has only %d volumes total (= new '-niter' value)",niter) ;
1303 else if( verb )
1304 INFO_message("-insdat had %d volumes = new '-niter' value",niter) ;
1305 if( do_ssave ){ WARNING_message("-insdat turns off -ssave"); do_ssave = 0; }
1306 if( do_acf ){ WARNING_message("-insdat turns off -acf") ; do_acf = 0; }
1307
1308 } else if( mask_dset != NULL ){
1309 nx = DSET_NX(mask_dset) ;
1310 ny = DSET_NY(mask_dset) ;
1311 nz = DSET_NZ(mask_dset) ;
1312 dx = fabsf(DSET_DX(mask_dset)) ;
1313 dy = fabsf(DSET_DY(mask_dset)) ;
1314 dz = fabsf(DSET_DZ(mask_dset)) ;
1315 if( do_ssave > 0 ){ /* 24 Apr 2014 */
1316 ssave_dset = EDIT_empty_copy( mask_dset ) ;
1317 EDIT_dset_items( ssave_dset ,
1318 ADN_nvals , 1 ,
1319 ADN_prefix , "RandomJunk" ,
1320 ADN_none ) ;
1321 }
1322 }
1323
1324 if( do_acf ){ /* 30 Nov 2015 */
1325 float val ; int_triple ijk ;
1326 if( fwhm_x > 0.0f )
1327 WARNING_message("-acf option ==> any -fwhm options are ignored!") ;
1328 if( nx < 4 || ny < 4 || nz < 4 )
1329 ERROR_exit("-acf option: minimum grid is 4x4x4, but you have %dx%dx%d",nx,ny,nz) ;
1330 fwhm_x = fwhm_y = fwhm_z = 0.0f ;
1331 acf_parm[0] = acf_a ;
1332 acf_parm[1] = acf_b ;
1333 acf_parm[2] = acf_c ;
1334 val = 2.0f * rfunc_inv(0.5f,acf_parm) ;
1335 ijk = get_random_field_size( nx,ny,nz , dx,dy,dz , acf_parm ) ;
1336 acf_nxx = ijk.i ; acf_nyy = ijk.j ; acf_nzz = ijk.k ;
1337 INFO_message("ACF(%.2f,%.2f,%.2f) => FWHM=%.2f => %dx%dx%d pads to %dx%dx%d",
1338 acf_a,acf_b,acf_c,val,
1339 nx,ny,nz , acf_nxx,acf_nyy,acf_nzz ) ;
1340 acf_wim = make_radial_weight( acf_nxx,acf_nyy,acf_nzz , dx,dy,dz , acf_parm ) ;
1341 acf_ntar = acf_nxx+acf_nyy+acf_nzz ;
1342 }
1343
1344 if( do_ssave > 0 && ssave_dset == NULL ){ /* 24 Apr 2014 */
1345 char gstr[128] ; float xorg,yorg,zorg ;
1346 xorg = -0.5*dx*(nx-1); yorg = -0.5*dy*(ny-1); zorg = -0.5*dz*(nz-1);
1347 sprintf(gstr,"RAI:D:%d,%f,%f,%d,%f,%f,%d,%f,%f",
1348 nx,xorg,dx , ny,yorg,dy , nz,zorg,dz ) ;
1349 ssave_dset = EDIT_geometry_constructor(gstr,"RandomJunque") ;
1350 }
1351
1352 nxy = nx*ny ; nxyz = nxy*nz ; nxyz1 = nxyz - nxy ;
1353 if( nxyz < 256 )
1354 ERROR_exit("Only %d voxels in simulation?! Need at least 256.",nxyz) ;
1355
1356 if( mask_vol == NULL && do_ball ){
1357 float xq,yq,zq , nx2,ny2,nz2 , nxq,nyq,nzq ; int ii,jj,kk ;
1358
1359 nx2 = 0.5f*(nx-1) ; nxq = nx2*nx2 + 0.501f ;
1360 ny2 = 0.5f*(ny-1) ; nyq = ny2*ny2 + 0.501f ;
1361 nz2 = 0.5f*(nz-1) ; nzq = nz2*nz2 + 0.501f ;
1362 mask_vol = (byte *)malloc(sizeof(byte)*nxyz) ;
1363 for( kk=0 ; kk < nz ; kk++ ){
1364 zq = kk - nz2 ; zq = (zq*zq) / nzq ;
1365 for( jj=0 ; jj < ny ; jj++ ){
1366 yq = jj - ny2 ; yq = (yq*yq) / nyq ;
1367 for( ii=0 ; ii < nx ; ii++ ){
1368 xq = ii - nx2 ; xq = (xq*xq) / nxq ;
1369 mask_vol[ii+jj*nx+kk*nxy] = ((xq+yq+zq) <= 1.0f) ;
1370 }}}
1371 mask_nvox = nxyz ;
1372 mask_ngood = THD_countmask( mask_nvox , mask_vol ) ;
1373 if( verb ) INFO_message("%d voxels in BALL mask (%.1f%% of total)",
1374 mask_ngood,100.0*mask_ngood/(double)mask_nvox) ;
1375 }
1376
1377 if( mask_ngood == 0 ) mask_ngood = nxyz ;
1378
1379 srand48(gseed) ; /* not really needed */
1380
1381 /*-- z-score thresholds for the various p-values --*/
1382
1383 zthr_1sid = (float *)malloc(sizeof(float)*npthr) ;
1384 zthr_2sid = (float *)malloc(sizeof(float)*npthr) ;
1385 for( ii=0 ; ii < npthr ; ii++ ){
1386 zthr_1sid[ii] = (float)zthresh( pthr[ii] ) ;
1387 zthr_2sid[ii] = (float)zthresh( 0.5*pthr[ii] ) ;
1388 }
1389
1390 /*-- check niter vs smallest alpha level --*/
1391
1392 if( niter*athr[nathr-1] < 10.0f )
1393 WARNING_message(
1394 "Smallest athr=%g ==> should have at least %d iterations instead of %d" ,
1395 athr[nathr-1] , (int)(10.0f/athr[nathr-1])+1 , niter ) ;
1396
1397 /*-- blurring stuff --*/
1398
1399 sigmax = FWHM_TO_SIGMA(fwhm_x) ; if( nx < 2 ) sigmax = 0.0f ;
1400 sigmay = FWHM_TO_SIGMA(fwhm_y) ; if( ny < 2 ) sigmay = 0.0f ;
1401 sigmaz = FWHM_TO_SIGMA(fwhm_z) ; if( nz < 2 ) sigmaz = 0.0f ;
1402
1403 do_blur = (sigmax > 0.0f || sigmay > 0.0f || sigmaz > 0.0f ) ;
1404
1405 nx_pad = nx ; ny_pad = ny ; nz_pad = nz ;
1406 if( num_inset == 0 ){
1407 if( do_acf ){ /* 30 Nov 2015 */
1408 ex_pad = (acf_nxx-nx)/2 ;
1409 ey_pad = (acf_nyy-ny)/2 ;
1410 ez_pad = (acf_nzz-nz)/2 ;
1411 nx_pad = acf_nxx ;
1412 ny_pad = acf_nyy ;
1413 nz_pad = acf_nzz ;
1414 } else if( do_blur && allow_padding ){ /* 12 May 2015 */
1415 ex_pad = (int)rintf(1.666f*fwhm_x/dx) ;
1416 ey_pad = (int)rintf(1.666f*fwhm_y/dy) ;
1417 ez_pad = (int)rintf(1.666f*fwhm_z/dz) ;
1418 do_pad = (ex_pad > 0) || (ey_pad > 0) || (ez_pad > 0) ;
1419 nx_pad = nx + 2*ex_pad ;
1420 ny_pad = ny + 2*ey_pad ;
1421 nz_pad = nz + 2*ez_pad ;
1422 if( do_pad )
1423 INFO_message(
1424 "Padding by %d x %d x %d slices to allow for edge effects of blurring" ,
1425 ex_pad,ey_pad,ez_pad) ;
1426 }
1427 }
1428 nxy_pad = nx_pad * ny_pad ;
1429 nxyz_pad = nxy_pad * nz_pad ;
1430
1431 EXRETURN ;
1432 }
1433
1434 /*---------------------------------------------------------------------------*/
1435
ssave_dataset(float * fim)1436 void ssave_dataset( float *fim ) /* 24 Apr 2014 */
1437 {
1438 #pragma omp critical
1439 { char *spr , sbuf[16] ;
1440 sprintf(sbuf,"_%06d",ssave_index) ; ssave_index++ ;
1441 spr = modify_afni_prefix(ssave_prefix,NULL,sbuf) ;
1442 EDIT_substitute_brick(ssave_dset,0,MRI_float,fim) ;
1443 EDIT_dset_items(ssave_dset,ADN_prefix,spr,ADN_none) ;
1444 ssave_dset->idcode = MCW_new_idcode() ;
1445 DSET_write(ssave_dset) ;
1446 DSET_NULL_ARRAY(ssave_dset,0) ;
1447 }
1448 return ;
1449 }
1450
1451 /*---------------------------------------------------------------------------*/
1452 /* Create the "functional" image, from the xinset struct [29 Aug 2016] */
1453
generate_fim_insdat(float * fim,int ival)1454 void generate_fim_insdat( float *fim , int ival )
1455 {
1456 load_from_Xdataset( xinset , ival , fim ) ;
1457 return ;
1458 }
1459
1460 /*---------------------------------------------------------------------------*/
1461 /* Create the "functional" image, from the inset dataset [02 Feb 2016] */
1462
generate_fim_inset(float * fim,int ival)1463 void generate_fim_inset( float *fim , int ival )
1464 {
1465 if( ival < 0 || ival >= niter ){ /* should not be possible */
1466 ERROR_message("inset[%d] == out of range!",ival) ;
1467 memset( fim , 0 , sizeof(float)*nxyz ) ;
1468 } else {
1469 float *bar=NULL ; int ii,qq,qval ;
1470 for( qval=ival,qq=0 ; qq < num_inset ; qq++ ){ /* find which */
1471 if( qval < nval_inset[qq] ) break ; /* inset[qq] to use */
1472 qval -= nval_inset[qq] ;
1473 }
1474 if( qq == num_inset ){ /* should not be possible */
1475 ERROR_message("inset[%d] == array overflow !!",ival) ;
1476 } else {
1477 bar = DSET_ARRAY(inset[qq],qval) ;
1478 if( bar == NULL ){
1479 ININFO_message("loading -inset '%s' with %d volumes",DSET_HEADNAME(inset[qq]),DSET_NVALS(inset[qq])) ;
1480 DSET_load(inset[qq]) ; bar = DSET_ARRAY(inset[qq],qval) ;
1481 }
1482 }
1483 if( bar != NULL ){
1484 for( ii=0 ; ii < nxyz ; ii++ ) fim[ii] = bar[ii] ; /* copy data */
1485 DSET_unload_one(inset[qq],qval) ;
1486 } else {
1487 ERROR_message("inset[%d] == NULL :-(",ival) ;
1488 memset( fim , 0 , sizeof(float)*nxyz ) ;
1489 }
1490 }
1491 }
1492
1493 /*---------------------------------------------------------------------------*/
1494 /* Create the "functional" image, the ACF way [30 Nov 2015]. */
1495
generate_fim_acf(float * fim,unsigned short xran[])1496 void generate_fim_acf( float *fim , unsigned short xran[] )
1497 {
1498 DECLARE_ithr ;
1499 int ii,jj,kk,pp,qq ;
1500 float *afim ; int do_a ;
1501
1502 /* the function generates images in pairs, but 3dClustSim only
1503 consumes them one at a time -- thus the artifice of the 'aim'
1504 and 'bim' images, which are created together, then aim is
1505 consumed right away and bim is consumed on the next call here. */
1506
1507 if( acf_aim[ithr] == NULL && acf_bim[ithr] == NULL ){ /* need new aim & bim */
1508 MRI_IMARR *imar ;
1509 imar = make_radial_random_field( acf_nxx , acf_nyy , acf_nzz ,
1510 acf_wim , acf_tar[ithr] , xran ) ;
1511 acf_aim[ithr] = IMARR_SUBIM(imar,0) ;
1512 acf_bim[ithr] = IMARR_SUBIM(imar,1) ;
1513 FREE_IMARR(imar) ;
1514 }
1515
1516 do_a = (acf_aim[ithr] != NULL) ;
1517 afim = (do_a) ? MRI_FLOAT_PTR(acf_aim[ithr]) : MRI_FLOAT_PTR(acf_bim[ithr]) ;
1518
1519 /* cut back to smaller unpadded volume */
1520
1521 for( pp=kk=0 ; kk < nz ; kk++ ){
1522 for( jj=0 ; jj < ny ; jj++ ){
1523 qq = ex_pad + (jj+ey_pad)*nx_pad + (kk+ez_pad)*nxy_pad ;
1524 for( ii=0 ; ii < nx ; ii++ ) fim[pp++] = afim[qq++] ;
1525 }}
1526
1527 /* free and clear whichever image we consumed here */
1528
1529 if( do_a ){
1530 mri_free(acf_aim[ithr]) ; acf_aim[ithr] = NULL ;
1531 } else {
1532 mri_free(acf_bim[ithr]) ; acf_bim[ithr] = NULL ;
1533 }
1534
1535 return ;
1536 }
1537
1538 /*---------------------------------------------------------------------------*/
1539 /* Create the "functional" image, with smoothing and padding [12 May 2015] */
1540
generate_fim_padded(float * fim,float * pfim,unsigned short xran[])1541 void generate_fim_padded( float *fim , float *pfim , unsigned short xran[] )
1542 {
1543 int ii,jj,kk,pp,qq ;
1544
1545 /* random N(0,1) stuff to create the larger (padded) volume */
1546
1547 for( ii=0 ; ii < nxyz_pad ; ii++ ) pfim[ii] = zgaussian_sss(xran) ;
1548
1549 /* smoothization */
1550
1551 if( do_blur )
1552 FIR_blur_volume_3d(nx_pad,ny_pad,nz_pad,dx,dy,dz,pfim,sigmax,sigmay,sigmaz) ;
1553
1554 /* cut back to smaller unpadded volume */
1555
1556 for( pp=kk=0 ; kk < nz ; kk++ ){
1557 for( jj=0 ; jj < ny ; jj++ ){
1558 qq = ex_pad + (jj+ey_pad)*nx_pad + (kk+ez_pad)*nxy_pad ;
1559 for( ii=0 ; ii < nx ; ii++ ) fim[pp++] = pfim[qq++] ;
1560 }}
1561
1562 return ;
1563 }
1564
1565 /*---------------------------------------------------------------------------*/
1566 /* Create the functional "image", with no padding and optional smoothing. */
1567
generate_fim_unpadded(float * fim,unsigned short xran[])1568 void generate_fim_unpadded( float *fim , unsigned short xran[] )
1569 {
1570 int ii ; float sum ;
1571
1572 /* random N(0,1) stuff */
1573
1574 for( ii=0 ; ii < nxyz ; ii++ ) fim[ii] = zgaussian_sss(xran) ;
1575
1576 /* smoothization */
1577
1578 if( do_blur )
1579 FIR_blur_volume_3d(nx,ny,nz,dx,dy,dz,fim,sigmax,sigmay,sigmaz) ;
1580
1581 return ;
1582 }
1583
1584 /*---------------------------------------------------------------------------*/
1585 /* Generate random smoothed masked image, with stdev=1. */
1586
generate_image(float * fim,float * pfim,unsigned short xran[],int iter)1587 void generate_image( float *fim , float *pfim , unsigned short xran[] , int iter )
1588 {
1589 register int ii ; register float sum ;
1590
1591 /* Outsource the creation of the smoothed random field [12 May 2015] */
1592
1593 if( inset != NULL )
1594 generate_fim_inset( fim , iter-1 ) ;
1595 else if( xinset != NULL )
1596 generate_fim_insdat( fim , iter-1 ) ;
1597 else if( do_acf )
1598 generate_fim_acf( fim , xran ) ;
1599 else if( do_pad )
1600 generate_fim_padded( fim , pfim , xran ) ;
1601 else
1602 generate_fim_unpadded( fim , xran ) ;
1603
1604 /* normalizing */
1605
1606 if( num_inset == 0 && xinset == NULL ){
1607 for( sum=0.0f,ii=0 ; ii < nxyz ; ii++ ) sum += fim[ii]*fim[ii] ;
1608 sum = sqrtf( nxyz / sum ) ; /* fix stdev back to 1 */
1609 for( ii=0 ; ii < nxyz ; ii++ ) fim[ii] *= sum ;
1610 }
1611
1612 /* save this volume? */
1613
1614 if( do_ssave == SSAVE_BLURRED ) ssave_dataset(fim) ;
1615
1616 /* maskizing? */
1617
1618 if( mask_vol != NULL && xinset == NULL ){
1619 for( ii=0 ; ii < nxyz ; ii++ ) if( !mask_vol[ii] ) fim[ii] = 0.0f ;
1620 }
1621
1622 if( num_inset == 0 && xinset == NULL && tdof > 0.0f ){
1623 float zfac = 1.0f/(1.0f-0.25f/tdof) ; /* 26 May 2015: secret stuff */
1624 float tfac = 0.5f/tdof ;
1625 float zhat , denom ;
1626 for( ii=0 ; ii < nxyz ; ii++ ){
1627 if( fim[ii] != 0.0f ){
1628 zhat = fim[ii]*zfac ;
1629 denom = 1.0f - zhat*zhat*tfac ; if( denom < 0.1f ) denom = 0.1f ;
1630 fim[ii] = zhat / sqrtf(denom) ;
1631 }
1632 }
1633 }
1634
1635 /* save THIS volume? */
1636
1637 if( do_ssave == SSAVE_MASKED ) ssave_dataset(fim) ;
1638
1639 #if 0
1640 /* absolution? */
1641
1642 if( do_2sid ){
1643 for( ii=0 ; ii < nxyz ; ii++ ) fim[ii] = fabsf(fim[ii]) ;
1644 }
1645 #endif
1646
1647 return ;
1648 }
1649
1650 #ifdef USE_SHAVE
1651 /*---------------------------------------------------------------------------*/
1652 /* 'shave' == 'short save' == saving generated images
1653 for re-use in the sumup phase of the program.
1654 Only the "in-mask" part is saved, and values are scale to shorts.
1655 Memory is malloc-ed for small needs, and is mmap-ed for large needs.
1656 *//*-------------------------------------------------------------------------*/
1657
1658 #include <sys/mman.h>
1659 #if !defined(MAP_ANON) && defined(MAP_ANONYMOUS)
1660 # define MAP_ANON MAP_ANONYMOUS
1661 #endif
1662
setup_shave(void)1663 void setup_shave(void)
1664 {
1665 int64_t twogig = 2ll * 1024ll * 1024ll * 1024ll ;
1666
1667 ENTRY("setup_shave") ;
1668
1669 shave_siz = (size_t)mask_ngood ; /* in units of sizeof(short) */
1670 shave_tot = shave_siz * (int64_t)niter * sizeof(short) ;
1671
1672 /* check is system is 32-bit and memory needed is over 2G */
1673
1674 if( shave_tot >= twogig &&
1675 ( sizeof(void *) < 8 || sizeof(size_t) < 8 ) )
1676 ERROR_exit("Total space needed for internal save of simulations\n"
1677 " exceeds 2 GB -- cannot proceed on a 32-bit system!") ;
1678
1679 #ifndef MAP_ANON
1680 do_shave = SHAVE_MALLOC
1681 #else
1682 do_shave = (shave_tot >= twogig) ? SHAVE_MMAP : SHAVE_MALLOC ;
1683 #endif
1684
1685 if( do_shave == SHAVE_MALLOC ){
1686 shave = (short *)malloc((size_t)shave_tot) ;
1687 } else {
1688 #ifdef MAP_ANON
1689 shave = mmap( (void *)0 , (size_t)shave_tot ,
1690 PROT_READ | PROT_WRITE , MAP_ANON | MAP_SHARED , -1,0 ) ;
1691 #endif
1692 }
1693
1694 if( shave == NULL )
1695 ERROR_exit("Cannot allocate space for internal save of simulations :-(") ;
1696
1697 if( verb )
1698 INFO_message("allocated %s (%s) bytes for sumup re-use",
1699 commaized_integer_string((long long)shave_tot) ,
1700 approximate_number_string((double)shave_tot) ) ;
1701
1702 EXRETURN ;
1703 }
1704
1705 /*---------------------------------------------------------------------------*/
1706 /* Trash the allocated shave space */
1707
destroy_shave(void)1708 void destroy_shave(void)
1709 {
1710 ENTRY("destroy_shave") ;
1711
1712 if( shave == NULL ) EXRETURN ;
1713
1714 if( do_shave == SHAVE_MALLOC ) free(shave) ;
1715 else munmap(shave,(size_t)shave_tot) ;
1716 shave = NULL ;
1717 EXRETURN ;
1718 }
1719
1720 /*---------------------------------------------------------------------------*/
1721 /* Actually save the iter-th image into the allocated space */
1722
1723 #define SHAVE_FAC 4681.0f /* 32767/7 */
1724 #define SHAVE_INV 0.0002136296f /* 7/32767 */
1725
fim_to_shave(float * fim,int iter)1726 void fim_to_shave( float *fim , int iter )
1727 {
1728 short *shar = shave + (size_t)(iter * shave_siz) ;
1729 int ii , jj ;
1730
1731 #if 0
1732 #pragma omp critical
1733 { fprintf(stderr,"f2s(%d) ",iter); }
1734 #endif
1735
1736 if( mask_vol != NULL ){
1737 for( jj=ii=0 ; ii < nxyz ; ii++ )
1738 if( mask_vol[ii] ) shar[jj++] = (short)(fim[ii]*SHAVE_FAC) ;
1739 } else {
1740 for( ii=0 ; ii < nxyz ; ii++ ) shar[ii] = (short)(fim[ii]*SHAVE_FAC) ;
1741 }
1742 return ;
1743 }
1744
1745 /*---------------------------------------------------------------------------*/
1746 /* Retrieve the iter-th image from the allocated space */
1747
shave_to_fim(float * fim,int iter)1748 void shave_to_fim( float *fim , int iter )
1749 {
1750 short *shar = shave + (size_t)(iter * shave_siz) ;
1751 int ii , jj ;
1752
1753 #if 0
1754 #pragma omp critical
1755 { fprintf(stderr,"s2f(%d) ",iter); }
1756 #endif
1757
1758 if( mask_vol != NULL ){
1759 for( jj=ii=0 ; ii < nxyz ; ii++ ){
1760 if( mask_vol[ii] ) fim[ii] = shar[jj++]*SHAVE_INV ;
1761 else fim[ii] = 0.0f ;
1762 }
1763 } else {
1764 for( ii=0 ; ii < nxyz ; ii++ ) fim[ii] = shar[ii]*SHAVE_INV ;
1765 }
1766 return ;
1767 }
1768 #endif
1769
1770 /*---------------------------------------------------------------------------*/
1771
1772 #define DALL MAX_CLUSTER_SIZE
1773
1774 /*! Put (i,j,k) into the current cluster, if it is nonzero. */
1775
1776 #define CPUT(i,j,k) \
1777 do{ ijk = THREE_TO_IJK(i,j,k,nx,nxy) ; \
1778 if( mmm[ijk] ){ \
1779 if( nnow == nall ){ /* increase array lengths */ \
1780 nall += DALL + nall/2 ; \
1781 inow = (short *) realloc(inow,sizeof(short)*nall) ; \
1782 jnow = (short *) realloc(jnow,sizeof(short)*nall) ; \
1783 know = (short *) realloc(know,sizeof(short)*nall) ; \
1784 } \
1785 inow[nnow] = i ; jnow[nnow] = j ; know[nnow] = k ; \
1786 nnow++ ; mmm[ijk] = 0 ; \
1787 } } while(0)
1788
1789 #define USE_MEMCHR
1790
1791 static int *nall_g = NULL ; /* per-thread workspaces */
1792 static short **inow_g = NULL ; /* for clusterizationing */
1793 static short **jnow_g = NULL ;
1794 static short **know_g = NULL ;
1795
1796 /*----------------------------------------------------------------------------*/
1797
find_largest_cluster_NN1(byte * mmm,int ithr)1798 int find_largest_cluster_NN1( byte *mmm , int ithr )
1799 {
1800 int max_size=0 ;
1801 int ii,jj,kk, icl , ijk , ijk_last ;
1802 int ip,jp,kp , im,jm,km ;
1803 int nnow , nall ; short *inow , *jnow , *know ;
1804 #ifdef USE_MEMCHR
1805 byte *mch ;
1806 #endif
1807
1808 /* get workspace for this thread */
1809
1810 nall = nall_g[ithr] ;
1811 inow = inow_g[ithr] ; jnow = jnow_g[ithr] ; know = know_g[ithr] ;
1812
1813 ijk_last = 0 ; /* start scanning at the start */
1814
1815 while(1) {
1816 /* find next nonzero point in mmm array */
1817
1818 #ifndef USE_MEMCHR
1819 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
1820 #else
1821 mch = memchr( mmm+ijk_last , 1 , nxyz-ijk_last ) ; /* quicker search */
1822 if( mch == NULL ) ijk = nxyz ;
1823 else ijk = mch - mmm ;
1824 #endif
1825 if( ijk == nxyz ) break ; /* didn't find any! */
1826 ijk_last = ijk+1 ; /* start here next time */
1827
1828 mmm[ijk] = 0 ; /* clear found point */
1829
1830 nnow = 1 ; /* # pts in cluster */
1831 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
1832
1833 /* loop over points in cluster, checking their neighbors,
1834 growing the cluster if we find any that belong therein */
1835
1836 for( icl=0 ; icl < nnow ; icl++ ){
1837 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
1838 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
1839 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
1840
1841 if( im >= 0 ) CPUT(im,jj,kk) ;
1842 if( ip < nx ) CPUT(ip,jj,kk) ;
1843 if( jm >= 0 ) CPUT(ii,jm,kk) ;
1844 if( jp < ny ) CPUT(ii,jp,kk) ;
1845 if( km >= 0 ) CPUT(ii,jj,km) ;
1846 if( kp < nz ) CPUT(ii,jj,kp) ;
1847 }
1848
1849 if( nnow > max_size ) max_size = nnow ;
1850 }
1851
1852 if( nall > nall_g[ithr] ){ /* probably won't happen */
1853 nall_g[ithr] = nall ;
1854 inow_g[ithr] = inow ; jnow_g[ithr] = jnow ; know_g[ithr] = know ;
1855 }
1856
1857 return max_size ;
1858 }
1859
1860 /*----------------------------------------------------------------------------*/
1861
find_largest_cluster_NN2(byte * mmm,int ithr)1862 int find_largest_cluster_NN2( byte *mmm , int ithr )
1863 {
1864 int max_size=0 ;
1865 int ii,jj,kk, icl , ijk , ijk_last ;
1866 int ip,jp,kp , im,jm,km ;
1867 int nnow , nall ; short *inow , *jnow , *know ;
1868 #ifdef USE_MEMCHR
1869 byte *mch ;
1870 #endif
1871
1872 /* get workspace for this thread */
1873
1874 nall = nall_g[ithr] ;
1875 inow = inow_g[ithr] ; jnow = jnow_g[ithr] ; know = know_g[ithr] ;
1876
1877 ijk_last = 0 ; /* start scanning at the start */
1878
1879 while(1) {
1880 /* find next nonzero point in mmm array */
1881
1882 #ifndef USE_MEMCHR
1883 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
1884 #else
1885 mch = memchr( mmm+ijk_last , 1 , nxyz-ijk_last ) ; /* quicker search */
1886 if( mch == NULL ) ijk = nxyz ;
1887 else ijk = mch - mmm ;
1888 #endif
1889 if( ijk == nxyz ) break ; /* didn't find any! */
1890 ijk_last = ijk+1 ; /* start here next time */
1891
1892 mmm[ijk] = 0 ; /* clear found point */
1893
1894 nnow = 1 ; /* # pts in cluster */
1895 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
1896
1897 /* loop over points in cluster, checking their neighbors,
1898 growing the cluster if we find any that belong therein */
1899
1900 for( icl=0 ; icl < nnow ; icl++ ){
1901 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
1902 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
1903 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
1904
1905 if( im >= 0 ){ CPUT(im,jj,kk) ;
1906 if( jm >= 0 ) CPUT(im,jm,kk) ; /* 2NN */
1907 if( jp < ny ) CPUT(im,jp,kk) ; /* 2NN */
1908 if( km >= 0 ) CPUT(im,jj,km) ; /* 2NN */
1909 if( kp < nz ) CPUT(im,jj,kp) ; /* 2NN */
1910 }
1911 if( ip < nx ){ CPUT(ip,jj,kk) ;
1912 if( jm >= 0 ) CPUT(ip,jm,kk) ; /* 2NN */
1913 if( jp < ny ) CPUT(ip,jp,kk) ; /* 2NN */
1914 if( km >= 0 ) CPUT(ip,jj,km) ; /* 2NN */
1915 if( kp < nz ) CPUT(ip,jj,kp) ; /* 2NN */
1916 }
1917 if( jm >= 0 ){ CPUT(ii,jm,kk) ;
1918 if( km >= 0 ) CPUT(ii,jm,km) ; /* 2NN */
1919 if( kp < nz ) CPUT(ii,jm,kp) ; /* 2NN */
1920 }
1921 if( jp < ny ){ CPUT(ii,jp,kk) ;
1922 if( km >= 0 ) CPUT(ii,jp,km) ; /* 2NN */
1923 if( kp < nz ) CPUT(ii,jp,kp) ; /* 2NN */
1924 }
1925 if( km >= 0 ) CPUT(ii,jj,km) ;
1926 if( kp < nz ) CPUT(ii,jj,kp) ;
1927 }
1928
1929 if( nnow > max_size ) max_size = nnow ;
1930 }
1931
1932 if( nall > nall_g[ithr] ){ /* probably won't happen */
1933 nall_g[ithr] = nall ;
1934 inow_g[ithr] = inow ; jnow_g[ithr] = jnow ; know_g[ithr] = know ;
1935 }
1936
1937 return max_size ;
1938 }
1939
1940 /*----------------------------------------------------------------------------*/
1941
find_largest_cluster_NN3(byte * mmm,int ithr)1942 int find_largest_cluster_NN3( byte *mmm , int ithr )
1943 {
1944 int max_size=0 ;
1945 int ii,jj,kk, icl , ijk , ijk_last ;
1946 int ip,jp,kp , im,jm,km ;
1947 int nnow , nall ; short *inow , *jnow , *know ;
1948 #ifdef USE_MEMCHR
1949 byte *mch ;
1950 #endif
1951
1952 /* get workspace for this thread */
1953
1954 nall = nall_g[ithr] ;
1955 inow = inow_g[ithr] ; jnow = jnow_g[ithr] ; know = know_g[ithr] ;
1956
1957 ijk_last = 0 ; /* start scanning at the start */
1958
1959 while(1) {
1960 /* find next nonzero point in mmm array */
1961
1962 #ifndef USE_MEMCHR
1963 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
1964 #else
1965 mch = memchr( mmm+ijk_last , 1 , nxyz-ijk_last ) ; /* quicker search */
1966 if( mch == NULL ) ijk = nxyz ;
1967 else ijk = mch - mmm ;
1968 #endif
1969 if( ijk == nxyz ) break ; /* didn't find any! */
1970 ijk_last = ijk+1 ; /* start here next time */
1971
1972 mmm[ijk] = 0 ; /* clear found point */
1973
1974 nnow = 1 ; /* # pts in cluster */
1975 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
1976
1977 /* loop over points in cluster, checking their neighbors,
1978 growing the cluster if we find any that belong therein */
1979
1980 for( icl=0 ; icl < nnow ; icl++ ){
1981 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
1982 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
1983 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
1984
1985 if( im >= 0 ){ CPUT(im,jj,kk) ;
1986 if( jm >= 0 ) CPUT(im,jm,kk) ; /* 2NN */
1987 if( jp < ny ) CPUT(im,jp,kk) ; /* 2NN */
1988 if( km >= 0 ) CPUT(im,jj,km) ; /* 2NN */
1989 if( kp < nz ) CPUT(im,jj,kp) ; /* 2NN */
1990 if( jm >= 0 && km >= 0 ) CPUT(im,jm,km) ; /* 3NN */
1991 if( jm >= 0 && kp < nz ) CPUT(im,jm,kp) ; /* 3NN */
1992 if( jp < ny && km >= 0 ) CPUT(im,jp,km) ; /* 3NN */
1993 if( jp < ny && kp < nz ) CPUT(im,jp,kp) ; /* 3NN */
1994 }
1995 if( ip < nx ){ CPUT(ip,jj,kk) ;
1996 if( jm >= 0 ) CPUT(ip,jm,kk) ; /* 2NN */
1997 if( jp < ny ) CPUT(ip,jp,kk) ; /* 2NN */
1998 if( km >= 0 ) CPUT(ip,jj,km) ; /* 2NN */
1999 if( kp < nz ) CPUT(ip,jj,kp) ; /* 2NN */
2000 if( jm >= 0 && km >= 0 ) CPUT(ip,jm,km) ; /* 3NN */
2001 if( jm >= 0 && kp < nz ) CPUT(ip,jm,kp) ; /* 3NN */
2002 if( jp < ny && km >= 0 ) CPUT(ip,jp,km) ; /* 3NN */
2003 if( jp < ny && kp < nz ) CPUT(ip,jp,kp) ; /* 3NN */
2004 }
2005 if( jm >= 0 ){ CPUT(ii,jm,kk) ;
2006 if( km >= 0 ) CPUT(ii,jm,km) ; /* 2NN */
2007 if( kp < nz ) CPUT(ii,jm,kp) ; /* 2NN */
2008 }
2009 if( jp < ny ){ CPUT(ii,jp,kk) ;
2010 if( km >= 0 ) CPUT(ii,jp,km) ; /* 2NN */
2011 if( kp < nz ) CPUT(ii,jp,kp) ; /* 2NN */
2012 }
2013 if( km >= 0 ) CPUT(ii,jj,km) ;
2014 if( kp < nz ) CPUT(ii,jj,kp) ;
2015 }
2016
2017 if( nnow > max_size ) max_size = nnow ;
2018 }
2019
2020 if( nall > nall_g[ithr] ){ /* probably won't happen */
2021 nall_g[ithr] = nall ;
2022 inow_g[ithr] = inow ; jnow_g[ithr] = jnow ; know_g[ithr] = know ;
2023 }
2024
2025 return max_size ;
2026 }
2027
2028 /*---------------------------------------------------------------------------*/
2029 /* Find clusters, save some info */
2030
gather_stats_NN1_1sid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2031 void gather_stats_NN1_1sid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2032 {
2033 register int ii ; register float thr ; int siz ;
2034
2035 thr = zthr_1sid[ipthr] ;
2036 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2037 siz = find_largest_cluster_NN1( bfim , ithr ) ;
2038 if( siz > max_cluster_size ) siz = max_cluster_size ;
2039 mtab[siz]++ ; csiz_1sid_NN1[ipthr][iter-1] = siz ;
2040
2041 return ;
2042 }
2043
2044 /*---------------------------------------------------------------------------*/
2045 /* Find clusters, save some info */
2046
gather_stats_NN2_1sid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2047 void gather_stats_NN2_1sid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2048 {
2049 register int ii ; register float thr ; int siz ;
2050
2051 thr = zthr_1sid[ipthr] ;
2052 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2053 siz = find_largest_cluster_NN2( bfim , ithr ) ;
2054 if( siz > max_cluster_size ) siz = max_cluster_size ;
2055 mtab[siz]++ ; csiz_1sid_NN2[ipthr][iter-1] = siz ;
2056
2057 return ;
2058 }
2059
2060 /*---------------------------------------------------------------------------*/
2061 /* Find clusters, save some info */
2062
gather_stats_NN3_1sid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2063 void gather_stats_NN3_1sid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2064 {
2065 register int ii ; register float thr ; int siz ;
2066
2067 thr = zthr_1sid[ipthr] ;
2068 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2069 siz = find_largest_cluster_NN3( bfim , ithr ) ;
2070 if( siz > max_cluster_size ) siz = max_cluster_size ;
2071 mtab[siz]++ ; csiz_1sid_NN3[ipthr][iter-1] = siz ;
2072
2073 return ;
2074 }
2075
2076 /*---------------------------------------------------------------------------*/
2077 /* Find clusters, save some info */
2078
gather_stats_NN1_2sid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2079 void gather_stats_NN1_2sid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2080 {
2081 register int ii ; register float thr ; int siz ;
2082
2083 thr = zthr_2sid[ipthr] ;
2084 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) || (fim[ii] < -thr) ;
2085 siz = find_largest_cluster_NN1( bfim , ithr ) ;
2086 if( siz > max_cluster_size ) siz = max_cluster_size ;
2087 mtab[siz]++ ; csiz_2sid_NN1[ipthr][iter-1] = siz ;
2088
2089 return ;
2090 }
2091
2092 /*---------------------------------------------------------------------------*/
2093 /* Find clusters, save some info */
2094
gather_stats_NN2_2sid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2095 void gather_stats_NN2_2sid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2096 {
2097 register int ii ; register float thr ; int siz ;
2098
2099 thr = zthr_2sid[ipthr] ;
2100 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) || (fim[ii] < -thr) ;
2101 siz = find_largest_cluster_NN2( bfim , ithr ) ;
2102 if( siz > max_cluster_size ) siz = max_cluster_size ;
2103 mtab[siz]++ ; csiz_2sid_NN2[ipthr][iter-1] = siz ;
2104
2105 return ;
2106 }
2107
2108 /*---------------------------------------------------------------------------*/
2109 /* Find clusters, save some info */
2110
gather_stats_NN3_2sid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2111 void gather_stats_NN3_2sid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2112 {
2113 register int ii ; register float thr ; int siz ;
2114
2115 thr = zthr_2sid[ipthr] ;
2116 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) || (fim[ii] < -thr) ;
2117 siz = find_largest_cluster_NN3( bfim , ithr ) ;
2118 if( siz > max_cluster_size ) siz = max_cluster_size ;
2119 mtab[siz]++ ; csiz_2sid_NN3[ipthr][iter-1] = siz ;
2120
2121 return ;
2122 }
2123
2124 /*---------------------------------------------------------------------------*/
2125 /* Find clusters, save some info */
2126
gather_stats_NN1_bsid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2127 void gather_stats_NN1_bsid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2128 {
2129 register int ii ; register float thr ; int siz_p , siz_m , siz ;
2130
2131 thr = zthr_2sid[ipthr] ;
2132 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2133 siz_p = find_largest_cluster_NN1( bfim , ithr ) ;
2134 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] < -thr) ;
2135 siz_m = find_largest_cluster_NN1( bfim , ithr ) ;
2136 siz = MAX(siz_p,siz_m) ;
2137 if( siz > max_cluster_size ) siz = max_cluster_size ;
2138 mtab[siz]++ ; csiz_bsid_NN1[ipthr][iter-1] = siz ;
2139
2140 return ;
2141 }
2142
2143 /*---------------------------------------------------------------------------*/
2144 /* Find clusters, save some info */
2145
gather_stats_NN2_bsid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2146 void gather_stats_NN2_bsid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2147 {
2148 register int ii ; register float thr ; int siz_p , siz_m , siz ;
2149
2150 thr = zthr_2sid[ipthr] ;
2151 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2152 siz_p = find_largest_cluster_NN2( bfim , ithr ) ;
2153 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] < -thr) ;
2154 siz_m = find_largest_cluster_NN2( bfim , ithr ) ;
2155 siz = MAX(siz_p,siz_m) ;
2156 if( siz > max_cluster_size ) siz = max_cluster_size ;
2157 mtab[siz]++ ; csiz_bsid_NN2[ipthr][iter-1] = siz ;
2158
2159 return ;
2160 }
2161
2162 /*---------------------------------------------------------------------------*/
2163 /* Find clusters, save some info */
2164
gather_stats_NN3_bsid(int ipthr,float * fim,byte * bfim,int * mtab,int ithr,int iter)2165 void gather_stats_NN3_bsid( int ipthr, float *fim, byte *bfim, int *mtab, int ithr,int iter )
2166 {
2167 register int ii ; register float thr ; int siz_p , siz_m , siz ;
2168
2169 thr = zthr_2sid[ipthr] ;
2170 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2171 siz_p = find_largest_cluster_NN3( bfim , ithr ) ;
2172 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] < -thr) ;
2173 siz_m = find_largest_cluster_NN3( bfim , ithr ) ;
2174 siz = MAX(siz_p,siz_m) ;
2175 if( siz > max_cluster_size ) siz = max_cluster_size ;
2176 mtab[siz]++ ; csiz_bsid_NN3[ipthr][iter-1] = siz ;
2177
2178 return ;
2179 }
2180
2181 /*---------------------------------------------------------------------------*/
2182
prob6(float p)2183 static char * prob6(float p) /* format p-value into 6 char */
2184 {
2185 static char str[16] ;
2186 if( p >= 0.00010f ){
2187 sprintf(str,"%7.5f",p) ;
2188 } else {
2189 int dec = (int)(0.9999-log10(p)) ;
2190 float man = p * pow( 10.0 , (double)dec ) ;
2191 sprintf(str,"%4.1fe-%1d",man,dec) ;
2192 }
2193 return (str+1) ;
2194 }
2195
2196 /*---------------------------------------------------------------------------*/
2197
prob9(float p)2198 static char * prob9(float p) /* format p-value into 9 char */
2199 {
2200 static char str[16] ;
2201 if( p >= 0.00010f ){
2202 sprintf(str,"%9.6f",p) ;
2203 } else {
2204 int dec = (int)(0.9999-log10(p)) ;
2205 float man = p * pow( 10.0 , (double)dec ) ;
2206 sprintf(str,"%6.3fe-%1d",man,dec) ;
2207 }
2208 return str ;
2209 }
2210
2211 /*---------------------------------------------------------------------------*/
2212
2213 #ifdef USE_SHAVE
2214 static int *fa_1sid_NN1, *fa_1sid_NN2, *fa_1sid_NN3 ;
2215 static int *fa_2sid_NN1, *fa_2sid_NN2, *fa_2sid_NN3 ;
2216 static int *fa_bsid_NN1, *fa_bsid_NN2, *fa_bsid_NN3 ;
2217
thresh_summer_fim(int iathr,int ipthr_bot,int ipthr_top,float * fim,byte * bfim,int ithr)2218 void thresh_summer_fim( int iathr, int ipthr_bot, int ipthr_top,
2219 float *fim, byte *bfim, int ithr )
2220 {
2221 register int ii ; register float thr ;
2222 int ipthr , siz ;
2223
2224 fa_1sid_NN1[ithr]=0; fa_1sid_NN2[ithr]=0; fa_1sid_NN3[ithr]=0;
2225 fa_2sid_NN1[ithr]=0; fa_2sid_NN2[ithr]=0; fa_2sid_NN3[ithr]=0;
2226 fa_bsid_NN1[ithr]=0; fa_bsid_NN2[ithr]=0; fa_bsid_NN3[ithr]=0;
2227
2228 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2229
2230 if( !fa_1sid_NN1[ithr] ){ /* 1-sided, NN1 */
2231 thr = zthr_1sid[ipthr] ;
2232 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2233 siz = find_largest_cluster_NN1(bfim,0) ;
2234 fa_1sid_NN1[ithr] = (siz >= clust_thresh_1sid_NN1[ipthr][iathr] ) ;
2235 }
2236
2237 if( !fa_1sid_NN2[ithr] ){ /* 1-sided, NN2 */
2238 thr = zthr_1sid[ipthr] ;
2239 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2240 siz = find_largest_cluster_NN2(bfim,0) ;
2241 fa_1sid_NN2[ithr] = (siz >= clust_thresh_1sid_NN2[ipthr][iathr] ) ;
2242 }
2243
2244 if( !fa_1sid_NN3[ithr] ){ /* 1-sided, NN3 */
2245 thr = zthr_1sid[ipthr] ;
2246 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2247 siz = find_largest_cluster_NN3(bfim,0) ;
2248 fa_1sid_NN3[ithr] = (siz >= clust_thresh_1sid_NN3[ipthr][iathr] ) ;
2249 }
2250
2251 if( do_athr_sum < 2 ) continue ; /* skip 2-sided and bi-sided */
2252
2253 if( !fa_2sid_NN1[ithr] ){ /* 2-sided, NN1 */
2254 thr = zthr_2sid[ipthr] ;
2255 for( ii=0 ; ii < nxyz ; ii++ )
2256 bfim[ii] = (fim[ii] > thr) || (fim[ii] < -thr) ;
2257 siz = find_largest_cluster_NN1(bfim,0) ;
2258 fa_2sid_NN1[ithr] = (siz >= clust_thresh_2sid_NN1[ipthr][iathr] ) ;
2259 }
2260
2261 if( !fa_2sid_NN2[ithr] ){ /* 2-sided, NN2 */
2262 thr = zthr_2sid[ipthr] ;
2263 for( ii=0 ; ii < nxyz ; ii++ )
2264 bfim[ii] = (fim[ii] > thr) || (fim[ii] < -thr) ;
2265 siz = find_largest_cluster_NN2(bfim,0) ;
2266 fa_2sid_NN2[ithr] = (siz >= clust_thresh_2sid_NN2[ipthr][iathr] ) ;
2267 }
2268
2269 if( !fa_2sid_NN3[ithr] ){ /* 2-sided, NN3 */
2270 thr = zthr_2sid[ipthr] ;
2271 for( ii=0 ; ii < nxyz ; ii++ )
2272 bfim[ii] = (fim[ii] > thr) || (fim[ii] < -thr) ;
2273 siz = find_largest_cluster_NN3(bfim,0) ;
2274 fa_2sid_NN3[ithr] = (siz >= clust_thresh_2sid_NN3[ipthr][iathr] ) ;
2275 }
2276
2277 if( !fa_bsid_NN1[ithr] ){ /* bi-sided, NN1 */
2278 thr = zthr_2sid[ipthr] ;
2279 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2280 siz = find_largest_cluster_NN1( bfim , 0 ) ;
2281 fa_bsid_NN1[ithr] = (siz >= clust_thresh_bsid_NN1[ipthr][iathr] ) ;
2282 if( !fa_bsid_NN1[ithr] ){
2283 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] < -thr) ;
2284 siz = find_largest_cluster_NN1( bfim , 0 ) ;
2285 fa_bsid_NN1[ithr] = (siz >= clust_thresh_bsid_NN1[ipthr][iathr] ) ;
2286 }
2287 }
2288
2289 if( !fa_bsid_NN2[ithr] ){ /* bi-sided, NN2 */
2290 thr = zthr_2sid[ipthr] ;
2291 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2292 siz = find_largest_cluster_NN2( bfim , 0 ) ;
2293 fa_bsid_NN2[ithr] = (siz >= clust_thresh_bsid_NN2[ipthr][iathr] ) ;
2294 if( !fa_bsid_NN2[ithr] ){
2295 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] < -thr) ;
2296 siz = find_largest_cluster_NN2( bfim , 0 ) ;
2297 fa_bsid_NN2[ithr] = (siz >= clust_thresh_bsid_NN2[ipthr][iathr] ) ;
2298 }
2299 }
2300
2301 if( !fa_bsid_NN3[ithr] ){ /* bi-sided, NN3 */
2302 thr = zthr_2sid[ipthr] ;
2303 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] > thr) ;
2304 siz = find_largest_cluster_NN3( bfim , 0 ) ;
2305 fa_bsid_NN3[ithr] = (siz >= clust_thresh_bsid_NN3[ipthr][iathr] ) ;
2306 if( !fa_bsid_NN3[ithr] ){
2307 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] < -thr) ;
2308 siz = find_largest_cluster_NN3( bfim , 0 ) ;
2309 fa_bsid_NN3[ithr] = (siz >= clust_thresh_bsid_NN3[ipthr][iathr] ) ;
2310 }
2311 }
2312
2313 }
2314
2315 return ;
2316 }
2317
2318 /*---------------------------------------------------------------------------*/
2319
2320 static float *rfa_1sid_NN1, *rfa_1sid_NN2, *rfa_1sid_NN3 ;
2321 static float *rfa_2sid_NN1, *rfa_2sid_NN2, *rfa_2sid_NN3 ;
2322 static float *rfa_bsid_NN1, *rfa_bsid_NN2, *rfa_bsid_NN3 ;
2323
thresh_summer_athr(int ipthr_bot,int ipthr_top)2324 void thresh_summer_athr( int ipthr_bot , int ipthr_top )
2325 {
2326 float const drfa = 1.0f/niter ;
2327
2328 ENTRY("thresh_summer_athr") ;
2329
2330 rfa_1sid_NN1 = (float *)calloc(sizeof(float),nathr) ;
2331 rfa_2sid_NN1 = (float *)calloc(sizeof(float),nathr) ;
2332 rfa_bsid_NN1 = (float *)calloc(sizeof(float),nathr) ;
2333 rfa_1sid_NN2 = (float *)calloc(sizeof(float),nathr) ;
2334 rfa_2sid_NN2 = (float *)calloc(sizeof(float),nathr) ;
2335 rfa_bsid_NN2 = (float *)calloc(sizeof(float),nathr) ;
2336 rfa_1sid_NN3 = (float *)calloc(sizeof(float),nathr) ;
2337 rfa_2sid_NN3 = (float *)calloc(sizeof(float),nathr) ;
2338 rfa_bsid_NN3 = (float *)calloc(sizeof(float),nathr) ;
2339
2340
2341 AFNI_OMP_START;
2342 #pragma omp parallel
2343 { DECLARE_ithr ;
2344 int iter , iathr , vstep , vii ;
2345 float *fim ; byte *bfim ;
2346
2347 #pragma omp master
2348 {
2349 #ifdef USE_OMP
2350 nthr = omp_get_num_threads() ;
2351 #else
2352 nthr = 1 ;
2353 #endif
2354 fa_1sid_NN1 = (int *)malloc(sizeof(int)*nthr) ;
2355 fa_2sid_NN1 = (int *)malloc(sizeof(int)*nthr) ;
2356 fa_bsid_NN1 = (int *)malloc(sizeof(int)*nthr) ;
2357 fa_1sid_NN2 = (int *)malloc(sizeof(int)*nthr) ;
2358 fa_2sid_NN2 = (int *)malloc(sizeof(int)*nthr) ;
2359 fa_bsid_NN2 = (int *)malloc(sizeof(int)*nthr) ;
2360 fa_1sid_NN3 = (int *)malloc(sizeof(int)*nthr) ;
2361 fa_2sid_NN3 = (int *)malloc(sizeof(int)*nthr) ;
2362 fa_bsid_NN3 = (int *)malloc(sizeof(int)*nthr) ;
2363 if( verb ){
2364 vstep = (int)( niter / (nthr*50.0f) + 0.901f) ;
2365 vii = 0 ; vstep_reset() ;
2366 fprintf(stderr,"Summing up:") ;
2367 }
2368 }
2369 #pragma omp barrier
2370
2371 fim = (float *)malloc(sizeof(float)*nxyz) ;
2372 bfim = (byte *)malloc(sizeof(byte) *nxyz) ;
2373
2374 #pragma omp for
2375 for( iter=1 ; iter <= niter ; iter++ ){
2376
2377 if( ithr==0 && verb ){
2378 vii++ ; if( vii%vstep == vstep/2 ) vstep_print() ;
2379 }
2380
2381 shave_to_fim(fim,iter-1) ;
2382
2383 for( iathr=0 ; iathr < nathr ; iathr++ ){
2384
2385 thresh_summer_fim( iathr,ipthr_bot,fim,bfim,ithr ) ;
2386
2387 #pragma omp critical
2388 {if( fa_1sid_NN1[ithr] ) rfa_1sid_NN1[iathr] += drfa ;
2389 if( fa_1sid_NN2[ithr] ) rfa_1sid_NN2[iathr] += drfa ;
2390 if( fa_1sid_NN3[ithr] ) rfa_1sid_NN3[iathr] += drfa ;
2391 if( do_athr_sum > 1 ){
2392 if( fa_2sid_NN1[ithr] ) rfa_2sid_NN1[iathr] += drfa ;
2393 if( fa_2sid_NN2[ithr] ) rfa_2sid_NN2[iathr] += drfa ;
2394 if( fa_2sid_NN3[ithr] ) rfa_2sid_NN3[iathr] += drfa ;
2395 if( fa_bsid_NN1[ithr] ) rfa_bsid_NN1[iathr] += drfa ;
2396 if( fa_bsid_NN2[ithr] ) rfa_bsid_NN2[iathr] += drfa ;
2397 if( fa_bsid_NN3[ithr] ) rfa_bsid_NN3[iathr] += drfa ;
2398 }
2399 }
2400 }
2401 }
2402
2403 free(bfim) ; free(fim) ;
2404
2405 }
2406 AFNI_OMP_END ;
2407 if( verb ) fprintf(stderr,"\n") ;
2408
2409 EXRETURN ;
2410 }
2411
2412 #else
2413
2414 /*---------------------------------------------------------------------------*/
2415
2416 static float *rfa_1sid_NN1, *rfa_1sid_NN2, *rfa_1sid_NN3 ;
2417 static float *rfa_2sid_NN1, *rfa_2sid_NN2, *rfa_2sid_NN3 ;
2418 static float *rfa_bsid_NN1, *rfa_bsid_NN2, *rfa_bsid_NN3 ;
2419
thresh_summer_onecurve(int iathr,int ipthr_bot,int ipthr_top)2420 void thresh_summer_onecurve( int iathr , int ipthr_bot , int ipthr_top )
2421 {
2422 float const drfa = 1.0f/niter ;
2423 int ipthr, siz , iter ;
2424
2425 ENTRY("thresh_summer_onecurve") ;
2426
2427 for( iter=0 ; iter < niter ; iter++ ){
2428
2429 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2430 siz = csiz_1sid_NN1[ipthr][iter] ;
2431 if( siz >= clust_thresh_1sid_NN1[ipthr][iathr] ){
2432 rfa_1sid_NN1[iathr] += drfa ; break ;
2433 }
2434 }
2435
2436 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2437 siz = csiz_1sid_NN2[ipthr][iter] ;
2438 if( siz >= clust_thresh_1sid_NN2[ipthr][iathr] ){
2439 rfa_1sid_NN2[iathr] += drfa ; break ;
2440 }
2441 }
2442
2443 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2444 siz = csiz_1sid_NN3[ipthr][iter] ;
2445 if( siz >= clust_thresh_1sid_NN3[ipthr][iathr] ){
2446 rfa_1sid_NN3[iathr] += drfa ; break ;
2447 }
2448 }
2449
2450 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2451 siz = csiz_2sid_NN1[ipthr][iter] ;
2452 if( siz >= clust_thresh_2sid_NN1[ipthr][iathr] ){
2453 rfa_2sid_NN1[iathr] += drfa ; break ;
2454 }
2455 }
2456
2457 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2458 siz = csiz_2sid_NN2[ipthr][iter] ;
2459 if( siz >= clust_thresh_2sid_NN2[ipthr][iathr] ){
2460 rfa_2sid_NN2[iathr] += drfa ; break ;
2461 }
2462 }
2463
2464 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2465 siz = csiz_2sid_NN3[ipthr][iter] ;
2466 if( siz >= clust_thresh_2sid_NN3[ipthr][iathr] ){
2467 rfa_2sid_NN3[iathr] += drfa ; break ;
2468 }
2469 }
2470
2471 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2472 siz = csiz_bsid_NN1[ipthr][iter] ;
2473 if( siz >= clust_thresh_bsid_NN1[ipthr][iathr] ){
2474 rfa_bsid_NN1[iathr] += drfa ; break ;
2475 }
2476 }
2477
2478 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2479 siz = csiz_bsid_NN2[ipthr][iter] ;
2480 if( siz >= clust_thresh_bsid_NN2[ipthr][iathr] ){
2481 rfa_bsid_NN2[iathr] += drfa ; break ;
2482 }
2483 }
2484
2485 for( ipthr=ipthr_bot ; ipthr <= ipthr_top ; ipthr++ ){
2486 siz = csiz_bsid_NN3[ipthr][iter] ;
2487 if( siz >= clust_thresh_bsid_NN3[ipthr][iathr] ){
2488 rfa_bsid_NN3[iathr] += drfa ; break ;
2489 }
2490 }
2491
2492 }
2493
2494 EXRETURN ;
2495 }
2496
2497 /*---------------------------------------------------------------------------*/
2498
thresh_summer_athr(int ipthr_bot,int ipthr_top)2499 void thresh_summer_athr( int ipthr_bot , int ipthr_top )
2500 {
2501 int iathr ;
2502
2503 if( ipthr_top < ipthr_bot ) ipthr_top = npthr-1 ;
2504
2505 ENTRY("thresh_summer_athr") ;
2506
2507 rfa_1sid_NN1 = (float *)calloc(sizeof(float),nathr) ;
2508 rfa_2sid_NN1 = (float *)calloc(sizeof(float),nathr) ;
2509 rfa_bsid_NN1 = (float *)calloc(sizeof(float),nathr) ;
2510 rfa_1sid_NN2 = (float *)calloc(sizeof(float),nathr) ;
2511 rfa_2sid_NN2 = (float *)calloc(sizeof(float),nathr) ;
2512 rfa_bsid_NN2 = (float *)calloc(sizeof(float),nathr) ;
2513 rfa_1sid_NN3 = (float *)calloc(sizeof(float),nathr) ;
2514 rfa_2sid_NN3 = (float *)calloc(sizeof(float),nathr) ;
2515 rfa_bsid_NN3 = (float *)calloc(sizeof(float),nathr) ;
2516
2517 for( iathr=0 ; iathr < nathr ; iathr++ ){
2518 thresh_summer_onecurve(iathr,ipthr_bot,ipthr_top) ;
2519 }
2520
2521 EXRETURN ;
2522 }
2523 #endif
2524
2525 /*---------------------------------------------------------------------------*/
2526 /* from the max_table,
2527 get the threshold cluster size for a given per-voxel threshold,
2528 for a given false alarm rate (alpha).
2529 *//*-------------------------------------------------------------------------*/
2530
get_one_clust_thresh(int ** mtnn,int ipthr,double aval)2531 double get_one_clust_thresh( int **mtnn , int ipthr , double aval )
2532 {
2533 static double *alpha=NULL ;
2534 double ahi,alo,jj ;
2535 float cmax=0.0f , **clust_thresh ;
2536 int ii , itop , iathr , mmm ;
2537
2538 if( alpha == NULL )
2539 alpha = (double *)malloc(sizeof(double)*(max_cluster_size+1)) ;
2540
2541 /* alpha[ii] = prob of getting clusters of exactly size ii */
2542
2543 for( itop=ii=1 ; ii <= max_cluster_size ; ii++ ){
2544 alpha[ii] = mtnn[ipthr][ii] / (double)niter ;
2545 if( alpha[ii] > 0.0 ) itop = ii ;
2546 }
2547
2548 /* alpha[ii] = prob of getting clusters of size ii or more */
2549
2550 for( ii=itop-1 ; ii >= 1 ; ii-- ) alpha[ii] += alpha[ii+1] ;
2551
2552 /* find ii that brackets the desired aval */
2553
2554 if( aval > alpha[1] ){ /* unpleasant situation */
2555 static int first=2 ; /* where there is no bracket */
2556 ii = 1 ;
2557 if( first ){
2558 WARNING_message(
2559 "pthr=%.6f ; desired alpha=%.6f -- but max simulated alpha=%.6f]\n" ,
2560 pthr[ipthr] , aval , alpha[1] ) ;
2561 first-- ;
2562 if( first == 0 )
2563 ININFO_message(" [further messages about alpha are suppressed]") ;
2564 }
2565 } else {
2566 for( ii=1 ; ii < itop ; ii++ ){
2567 if( alpha[ii] >= aval && alpha[ii+1] <= aval ) break ;
2568 }
2569 }
2570
2571 /* inverse interpolate to find the index jj where alpha[jj]=aval
2572 -- except, of course, jj is a floating point value between ii and ii+1 */
2573
2574 alo=alpha[ii] ; ahi=alpha[ii+1] ;
2575
2576 if( alo >= 1.0 ) alo = 1.0 - 0.1/niter ; /* unlikely */
2577 if( ahi <= 0.0 ) ahi = 0.1/niter ; /* should not be possible */
2578 if( ahi >= alo ) ahi = 0.1*alo ; /* should not be possible */
2579 jj = log(-log(1.0-aval)) ;
2580 alo = log(-log(1.0-alo)) ;
2581 ahi = log(-log(1.0-ahi)) ;
2582 jj = ii + (jj-alo)/(ahi-alo) ;
2583 if( jj < 1.0 ) jj = 1.0 ;
2584
2585 return jj ;
2586 }
2587
2588 /*===========================================================================*/
2589
main(int argc,char ** argv)2590 int main( int argc , char **argv )
2591 {
2592 int nnn , ipthr , first_mask=1 ;
2593 char *refit_cmd = NULL ;
2594 double ct ;
2595 #ifdef USE_OMP
2596 int ***mtab_1sid[4] , ***mtab_2sid[4] , ***mtab_bsid[4] ;
2597 #endif
2598
2599 /*----- does user request help menu? -----*/
2600
2601 AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */
2602 (void)COX_clock_time() ;
2603
2604 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) display_help_menu() ;
2605
2606 /*----- get the list of things to do -----*/
2607
2608 mainENTRY("3dClustSim"); machdep();
2609 AFNI_logger("3dClustSim",argc,argv);
2610 PRINT_VERSION("3dClustSim"); AUTHOR("RW Cox and BD Ward");
2611 THD_check_AFNI_version("3dClustSim") ;
2612
2613 get_options( argc , argv ) ;
2614
2615 /*----- create some space for the results -----*/
2616
2617 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2618 max_table_1sid[nnn] = (int **)malloc(sizeof(int *)*npthr) ; /* array of tables */
2619 max_table_2sid[nnn] = (int **)malloc(sizeof(int *)*npthr) ; /* array of tables */
2620 max_table_bsid[nnn] = (int **)malloc(sizeof(int *)*npthr) ; /* array of tables */
2621 cmx_table_1sid[nnn] = (int * )malloc(sizeof(int )*npthr) ; /* 23 Dec 2015 */
2622 cmx_table_2sid[nnn] = (int * )malloc(sizeof(int )*npthr) ;
2623 cmx_table_bsid[nnn] = (int * )malloc(sizeof(int )*npthr) ;
2624 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){ /* create tables */
2625 max_table_1sid[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
2626 max_table_2sid[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
2627 max_table_bsid[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
2628 }
2629 }
2630
2631 csiz_1sid_NN1 = (int **)malloc(sizeof(int *)*npthr) ; /* 23 Dec 2015 */
2632 csiz_1sid_NN2 = (int **)malloc(sizeof(int *)*npthr) ;
2633 csiz_1sid_NN3 = (int **)malloc(sizeof(int *)*npthr) ;
2634 csiz_2sid_NN1 = (int **)malloc(sizeof(int *)*npthr) ;
2635 csiz_2sid_NN2 = (int **)malloc(sizeof(int *)*npthr) ;
2636 csiz_2sid_NN3 = (int **)malloc(sizeof(int *)*npthr) ;
2637 csiz_bsid_NN1 = (int **)malloc(sizeof(int *)*npthr) ;
2638 csiz_bsid_NN2 = (int **)malloc(sizeof(int *)*npthr) ;
2639 csiz_bsid_NN3 = (int **)malloc(sizeof(int *)*npthr) ;
2640 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2641 csiz_1sid_NN1[ipthr] = (int *)calloc(sizeof(int),niter) ;
2642 csiz_1sid_NN2[ipthr] = (int *)calloc(sizeof(int),niter) ;
2643 csiz_1sid_NN3[ipthr] = (int *)calloc(sizeof(int),niter) ;
2644 csiz_2sid_NN1[ipthr] = (int *)calloc(sizeof(int),niter) ;
2645 csiz_2sid_NN2[ipthr] = (int *)calloc(sizeof(int),niter) ;
2646 csiz_2sid_NN3[ipthr] = (int *)calloc(sizeof(int),niter) ;
2647 csiz_bsid_NN1[ipthr] = (int *)calloc(sizeof(int),niter) ;
2648 csiz_bsid_NN2[ipthr] = (int *)calloc(sizeof(int),niter) ;
2649 csiz_bsid_NN3[ipthr] = (int *)calloc(sizeof(int),niter) ;
2650 }
2651
2652 #ifdef USE_SHAVE
2653 if( do_athr_sum ) setup_shave() ;
2654 #endif
2655
2656 if( verb )
2657 INFO_message("Startup clock time = %.1f s",COX_clock_time()) ;
2658
2659 AFNI_OMP_START ;
2660 #pragma omp parallel
2661 {
2662 DECLARE_ithr ;
2663 int iter, ipthr, **mt_1sid[4],**mt_2sid[4],**mt_bsid[4] , nnn ;
2664 float *fim ; byte *bfim ; unsigned short xran[3] ;
2665 float *pfim ;
2666 int vstep , vii ;
2667
2668 /* create separate tables for each thread, if using OpenMP */
2669 #ifdef USE_OMP
2670 #pragma omp master /* only in the master thread */
2671 {
2672 nthr = omp_get_num_threads() ;
2673 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2674 mtab_1sid[nnn] = (int ***)malloc(sizeof(int **) *nthr) ;
2675 mtab_2sid[nnn] = (int ***)malloc(sizeof(int **) *nthr) ;
2676 mtab_bsid[nnn] = (int ***)malloc(sizeof(int **) *nthr) ;
2677 }
2678
2679 nall_g = (int *) malloc(sizeof(int) *nthr) ; /* workspaces for */
2680 inow_g = (short **)malloc(sizeof(short *)*nthr) ; /* find_largest_cluster() */
2681 jnow_g = (short **)malloc(sizeof(short *)*nthr) ;
2682 know_g = (short **)malloc(sizeof(short *)*nthr) ;
2683 if( do_acf ){
2684 acf_aim = (MRI_IMAGE **)calloc(sizeof(MRI_IMAGE *),nthr) ;
2685 acf_bim = (MRI_IMAGE **)calloc(sizeof(MRI_IMAGE *),nthr) ;
2686 acf_tar = (complex **) malloc(sizeof(complex *) *nthr) ;
2687 }
2688 if( verb ) INFO_message("Using %d OpenMP threads",nthr) ;
2689 }
2690 #pragma omp barrier /* all threads wait until the above is finished */
2691 /* create tables for each thread separately */
2692 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2693 mtab_1sid[nnn][ithr] = mt_1sid[nnn] = (int **)malloc(sizeof(int *)*npthr) ;
2694 mtab_2sid[nnn][ithr] = mt_2sid[nnn] = (int **)malloc(sizeof(int *)*npthr) ;
2695 mtab_bsid[nnn][ithr] = mt_bsid[nnn] = (int **)malloc(sizeof(int *)*npthr) ;
2696 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2697 mt_1sid[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
2698 mt_2sid[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
2699 mt_bsid[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
2700 }
2701 }
2702
2703 /* create workspace for find_largest_cluster(), for this thread */
2704
2705 nall_g[ithr] = DALL ;
2706 inow_g[ithr] = (short *) malloc(sizeof(short)*DALL) ;
2707 jnow_g[ithr] = (short *) malloc(sizeof(short)*DALL) ;
2708 know_g[ithr] = (short *) malloc(sizeof(short)*DALL) ;
2709
2710 if( do_acf ){
2711 acf_tar[ithr] = (complex *)malloc(sizeof(complex)*acf_ntar) ;
2712 }
2713
2714 /* initialize random seed array for each thread separately */
2715 xran[2] = ( gseed & 0xffff) + (unsigned short)ithr ;
2716 xran[1] = ((gseed >> 16) & 0xffff) - (unsigned short)ithr ;
2717 xran[0] = 0x330e + (unsigned short)ithr ;
2718
2719 #else /* not OpenMP ==> only one set of tables */
2720 xran[2] = ( gseed & 0xffff) ;
2721 xran[1] = ((gseed >> 16) & 0xffff) ;
2722 xran[0] = 0x330e ;
2723 nall_g = (int *) malloc(sizeof(int) *nthr) ;
2724 inow_g = (short **)malloc(sizeof(short *)*nthr) ;
2725 jnow_g = (short **)malloc(sizeof(short *)*nthr) ;
2726 know_g = (short **)malloc(sizeof(short *)*nthr) ;
2727 nall_g[ithr] = DALL ;
2728 inow_g[ithr] = (short *) malloc(sizeof(short)*DALL) ;
2729 jnow_g[ithr] = (short *) malloc(sizeof(short)*DALL) ;
2730 know_g[ithr] = (short *) malloc(sizeof(short)*DALL) ;
2731 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2732 mt_1sid[nnn] = max_table_1sid[nnn] ;
2733 mt_2sid[nnn] = max_table_2sid[nnn] ;
2734 mt_bsid[nnn] = max_table_bsid[nnn] ;
2735 }
2736 if( do_acf ){
2737 acf_aim = (MRI_IMAGE **)calloc(sizeof(MRI_IMAGE *),nthr) ;
2738 acf_bim = (MRI_IMAGE **)calloc(sizeof(MRI_IMAGE *),nthr) ;
2739 acf_tar = (complex **) calloc(sizeof(complex *) ,nthr) ;
2740 acf_tar[ithr] = (complex *)malloc(sizeof(complex)*acf_ntar) ;
2741 }
2742 #endif
2743
2744 fim = (float *)malloc(sizeof(float)*nxyz) ; /* image space */
2745 bfim = (byte * )malloc(sizeof(byte) *nxyz) ;
2746 if( do_pad ) pfim = (float *)malloc(sizeof(float)*nxyz_pad); /* 12 May 2015 */
2747 else pfim = NULL ;
2748
2749 vstep = (int)( niter / (nthr*50.0f) + 0.901f) ;
2750 vii = 0 ;
2751 if( ithr == 0 && verb ) fprintf(stderr,"Simulating:") ;
2752
2753 /*----- Monte Carlo iterations -----*/
2754
2755 #pragma omp for
2756 for( iter=1 ; iter <= niter ; iter++ ){
2757
2758 if( verb && ithr == 0 ){
2759 vii++ ; if( vii%vstep == 2 ) vstep_print() ;
2760 }
2761
2762 generate_image( fim , pfim , xran , iter ) ;
2763
2764 #ifdef USE_SHAVE
2765 if( do_shave ) fim_to_shave(fim,iter-1) ;
2766 #endif
2767
2768 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2769 gather_stats_NN1_1sid( ipthr, fim, bfim, mt_1sid[1][ipthr], ithr,iter ) ;
2770 gather_stats_NN2_1sid( ipthr, fim, bfim, mt_1sid[2][ipthr], ithr,iter ) ;
2771 gather_stats_NN3_1sid( ipthr, fim, bfim, mt_1sid[3][ipthr], ithr,iter ) ;
2772
2773 gather_stats_NN1_2sid( ipthr, fim, bfim, mt_2sid[1][ipthr], ithr,iter ) ;
2774 gather_stats_NN2_2sid( ipthr, fim, bfim, mt_2sid[2][ipthr], ithr,iter ) ;
2775 gather_stats_NN3_2sid( ipthr, fim, bfim, mt_2sid[3][ipthr], ithr,iter ) ;
2776
2777 gather_stats_NN1_bsid( ipthr, fim, bfim, mt_bsid[1][ipthr], ithr,iter ) ;
2778 gather_stats_NN2_bsid( ipthr, fim, bfim, mt_bsid[2][ipthr], ithr,iter ) ;
2779 gather_stats_NN3_bsid( ipthr, fim, bfim, mt_bsid[3][ipthr], ithr,iter ) ;
2780 }
2781
2782 } /* end of simulation loop */
2783
2784 free(fim) ; free(bfim) ; if( pfim != NULL ) free(pfim) ;
2785
2786 if( !do_athr_sum ){
2787 free(inow_g[ithr]) ; free(jnow_g[ithr]) ; free(know_g[ithr]) ;
2788 if( do_acf ){
2789 free(acf_tar[ithr]) ;
2790 if( acf_bim[ithr] != NULL ) mri_free(acf_bim[ithr]) ;
2791 }
2792 }
2793
2794 if( ithr == 0 && verb ) fprintf(stderr,"!\n") ;
2795
2796 } /* end OpenMP parallelization */
2797 AFNI_OMP_END ;
2798
2799 if( verb )
2800 INFO_message("Clock time now = %.1f s",COX_clock_time()) ;
2801
2802 if( do_acf && !do_athr_sum ){
2803 free(acf_aim) ; free(acf_bim) ; free(acf_tar) ;
2804 }
2805
2806 if( do_ssave > 0 ) DSET_delete(ssave_dset) ; /* 24 Apr 2014 */
2807
2808 /*-------- sum tables from various threads into one result ----------*/
2809
2810 #ifdef USE_OMP
2811 { int ithr , ii , ipthr , **mt_1sid,**mt_2sid,**mt_bsid , *mth_1sid,*mth_2sid,*mth_bsid ;
2812 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2813 for( ithr=0 ; ithr < nthr ; ithr++ ){
2814 mt_1sid = mtab_1sid[nnn][ithr] ;
2815 mt_2sid = mtab_2sid[nnn][ithr] ;
2816 mt_bsid = mtab_bsid[nnn][ithr] ;
2817 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2818 mth_1sid = mt_1sid[ipthr] ;
2819 mth_2sid = mt_2sid[ipthr] ;
2820 mth_bsid = mt_bsid[ipthr] ;
2821 for( ii=1 ; ii <= max_cluster_size ; ii++ ){
2822 max_table_1sid[nnn][ipthr][ii] += mth_1sid[ii] ;
2823 max_table_2sid[nnn][ipthr][ii] += mth_2sid[ii] ;
2824 max_table_bsid[nnn][ipthr][ii] += mth_bsid[ii] ;
2825 }
2826 }
2827 }
2828 }
2829 }
2830 #endif
2831
2832 /* For each table, find the largest cluster size
2833 cc such that max_table_xxxx[nnn][ipthr][cc] > 0 [23 Dec 2015] */
2834
2835 { int cc ;
2836 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2837 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2838 for( cc=max_cluster_size ;
2839 cc >= 1 && max_table_1sid[nnn][ipthr][cc]==0 ; cc-- ) ; /*nada*/
2840 cmx_table_1sid[nnn][ipthr] = cc ;
2841 for( cc=max_cluster_size ;
2842 cc >= 1 && max_table_2sid[nnn][ipthr][cc]==0 ; cc-- ) ; /*nada*/
2843 cmx_table_2sid[nnn][ipthr] = cc ;
2844 for( cc=max_cluster_size ;
2845 cc >= 1 && max_table_bsid[nnn][ipthr][cc]==0 ; cc-- ) ; /*nada*/
2846 cmx_table_bsid[nnn][ipthr] = cc ;
2847 }
2848 }
2849 }
2850
2851 #if 1
2852 if( AFNI_yesenv("AFNI_CLUSTSIM_SAVE") ){
2853 FILE *fp; char fname[128]; static char *lnn[3] = { "NN1","NN2","NN3" }; int cc;
2854 for( nnn=1 ; nnn <= 3 ; nnn++ ){
2855 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2856 sprintf(fname,"mt.1sid.%s.%.5f.1D",lnn[nnn-1],pthr[ipthr]) ;
2857 fp = fopen(fname,"w") ;
2858 for( cc=1 ; cc <= cmx_table_1sid[nnn][ipthr] ; cc++ )
2859 fprintf(fp,"%d %d\n",cc,max_table_1sid[nnn][ipthr][cc]) ;
2860 fclose(fp) ;
2861 sprintf(fname,"mt.2sid.%s.%.5f.1D",lnn[nnn-1],pthr[ipthr]) ;
2862 fp = fopen(fname,"w") ;
2863 for( cc=1 ; cc <= cmx_table_2sid[nnn][ipthr] ; cc++ )
2864 fprintf(fp,"%d %d\n",cc,max_table_2sid[nnn][ipthr][cc]) ;
2865 fclose(fp) ;
2866 sprintf(fname,"mt.bsid.%s.%.5f.1D",lnn[nnn-1],pthr[ipthr]) ;
2867 fp = fopen(fname,"w") ;
2868 for( cc=1 ; cc <= cmx_table_bsid[nnn][ipthr] ; cc++ )
2869 fprintf(fp,"%d %d\n",cc,max_table_bsid[nnn][ipthr][cc]) ;
2870 fclose(fp) ;
2871 }
2872 }
2873 }
2874 #endif
2875
2876 #if 0
2877 enable_mcw_malloc() ;
2878 #endif
2879
2880 /*---------- compute and print the output tables ----------*/
2881
2882 { double *alpha , aval , ahi,alo ;
2883 float cmax=0.0f , **clust_thresh ;
2884 int ii , itop , iathr , mmm, ***mtt , **mtnn ;
2885 char *commandline = tross_commandline("3dClustSim",argc,argv) ;
2886 char fname[THD_MAX_NAME] , pname[THD_MAX_NAME] ;
2887 char *amesg = NULL , *mlab = NULL , *mlll = NULL ;
2888
2889 alpha = (double *)malloc(sizeof(double)*(max_cluster_size+1)) ;
2890
2891 for( mmm=1 ; mmm <= 3 ; mmm++ ){ /* loop over sidedness */
2892
2893 if( mmm == 1 ){
2894 mtt = max_table_1sid ; mlab = "1-sided" ; mlll = "1sided" ;
2895 } else if( mmm == 2 ){
2896 mtt = max_table_2sid ; mlab = "2-sided" ; mlll = "2sided" ;
2897 } else {
2898 mtt = max_table_bsid ; mlab = "bi-sided"; mlll = "bisided";
2899 }
2900
2901 amesg = NULL ; cmax = 0.0f ;
2902 for( nnn=1 ; nnn <= 3 ; nnn++ ){ /* loop over NN level */
2903
2904 clust_thresh = (float **)malloc(sizeof(float *)*npthr) ;
2905 for( ipthr=0 ; ipthr < npthr ; ipthr++ )
2906 clust_thresh[ipthr] = (float *)malloc(sizeof(float)*nathr) ;
2907
2908 mtnn = mtt[nnn] ;
2909
2910 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2911 #if 1
2912 for( iathr=0 ; iathr < nathr ; iathr++ ){
2913 clust_thresh[ipthr][iathr] =
2914 get_one_clust_thresh( mtnn,ipthr,athr[iathr]) ;
2915 if( nodec ) clust_thresh[ipthr][iathr] =
2916 (int)(clust_thresh[ipthr][iathr]+0.951f) ;
2917 if( clust_thresh[ipthr][iathr] > cmax ) cmax = clust_thresh[ipthr][iathr] ;
2918 }
2919 #else
2920 for( itop=ii=1 ; ii <= max_cluster_size ; ii++ ){
2921 alpha[ii] = mtt[nnn][ipthr][ii] / (double)niter ;
2922 if( alpha[ii] > 0.0 ) itop = ii ;
2923 }
2924 for( ii=itop-1 ; ii >= 1 ; ii-- ) alpha[ii] += alpha[ii+1] ;
2925 for( iathr=0 ; iathr < nathr ; iathr++ ){
2926 aval = athr[iathr] ;
2927 if( aval > alpha[1] ){ /* unpleasant situation */
2928 ii = 1 ;
2929 amesg = THD_zzprintf(
2930 amesg ,
2931 " %s NN=%d pthr=%9.6f alpha=%6.3f [max simulated alpha=%6.3f]\n" ,
2932 mlab , nnn , pthr[ipthr] , aval , alpha[1] ) ;
2933 } else {
2934 for( ii=1 ; ii < itop ; ii++ ){
2935 if( alpha[ii] >= aval && alpha[ii+1] <= aval ) break ;
2936 }
2937 }
2938
2939 alo=alpha[ii] ; ahi=alpha[ii+1] ;
2940 if( do_lohi ){
2941 aval = ii ; /* for debugging */
2942 } else {
2943 if( alo >= 1.0 ) alo = 1.0 - 0.1/niter ;
2944 if( ahi <= 0.0 ) ahi = 0.1/niter ;
2945 if( ahi >= alo ) ahi = 0.1*alo ;
2946 aval = log(-log(1.0-aval)) ;
2947 alo = log(-log(1.0-alo)) ;
2948 ahi = log(-log(1.0-ahi)) ;
2949 aval = ii + (aval-alo)/(ahi-alo) ;
2950 if( aval < 1.0 ) aval = 1.0 ;
2951 else if( nodec ) aval = (int)(aval+0.951) ;
2952 }
2953 clust_thresh[ipthr][iathr] = aval ;
2954
2955 if( clust_thresh[ipthr][iathr] > cmax ) cmax = clust_thresh[ipthr][iathr] ;
2956 }
2957 #endif
2958 } /* end of loop over ipthr */
2959
2960 if( do_lohi == 0 ){
2961 /* edit each column to increase as pthr increases [shouldn't be needed] */
2962
2963 for( iathr=0 ; iathr < nathr ; iathr++ ){
2964 for( ipthr=npthr-2 ; ipthr >= 0 ; ipthr-- ){
2965 if( clust_thresh[ipthr][iathr] < clust_thresh[ipthr+1][iathr] )
2966 clust_thresh[ipthr][iathr] = clust_thresh[ipthr+1][iathr] ;
2967 }
2968 }
2969
2970 /* edit each row to increase as athr decreases [shouldn't be needed] */
2971
2972 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
2973 for( iathr=1 ; iathr < nathr ; iathr++ ){
2974 if( clust_thresh[ipthr][iathr] < clust_thresh[ipthr][iathr-1] )
2975 clust_thresh[ipthr][iathr] = clust_thresh[ipthr][iathr-1] ;
2976 }
2977 }
2978 }
2979
2980 #if 0
2981 if( !nodec && !do_niml && cmax > 9999.9f ){ /* if largest is way big, */
2982 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){ /* then truncate to ints. */
2983 for( iathr=0 ; iathr < nathr ; iathr++ ){
2984 aval = clust_thresh[ipthr][iathr] ;
2985 aval = (int)(aval+0.951) ;
2986 clust_thresh[ipthr][iathr] = aval ;
2987 }
2988 }
2989 nodec = 1 ;
2990 }
2991 #endif
2992
2993 MPROBE ;
2994
2995 if( prefix != NULL ){
2996 sprintf(pname,"%s.NN%d_%s.",prefix,nnn,mlll) ;
2997 } else {
2998 fflush(stderr) ; fflush(stdout) ;
2999 }
3000
3001 if( do_1D ){ /* output in 1D format */
3002 FILE *fp = stdout ;
3003 if( prefix != NULL ){
3004 strcpy(fname,pname) ; strcat(fname,"1D") ; fp = fopen(fname,"w") ;
3005 if( fp == NULL ){
3006 ERROR_message("Can't open file %s -- using stdout",fname) ;
3007 fp = stdout ;
3008 }
3009 }
3010 fprintf(fp,
3011 "# %s\n"
3012 "# %s thresholding\n"
3013 "# Grid: %dx%dx%d %.2fx%.2fx%.2f mm^3 (%d voxels%s)\n"
3014 "#\n"
3015 "# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels\n"
3016 "# -NN %d | alpha = Prob(Cluster >= given size)\n"
3017 "# pthr |" ,
3018 commandline , mlab ,
3019 nx,ny,nz , dx,dy,dz ,
3020 mask_ngood , (mask_ngood < nxyz) ? " in mask" : "\0" , nnn ) ;
3021 for( iathr=0 ; iathr < nathr ; iathr++ ) fprintf(fp," %s",prob6(athr[iathr])) ;
3022 fprintf(fp,"\n"
3023 "# ------ |" ) ;
3024 for( iathr=0 ; iathr < nathr ; iathr++ ) fprintf(fp," ------") ;
3025 fprintf(fp,"\n") ;
3026 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
3027 fprintf(fp,"%s ",prob9(pthr[ipthr])) ;
3028 for( iathr=0 ; iathr < nathr ; iathr++ ){
3029 if( nodec )
3030 fprintf(fp,"%7d" ,(int)clust_thresh[ipthr][iathr]) ;
3031 else if( clust_thresh[ipthr][iathr] <= 9999.9f )
3032 fprintf(fp,"%7.1f", clust_thresh[ipthr][iathr]) ;
3033 else
3034 fprintf(fp,"%7.0f", clust_thresh[ipthr][iathr]) ;
3035 }
3036 fprintf(fp,"\n") ;
3037 }
3038 if( fp != stdout ) fclose(fp) ;
3039 }
3040
3041 if( do_niml ){ /* output in NIML format */
3042 NI_element *nel ; float *vec ; char buf[1024] , *bbb ; NI_float_array nfar ;
3043 sprintf(buf,"3dClustSim_NN%d",nnn) ;
3044 nel = NI_new_data_element( buf , npthr ) ;
3045 vec = (float *)malloc(sizeof(float)*MAX(npthr,nathr)) ;
3046 for( iathr=0 ; iathr < nathr ; iathr++ ){
3047 for( ipthr=0 ; ipthr < npthr ; ipthr++ )
3048 vec[ipthr] = clust_thresh[ipthr][iathr] ;
3049 NI_add_column( nel , NI_FLOAT , vec ) ;
3050 }
3051 NI_set_attribute( nel , "commandline" , commandline ) ;
3052 NI_set_attribute( nel , "thresholding" , mlab ) ;
3053 sprintf(buf,"%d,%d,%d",nx,ny,nz) ;
3054 NI_set_attribute(nel,"nxyz",buf) ;
3055 sprintf(buf,"%.3f,%.3f,%.3f",dx,dy,dz) ;
3056 NI_set_attribute(nel,"dxyz",buf) ;
3057 sprintf(buf,"%.2f,%.2f,%.2f",fwhm_x,fwhm_y,fwhm_z) ;
3058 NI_set_attribute(nel,"fwhmxyz",buf) ;
3059 sprintf(buf,"%d",niter) ;
3060 NI_set_attribute(nel,"iter",buf) ;
3061 for( ipthr=0 ; ipthr < npthr ; ipthr++ ) vec[ipthr] = pthr[ipthr] ;
3062 nfar.num = npthr ; nfar.ar = vec ; bbb = NI_encode_float_list(&nfar,",") ;
3063 NI_set_attribute(nel,"pthr",bbb) ; NI_free(bbb) ;
3064 for( iathr=0 ; iathr < nathr ; iathr++ ) vec[iathr] = athr[iathr] ;
3065 nfar.num = nathr ; nfar.ar = vec ; bbb = NI_encode_float_list(&nfar,",") ;
3066 NI_set_attribute(nel,"athr",bbb) ; NI_free(bbb) ;
3067 if( mask_dset != NULL ){
3068 NI_set_attribute(nel,"mask_dset_idcode",DSET_IDCODE_STR(mask_dset)) ;
3069 NI_set_attribute(nel,"mask_dset_name" ,DSET_HEADNAME(mask_dset)) ;
3070 sprintf(buf,"%d",mask_ngood) ;
3071 NI_set_attribute(nel,"mask_count",buf) ;
3072 }
3073 if( prefix != NULL ){ strcpy(fname,pname) ; strcat(fname,"niml") ; }
3074 else strcpy(fname,"stdout:") ;
3075 NI_write_element_tofile( fname , nel , NI_TEXT_MODE ) ;
3076 NI_free_element( nel ) ;
3077 if( prefix != NULL ){
3078 if( refit_cmd == NULL )
3079 refit_cmd = THD_zzprintf( refit_cmd , "3drefit " ) ;
3080 refit_cmd = THD_zzprintf( refit_cmd ,
3081 "-atrstring AFNI_CLUSTSIM_NN%d_%s file:%s " ,
3082 nnn , mlll , fname ) ;
3083 }
3084 if( prefix != NULL && mask_vol != NULL && first_mask ){
3085 bbb = mask_to_b64string( mask_nvox , mask_vol ) ; first_mask = 0 ;
3086 if( bbb != NULL ){
3087 FILE *fp ;
3088 sprintf(fname,"%s.mask",prefix) ; fp = fopen(fname,"w") ;
3089 if( fp != NULL ){ fprintf(fp,"%s\n",bbb); fclose(fp); }
3090 free(bbb) ;
3091 refit_cmd = THD_zzprintf( refit_cmd ,
3092 "-atrstring AFNI_CLUSTSIM_MASK file:%s " ,
3093 fname) ;
3094 }
3095 }
3096 } /* end of NIML output */
3097
3098 #if 0
3099 switch( nnn*10 + mmm ){
3100 case 11: clust_thresh_1sid_NN1 = clust_thresh ; break ;
3101 case 12: clust_thresh_2sid_NN1 = clust_thresh ; break ;
3102 case 13: clust_thresh_bsid_NN1 = clust_thresh ; break ;
3103 case 21: clust_thresh_1sid_NN2 = clust_thresh ; break ;
3104 case 22: clust_thresh_2sid_NN2 = clust_thresh ; break ;
3105 case 23: clust_thresh_bsid_NN2 = clust_thresh ; break ;
3106 case 31: clust_thresh_1sid_NN3 = clust_thresh ; break ;
3107 case 32: clust_thresh_2sid_NN3 = clust_thresh ; break ;
3108 case 33: clust_thresh_bsid_NN3 = clust_thresh ; break ;
3109 default: /* should never transpire */
3110 ERROR_exit("nnn=%d mmm=%d :: This should never happen!",nnn,mmm) ;
3111 break ;
3112 }
3113 #else
3114 if( nnn*10+mmm == 11 ) clust_thresh_1sid_NN1 = clust_thresh ;
3115 if( nnn*10+mmm == 12 ) clust_thresh_2sid_NN1 = clust_thresh ;
3116 if( nnn*10+mmm == 13 ) clust_thresh_bsid_NN1 = clust_thresh ;
3117 if( nnn*10+mmm == 21 ) clust_thresh_1sid_NN2 = clust_thresh ;
3118 if( nnn*10+mmm == 22 ) clust_thresh_2sid_NN2 = clust_thresh ;
3119 if( nnn*10+mmm == 23 ) clust_thresh_bsid_NN2 = clust_thresh ;
3120 if( nnn*10+mmm == 31 ) clust_thresh_1sid_NN3 = clust_thresh ;
3121 if( nnn*10+mmm == 32 ) clust_thresh_2sid_NN3 = clust_thresh ;
3122 if( nnn*10+mmm == 33 ) clust_thresh_bsid_NN3 = clust_thresh ;
3123 #endif
3124
3125 } /* end of loop over nnn = NN degree */
3126
3127 if( amesg != NULL ){
3128 /* WARNING_message cannot handle "%s" with a very long string (>16K),
3129 so break this up a little 20 Feb 2014 [rickr] */
3130 /* option: change to vsnprintf() in output_message()? */
3131 WARNING_message("Simulation not effective for these cases:\n\n");
3132 fprintf(stderr, "%s\n", amesg);
3133 fprintf(stderr,
3134 "*+ This means that not enough clusters, of any size, +*\n"
3135 " of voxels at or below each pthr threshold, were +*\n"
3136 " found to estimate at each alpha level. +*\n"
3137 "*+ In other words, the probability that noise-only +*\n"
3138 " data (of the given smoothness) will cause +*\n"
3139 " above-threshold (at the given pthr) clusters is +*\n"
3140 " smaller than the desired alpha levels. +*\n"
3141 "*+ This problem can arise when the masked region +*\n"
3142 " being simulated is small and at the same time +*\n"
3143 " the smoothness (FWHM) is large. +*\n"
3144 "*+ Read the 'CAUTION and CAVEAT' section at the end +*\n"
3145 " of the '-help' output for a longer explanation. +*\n\n");
3146 free(amesg) ; amesg = NULL ;
3147 }
3148
3149 } /* end of loop over mmm == sidedness */
3150
3151 if( do_athr_sum ){
3152 int iathr ; FILE *fp ; char fname[256] ;
3153
3154 thresh_summer_athr(athr_sum_bot,athr_sum_top) ;
3155
3156 sprintf(fname,"%s.sumup.1D",(prefix!=NULL)?prefix:"ClustSim") ;
3157 fp = fopen(fname,"w") ;
3158 fprintf(fp,"# %s\n",commandline) ;
3159 fprintf(fp,"# "
3160 "1dplot -one -box -xaxis 0:0.1:10:5 -yaxis 0:0.32:16:2"
3161 " -xlabel 'nominal \\alpha at fixed p'"
3162 " -ylabel 'integrated FAR over p\\in[%g,%g]' "
3163 " -ynames 1s:NN1 1s:NN2 1s:NN3 2s:NN1 2s:NN2 2s:NN3 bs:NN1 bs:NN2 bs:NN3"
3164 " -x %s'[0]' -plabel '\\noesc %s' %s'[1..$]'\n" ,
3165 pthr[athr_sum_top] , pthr[athr_sum_bot] , fname,fname,fname ) ;
3166 fprintf(fp,"# alpha ") ;
3167 fprintf(fp," 1s:NN1 1s:NN2 1s:NN3") ;
3168 if( do_athr_sum > 1 ){
3169 fprintf(fp," 2s:NN1 2s:NN2 2s:NN3") ;
3170 fprintf(fp," bs:NN1 bs:NN2 bs:NN3") ;
3171 }
3172 fprintf(fp,"\n") ;
3173 for( iathr=0 ; iathr < nathr ; iathr++ ){
3174 fprintf(fp," %.4f ",athr[iathr]) ;
3175 fprintf(fp," %.4f %.4f %.4f",
3176 rfa_1sid_NN1[iathr],
3177 rfa_1sid_NN2[iathr], rfa_1sid_NN3[iathr]) ;
3178 if( do_athr_sum > 1 ){
3179 fprintf(fp," %.4f %.4f %.4f",
3180 rfa_2sid_NN1[iathr],
3181 rfa_2sid_NN2[iathr], rfa_2sid_NN3[iathr]) ;
3182 fprintf(fp," %.4f %.4f %.4f",
3183 rfa_bsid_NN1[iathr],
3184 rfa_bsid_NN2[iathr], rfa_bsid_NN3[iathr]) ;
3185 }
3186 fprintf(fp,"\n") ;
3187 }
3188 fclose(fp) ;
3189 }
3190
3191 #ifdef USE_SHAVE
3192 destroy_shave() ;
3193 #endif
3194
3195 #if 0
3196 if( verb )
3197 INFO_message("Clock time now = %.1f s",COX_clock_time()) ;
3198 #endif
3199
3200 } /* end of outputizationing */
3201
3202 /*------- a minor aid for the pitiful helpless user [e.g., me] -------*/
3203
3204 if( refit_cmd != NULL ){
3205 FILE *fp ;
3206 INFO_message("Command fragment to put cluster results into a dataset header;\n"
3207 " + (also echoed to file %s for your scripting pleasure)\n"
3208 " + Append the name of the datasets to be patched to this command:\n"
3209 " %s" , cmd_fname , refit_cmd ) ;
3210 fp = fopen(cmd_fname,"w") ;
3211 if( fp == NULL ){
3212 ERROR_message("Can't write 3drefit command fragment to file %s :-(",cmd_fname) ;
3213 } else {
3214 fprintf(fp,"%s\n",refit_cmd) ; fclose(fp) ;
3215 }
3216 }
3217
3218 /*-------- run away screaming into the night ----- AAUUGGGHHH!!! --------*/
3219
3220 exit(0);
3221 }
3222