1 #include "mrilib.h"
2 #include <ctype.h>
3
4 /*----------------------------------------------------------------------------*/
5 /***** This is where the ETAC algorithm is implemented (God help us all). *****/
6 /*----------------------------------------------------------------------------*/
7
8 #undef DEBUG_MEM
9 #if 0 && defined(ALLOW_MCW_MALLOC)
10 # define DEBUG_MEM
11 #endif
12
13 #ifdef USE_OMP
14 # include <omp.h>
15 #endif
16
17 #ifdef DEBUG_MEM
18 # include "mcw_malloc.c"
19 # define malloc(a) mcw_malloc((a),__FILE__,__LINE__)
20 # define realloc(a,b) mcw_realloc((a),(b),__FILE__,__LINE__)
21 # define calloc(a,b) mcw_calloc((a),(b),__FILE__,__LINE__)
22 # define free(a) mcw_free((a),__FILE__,__LINE__)
23 # undef strdup
24 # define strdup(a) mcw_strdup((a),__FILE__,__LINE__)
25 #endif
26
27 #include "mri_threshX.c" /* important stuff: clustering, thresholding */
28 #include "thd_Xdataset.c" /* input dataset format and funcs */
29
30 #undef MEMORY_CHECK
31 #ifdef DEBUG_MEM
32 # define MEMORY_CHECK(mm) \
33 do{ if( verb > 1 ) \
34 ININFO_message("Memory %s: %s",mm,mcw_malloc_status("3dXClustSim.c",__LINE__)) ; \
35 } while(0)
36 #else
37 # define MEMORY_CHECK(mm) /*nada*/
38 #endif
39
40 /*---------------------------------------------------------------------------*/
41 /*--- Global data ---*/
42
43 static THD_3dim_dataset *mask_dset = NULL ; /* mask dataset (required) */
44 static byte *mask_vol = NULL; /* mask dataset volume */
45 static int mask_nvox = 0, mask_ngood = 0; /* number of good voxels in mask volume */
46
47 #define INMASK(ijk) (mask_vol[ijk]!=0) /* is a voxel in the mask? */
48
49 /* 3D indexes for each point in the mask */
50
51 static ind_t *ipmask=NULL , *jpmask=NULL , *kpmask=NULL ;
52 static int *ijkmask=NULL ;
53
54 /* map from 1D index in volume to points in the mask */
55
56 static int *ijk_to_vec=NULL ;
57
58 static Xdataset *xinset=NULL ; /* global struct of input dataset(s) */
59
60 static int nx ; /* 3D grid stuff */
61 static int ny ;
62 static int nz ;
63 static int nxy ;
64 static int nxyz ;
65 static int nxyz1 ;
66 static float dx ;
67 static float dy ;
68 static float dz ;
69
70 static int niter = 0 ; /* number of iterations (realizations per case) */
71 static float nit33 = 0.0f ;
72
73 static int nsplit = 0 ; /* 24 Jan 2019 */
74 static float split_frac = 0.0f ;
75 static int niter_clust = 0 ;
76 static int niter_tbot = 0 ;
77 static int niter_ttop = 0 ;
78 static int niter_test = 0 ;
79
80 static int ncase = 1 ; /* number of cases (e.g., blurs) */
81 static char **lcase = NULL ; /* labels for cases */
82 static int nsim = 0 ; /* ncase * niter_clust = number of sim volumes */
83
84 /*-- stuff for extra outputs (don't enable this without good reason) --*/
85
86 #undef ALLOW_EXTRAS
87
88 #ifdef ALLOW_EXTRAS
89 static int do_FOMcount = 0 ;
90 static int do_FARvox = 0 ;
91 #else
92 # define do_FOMcount 0
93 # define do_FARvox 0
94 #endif
95
96 /* p-value thresholds, etc. */
97
98 #define PMAX 0.5
99
100 static int npthr = 10 ;
101 static double *pthr = NULL ;
102 static double pthr_init[10] = { 0.010, 0.009, 0.008, 0.007, 0.006, 0.005, 0.004, 0.003, 0.002, 0.001 } ;
103
104 static float *zthr_1sid = NULL ;
105 static float *zthr_2sid = NULL ;
106 static float *zthr_used = NULL ;
107
108 /* hpow values: FOM = sum |z|^H for H=0, 1, and/or 2 */
109
110 static int do_hpow0 = 0 ; /* 1 or 0 */
111 static int do_hpow1 = 0 ;
112 static int do_hpow2 = 1 ;
113 static int nhpow = 1 ; /* sum of the above */
114
115 /** athr would be for computing threshold for more than one goal **/
116
117 #if 0
118 static int nathr = 5 ;
119 static double *athr = NULL ;
120 static double athr_init[5] = { 0.05 , 0.04 , 0.03 , 0.02 , 0.01 } ;
121 #endif
122
123 static int verb = 1 ;
124 static int nthr = 1 ; /* default number of threads */
125
126 /* These defaults are also set in 3dttest++.c and
127 if changed in one program, must be changed in the other! */
128
129 static int nnlev = 2 ; /* NN = 1 or 2 or 3 */
130 static int nnsid = 2 ; /* sidedness of test = 1 or 2 */
131
132 /* macro for use inside parallel regions,
133 to declare and get ithr = thread number */
134
135 #undef DECLARE_ithr
136 #ifdef USE_OMP
137 # define DECLARE_ithr const int ithr=omp_get_thread_num()
138 #else
139 # define DECLARE_ithr const int ithr=0
140 #endif
141
142 static const int min_mask = 1024 ; /* min # voxels */
143 static const int min_nvol = 10000 ; /* min # realizations */
144
145 static float dilate_fac = 0.0666f ; /* for STEP 2 */
146
147 static char *prefix = "Xsim.nii" ;
148
149 static MRI_IMAGE *imtemplate = NULL ;
150
151 #undef PSMALL
152 #define PSMALL 1.e-15
153
154 #define FARP_GOAL 5.00f /* 5 percent */
155 #define FGFAC 1.00f /* fudge factor (1 = no fudge for you) */
156
157 static float fgfac = FGFAC ;
158 static float farp_goal = FARP_GOAL ;
159
160 /* lines directly below are also in 3dttest++.c
161 only change them here if you change them there as well! */
162 #define NFARP 9
163 #define NMUCHO 25 /* expand gamut [15 Jan 2021] */
164 static float farplist[NMUCHO] =
165 { 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f,
166 10.f, 11.f, 12.f, 13.f, 14.f, 15.f, 16.f, 17.f, 18.f, 19.f,
167 20.f, 21.f, 22.f, 23.f, 24.f, 25.f } ;
168 static float min_fgoal = 1.0f ;
169 static float max_fgoal = 25.0f ;
170 static char *abcd[NMUCHO] = /* just for fun */
171 { "a", "b", "c", "d", "e", "f", "g", "h", "i",
172 "j", "k", "l", "m", "n", "o", "p", "q", "r", "s",
173 "t", "u", "v", "w", "x", "y" } ;
174
175 static int do_multifarp = 0 ;
176 static int numfarp = 1 ;
177
178 #define FG_GOAL (farp_goal*fgfac)
179 #define MAXITE 11
180
181 static int do_global_etac = 1 ; /* Sep 2018 */
182 static int do_local_etac = 0 ;
183
184 #define TOPFRAC 0.2468f
185
186 /* [08 Mar 2019] eliminate multiple clusters from same indx by sorting */
187 static int use_indx = 1 ;
188
189 void qsort_IFFF( int n, int *iar, float *xar, float *yar, float *zar ) ;
190 int ixmax_IFFF( int n, int *iar, float *xar, float *yar, float *zar,
191 float *xout, float *yout, float *zout ) ;
192
193 /*----------------------------------------------------------------------------*/
194 /*! Threshold for upper tail probability of N(0,1) = 1-inverseCDF(pval) */
195
zthresh(double pval)196 double zthresh( double pval )
197 {
198 if( pval <= 0.0 ) pval = PSMALL ;
199 else if( pval >= 1.0 ) pval = 1.0 - PSMALL ;
200 return qginv(pval) ; /* mri_stats.c */
201 }
202
203 /*---------------------------------------------------------------------------*/
204 /* For use in progress counts [currently there is no progress] */
205
206 #if 0
207 static int vsnn=0 ;
208
209 static void vstep_reset(void){ vsnn=0; }
210
211 static void vstep_print(void)
212 {
213 static char xx[10] = "0123456789" ;
214 fprintf(stderr , "%c" , xx[vsnn%10] ) ;
215 if( vsnn%10 == 9) fprintf(stderr,".") ;
216 vsnn++ ;
217 }
218 #endif
219
220 /*---------------------------------------------------------------------------*/
221 /* Routine to initialize the input options (values are in global variables). */
222
get_options(int argc,char ** argv)223 void get_options( int argc , char **argv )
224 {
225 char *ep , msg[1024] ;
226 int nopt=1 , ii ;
227
228 ENTRY("get_options") ;
229
230 while( nopt < argc ){ /* scan until all args processed */
231
232 /*----- -ncase N lab1 lab2 ... labN -----*/
233
234 if( strcasecmp(argv[nopt],"-ncase") == 0 ){ /* for multiple blurs (e.g.) */
235 char *thisopt=argv[nopt] ;
236 if( ++nopt >= argc )
237 ERROR_exit("No argument after '%s' :(",thisopt) ;
238 ncase = (int)strtod(argv[nopt],NULL) ;
239 if( ncase < 1 )
240 ERROR_exit("Value '%s' after option '%s' is illegal :(",argv[nopt],thisopt) ;
241 if( ncase > 9 )
242 WARNING_message("Value %d after option '%s' is very big (> 9)",ncase,thisopt) ;
243 if( nopt+ncase >= argc )
244 ERROR_exit("Need %d label arguments after option '%s %d'",ncase,thisopt,ncase) ;
245 lcase = (char **)malloc(sizeof(char *)*ncase) ;
246 for( ii=0 ; ii < ncase ; ii++ ){
247 lcase[ii] = strdup(argv[nopt+ii+1]) ;
248 if( lcase[ii][0] == '-' )
249 WARNING_message("label '%s' after option '%s' looks like another option :(",
250 lcase[ii],thisopt) ;
251 else if( lcase[ii][0] == '\0' || isspace(lcase[ii][0]) ){
252 WARNING_message("label #%d after option '%s' is empty :(",
253 ii+1,thisopt) ;
254 }
255 }
256 nopt += ncase+1 ; continue ;
257 }
258
259 /*----- -dilfac ddd -----*/
260
261 if( strcasecmp(argv[nopt],"-dilfac") == 0 ){ /* dilation factor */
262 if( ++nopt >= argc )
263 ERROR_exit("You need 1 argument after option '-dilfac'") ;
264 dilate_fac = (float)strtod(argv[nopt],NULL) ;
265 nopt++ ; continue ;
266 }
267
268 /*----- -NN 1 or 2 or 3 -----*/
269
270 if( strcasecmp(argv[nopt],"-NN") == 0 ){
271 if( ++nopt >= argc )
272 ERROR_exit("You need 1 argument after option '-NN'") ;
273 nnlev = (int)strtod(argv[nopt],NULL) ;
274 if( nnlev < 1 || nnlev > 3 )
275 ERROR_exit("-NN must be 1 or 2 or 3 :(") ;
276 nopt++ ; continue ;
277 }
278 if( strcasecmp(argv[nopt],"-NN1") == 0 || -strcasecmp(argv[nopt],"-1NN") == 0 ){
279 nnlev = 1 ; nopt++ ; continue ;
280 }
281 if( strcasecmp(argv[nopt],"-NN2") == 0 || -strcasecmp(argv[nopt],"-2NN") == 0 ){
282 nnlev = 2 ; nopt++ ; continue ;
283 }
284 if( strcasecmp(argv[nopt],"-NN3") == 0 || -strcasecmp(argv[nopt],"-3NN") == 0 ){
285 nnlev = 3 ; nopt++ ; continue ;
286 }
287
288 /*----- -sid 1 or 2 -----*/
289
290 if( strcasecmp(argv[nopt],"-sid") == 0 ){
291 if( ++nopt >= argc )
292 ERROR_exit("You need 1 argument after option '-sid'") ;
293 nnsid = (int)strtod(argv[nopt],NULL) ;
294 if( nnsid < 1 || nnsid > 2 )
295 ERROR_exit("-sid must be 1 or 2 :(") ;
296 nopt++ ; continue ;
297 }
298 if( strcasecmp(argv[nopt],"-sid1") == 0 || strcasecmp(argv[nopt],"-1sid") == 0 ){
299 nnsid = 1 ; nopt++ ; continue ;
300 }
301 if( strcasecmp(argv[nopt],"-sid2") == 0 || strcasecmp(argv[nopt],"-2sid") == 0 ){
302 nnsid = 2 ; nopt++ ; continue ;
303 }
304
305 /*----- -inset mask sdata ... -----*/
306
307 if( strcasecmp(argv[nopt],"-inset" ) == 0 || /* no dataset inputs */
308 strcasecmp(argv[nopt],"-insdat") == 0 ){ /* just sdat files */
309 int nfile ;
310 if( xinset != NULL )
311 ERROR_exit("You can't use option '%s' more than once!",argv[nopt]) ;
312 if( ++nopt >= argc-1 )
313 ERROR_exit("You need at least 2 arguments after option '%s'",argv[nopt-1]) ;
314
315 for( ii=nopt ; ii < argc && argv[ii][0] != '-' ; ii++ ) ; /*nada*/
316 nfile = ii-nopt ;
317
318 if( verb )
319 INFO_message("Loading %s datasets",argv[nopt-1]) ;
320
321 /* see thd_Xdataset.c */
322 xinset = open_Xdataset( argv[nopt], nfile-1,argv+(nopt+1) ) ;
323
324 if( xinset->ngood < min_mask )
325 ERROR_exit("mask has %d good voxels; minimum allowed is %d",
326 xinset->ngood , min_mask ) ;
327
328 if( xinset->nvtot < min_nvol )
329 WARNING_message("only %d input volumes, less than minimum of %d",
330 xinset->nvtot,min_nvol) ;
331
332 nopt += nfile ; continue ;
333 }
334
335 /*----- -prefix -----*/
336
337 if( strcasecmp(argv[nopt],"-prefix") == 0 ){
338 nopt++ ; if( nopt >= argc ) ERROR_exit("need argument after -prefix!") ;
339 prefix = strdup(argv[nopt]) ;
340 if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix option!") ;
341 nopt++ ; continue ;
342 }
343
344 /*---- -quiet and -verb ----*/
345
346 if( strcasecmp(argv[nopt],"-quiet") == 0 ){ /* why would you want this? */
347 verb = 0 ; nopt++ ; continue ;
348 }
349
350 if( strcasecmp(argv[nopt],"-verb") == 0 ){ /* more fun fun fun! */
351 verb++ ; nopt++ ; continue ;
352 }
353
354 /*----- -sorter and -nosorter [HIDDEN] -----*/
355
356 if( strcasecmp(argv[nopt],"-sorter") == 0 ){
357 use_indx = 1 ; nopt++ ; continue ;
358 }
359
360 if( strcasecmp(argv[nopt],"-nosorter") == 0 ){
361 use_indx = 0 ; nopt++ ; continue ;
362 }
363
364 #ifdef ALLOW_EXTRAS
365 /*---- -FOMcount ----*/
366
367 if( strcasecmp(argv[nopt],"-FOMcount") == 0 ){
368 do_FOMcount = 1 ; nopt++ ; continue ;
369 }
370
371 /*---- -noFARvox ----*/
372
373 if( strcasecmp(argv[nopt],"-noFARvox") == 0 ){
374 do_FARvox = 0 ; nopt++ ; continue ;
375 }
376
377 /*---- -FARvox ----*/
378
379 if( strcasecmp(argv[nopt],"-FARvox") == 0 ){
380 do_FARvox = 1 ; nopt++ ; continue ;
381 }
382 #endif
383
384 /*----- -pthr p p p -----*/
385
386 #define IS_NUMERIC(sss) \
387 ( ( isdigit((sss)[0]) ) || \
388 ( (sss)[0] == '.' && isdigit((sss)[1]) ) )
389
390
391 if( strcmp(argv[nopt],"-pthr") == 0 || strcmp(argv[nopt],"-pval") == 0 ){
392
393 if( pthr != NULL ) ERROR_exit("you can't use -pthr twice!") ;
394 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
395
396 /* scan to see how many numeric args follow */
397
398 for( ii=nopt ; ii < argc && IS_NUMERIC(argv[ii]) ; ii++ ) ; /*nada*/
399 npthr = ii-nopt ;
400 if( npthr <= 0 ) ERROR_exit("No positive values found after %s",argv[nopt-1]) ;
401
402 /* re-scan to load numeric values */
403
404 pthr = (double *)realloc(pthr,sizeof(double)*npthr) ;
405 for( ii=0 ; ii < npthr ; ii++ ){
406 pthr[ii] = strtod(argv[nopt+ii],NULL) ;
407 if( pthr[ii] <= 0.0 || pthr[ii] > PMAX )
408 ERROR_exit("value '%s' after '%s' is illegal!",argv[nopt+ii],argv[nopt-1]) ;
409 }
410
411 if( npthr > 1 ){
412 for( ii=0 ; ii < npthr ; ii++ ) pthr[ii] = -pthr[ii] ; /* sort into */
413 qsort_double( npthr , pthr ) ; /* descending */
414 for( ii=0 ; ii < npthr ; ii++ ) pthr[ii] = -pthr[ii] ; /* order */
415 for( ii=1 ; ii < npthr ; ii++ ){
416 if( pthr[ii] == pthr[ii-1] )
417 WARNING_message("duplicate value %g after '%s'",pthr[ii],argv[nopt-1]) ;
418 }
419 }
420
421 nopt += npthr ;
422 continue ;
423 }
424
425 /*----- -hpow h h h [28 Mar 2017] -----*/
426
427 if( strcmp(argv[nopt],"-hpow") == 0 ){
428 int jj ;
429 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
430
431 /* scan for integers 0 and/or 1 and/or 2 */
432
433 do_hpow0 = do_hpow1 = do_hpow2 = 0 ;
434 for( ; nopt < argc && isdigit(argv[nopt][0]) ; nopt++ ){
435 switch( (int)strtod(argv[nopt],NULL) ){
436 case 0: do_hpow0 = 1 ; break ;
437 case 1: do_hpow1 = 1 ; break ;
438 case 2: do_hpow2 = 1 ; break ;
439 default:
440 ERROR_exit("value '%s' after '-hpow' is illegal (only 0, 1, and/or 2 allowed)",argv[nopt]) ;
441 }
442 }
443
444 nhpow = do_hpow0 + do_hpow1 + do_hpow2 ;
445 if( nhpow < 1 )
446 ERROR_exit("No valid powers (0 and/or 1 and/or 2) found after '-hpow' :(") ;
447
448 continue ;
449 }
450
451 /*----- -multiFPR [23 Aug 2017] -----*/
452
453 if( strcasecmp(argv[nopt],"-multiFPR") == 0 ){
454 do_multifarp = 1 ;
455 nopt++ ; continue ;
456 }
457
458 if( strcasecmp(argv[nopt],"-muchoFPR") == 0 ){
459 do_multifarp = 2 ;
460 nopt++ ; continue ;
461 }
462
463 /*----- -minclust [21 Sep 2017] -----*/
464
465 if( strcasecmp(argv[nopt],"-minclust") == 0 ){
466 int mc ;
467 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
468 mc = (int)strtod(argv[nopt],NULL) ;
469 if( mc < 2 )
470 WARNING_message("value %d after %s is illegal",mc,argv[nopt-1]) ;
471 else
472 set_Xmin_clust(mc) ;
473
474 nopt++ ; continue ;
475 }
476
477
478 /*----- -FPR xx [23 Aug 2017] -----*/
479
480 if( strcasecmp(argv[nopt],"-FPR") == 0 ){
481 float fgoal ;
482 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
483 do_multifarp = 0 ;
484 fgoal = (float)rint(strtod(argv[nopt],NULL)) ;
485 if( fgoal < min_fgoal ){
486 WARNING_message("fpr=%.1f%% too small : setting fpr=%.1f%%",fgoal,min_fgoal) ;
487 fgoal = min_fgoal ;
488 } else if( fgoal > max_fgoal ){
489 WARNING_message("fpr=%.1f%% too large : setting fpr=%.1f%%",fgoal,max_fgoal) ;
490 fgoal = max_fgoal ;
491 }
492 farp_goal = fgoal ;
493 nopt++ ; continue ;
494 }
495
496 /*----- global and local [Sep 2018] -----*/
497
498 if( strcasecmp(argv[nopt],"-noglobal") == 0 ){ /* Sep 2018 */
499 do_global_etac = 0 ; nopt++ ; continue ;
500 }
501 if( strcasecmp(argv[nopt],"-nolocal") == 0 ){
502 do_local_etac = 0 ; nopt++ ; continue ;
503 }
504 if( strcasecmp(argv[nopt],"-global") == 0 ){
505 do_global_etac = 1 ; nopt++ ; continue ;
506 }
507 if( strcasecmp(argv[nopt],"-local") == 0 ){
508 do_local_etac = 1 ; nopt++ ; continue ;
509 }
510
511 /*----- split fraction [24 Jan 2019] -----*/
512
513 if( strcasecmp(argv[nopt],"-splitfrac") == 0 ){
514 nopt++; if( nopt >= argc ) ERROR_exit("need argument after %s",argv[nopt-1]);
515 split_frac = (float)strtod(argv[nopt],NULL) ;
516 nopt++ ; continue ;
517 }
518
519 /*----- unknown option -----*/
520
521 ERROR_exit("3dXClustSim -- unknown option '%s' :(",argv[nopt]) ;
522
523 } /* end loop over command line args */
524
525 if( !do_global_etac && !do_local_etac )
526 ERROR_exit("3dXClustSim: global and local calculations both turned off?! :(") ;
527
528 /*-- sneaky way to change default parameters --*/
529
530 fgfac = AFNI_numenv("AFNI_XCLUSTSIM_FGFAC") ;
531 if( fgfac < 0.1f || fgfac > 1.0f ) fgfac = FGFAC ;
532
533 switch( do_multifarp ){
534 default:
535 case 0:
536 numfarp = 1 ; farplist[0] = farp_goal ;
537 INFO_message("Single FPR goal: %.1f%%",farp_goal) ;
538 break ;
539
540 case 1:
541 case 2:
542 numfarp = (do_multifarp==1) ? NFARP : NMUCHO ; /* 15 Jan 2021 */
543 strcpy(msg,"Multiple FPR goals:") ;
544 for( ii=0 ; ii < numfarp ; ii++ ) /* quadruple % reduced to double % below */
545 sprintf( msg+strlen(msg) , " %.1f%%%%" , farplist[ii] ) ;
546 INFO_message(msg) ; /* and then to single % in this output */
547 break ;
548 }
549
550 /*------- finalize some simple setup stuff --------*/
551
552 if( xinset == NULL ) ERROR_exit("-inset or -insdat option is mandatory :(") ;
553
554 /* create default pthr array, if none was given */
555
556 if( pthr == NULL ){
557 pthr = (double *)malloc(sizeof(double)*npthr) ;
558 AAmemcpy( pthr , pthr_init , sizeof(double)*npthr ) ;
559 if( verb > 1 )
560 INFO_message("Using default %d p-value thresholds",npthr) ;
561 }
562
563 /* Show the p-value thresholds [08 Mar 2019] */
564
565 strcpy(msg,"p-value thresholds:") ;
566 for( ii=0 ; ii < npthr ; ii++ )
567 sprintf( msg+strlen(msg) , " %.4f" , pthr[ii] ) ;
568 INFO_message(msg) ;
569 if( npthr == 1 ) /* should not happen, but you never know about users */
570 WARNING_message("only 1 p-value threshold?!") ;
571
572 /* another tidbit of info */
573
574 INFO_message("min cluster size : %d voxels",min_clust) ; /* 21 Sep 2017 */
575
576 /* create the default label for ncase=1 */
577
578 if( ncase == 1 && lcase == NULL ){
579 lcase = (char **)malloc(sizeof(char *)) ;
580 lcase[0] = strdup("A") ;
581 if( verb > 1 )
582 INFO_message("Using default case label 'A'") ;
583 }
584
585 /* check to see if -ncase is used correctly; set number of iterations */
586
587 if( ncase > 1 ){
588 if( xinset->nvtot % ncase != 0 ){
589 ERROR_exit("number of input volumes=%d not evenly divisible by ncase=%d",
590 xinset->nvtot,ncase) ;
591 } else {
592 niter = xinset->nvtot / ncase ;
593 if( verb > 1 ){
594 INFO_message("number of input volumes = %d",xinset->nvtot) ;
595 ININFO_message("number of cases = %d",ncase) ;
596 ININFO_message("number of vols/case = %d / %d = %d",xinset->nvtot,ncase,niter) ;
597 }
598 }
599 } else {
600 niter = xinset->nvtot ;
601 if( verb > 1 )
602 INFO_message("number of input volumes = %d",niter) ;
603 }
604
605 if( split_frac == 0.0f )
606 split_frac = AFNI_numenv("AFNI_XCLUSTSIM_SPLITFRAC") ;
607
608 if( split_frac >= 0.2f && split_frac <= 0.8f )
609 nsplit = (int)rintf(niter*split_frac) ;
610
611 if( nsplit > 0 ){ /* split clustering and testing [24 Jan 2019] */
612 niter_clust = nsplit ; /* how many to use for clustering */
613 niter_tbot = nsplit ; /* bot and top indexes to use for testing */
614 niter_ttop = niter ;
615 if( verb > 1 )
616 ININFO_message("split testing at volume = %d",nsplit) ;
617 } else { /* the olden way */
618 niter_clust = niter ;
619 niter_tbot = 0 ;
620 niter_ttop = niter ;
621 if( verb > 1 )
622 ININFO_message("split testing = OFF") ;
623 }
624 niter_test = niter_ttop - niter_tbot ; /* number used for testing */
625
626 nsim = niter_clust * ncase ; /* number of simulation volumes */
627
628 #if 0
629 if( athr == NULL ){
630 athr = (double *)malloc(sizeof(double)*nathr) ;
631 AAmemcpy( athr , athr_init , sizeof(double)*nathr ) ;
632 }
633 #endif
634
635 /* load the geometry variables */
636
637 nx = DSET_NX(xinset->mask_dset) ;
638 ny = DSET_NY(xinset->mask_dset) ;
639 nz = DSET_NZ(xinset->mask_dset) ;
640 dx = fabsf(DSET_DX(xinset->mask_dset)) ;
641 dy = fabsf(DSET_DY(xinset->mask_dset)) ;
642 dz = fabsf(DSET_DZ(xinset->mask_dset)) ;
643
644 /* copy some variables because xinset was a late addition to the code */
645
646 mask_dset = xinset->mask_dset ;
647 imtemplate = DSET_BRICK(xinset->mask_dset,0) ; /* sample 3D volume */
648 mask_ngood = xinset->ngood ;
649 mask_vol = xinset->mask_vol ;
650 if( verb > 1 )
651 INFO_message("mask has %d points",mask_ngood) ;
652
653 nxy = nx*ny ; nxyz = nxy*nz ; nxyz1 = nxyz - nxy ;
654 if( nxyz < 1024 )
655 ERROR_exit("Only %d voxels in simulation?! Need at least 1024.",nxyz) ;
656
657 /* make a list of the i,j,k coordinates of each point in the mask */
658
659 { int pp,qq , xx,yy,zz ;
660 ipmask = xinset->ipmask = (ind_t *)malloc(sizeof(ind_t)*mask_ngood) ;
661 jpmask = xinset->jpmask = (ind_t *)malloc(sizeof(ind_t)*mask_ngood) ;
662 kpmask = xinset->kpmask = (ind_t *)malloc(sizeof(ind_t)*mask_ngood) ;
663 ijkmask= xinset->ijkmask ;
664 for( pp=qq=0 ; qq < nxyz ; qq++ ){
665 if( INMASK(qq) ){
666 IJK_TO_THREE(qq,xx,yy,zz,nx,nxy) ;
667 ipmask[pp] = (ind_t)xx; jpmask[pp] = (ind_t)yy; kpmask[pp] = (ind_t)zz;
668 pp++ ;
669 }
670 }
671 ijk_to_vec = xinset->ijk_to_vec = (int *)malloc(sizeof(int)*nxyz) ;
672 for( pp=qq=0 ; qq < nxyz ; qq++ )
673 ijk_to_vec[qq] = INMASK(qq) ? pp++ : -1 ;
674 }
675
676 /*-- z-score thresholds for the various p-values --*/
677
678 zthr_1sid = (float *)malloc(sizeof(float)*npthr) ;
679 zthr_2sid = (float *)malloc(sizeof(float)*npthr) ;
680 for( ii=0 ; ii < npthr ; ii++ ){
681 zthr_1sid[ii] = (float)zthresh( pthr[ii] ) ;
682 zthr_2sid[ii] = (float)zthresh( 0.5*pthr[ii] ) ;
683 }
684 zthr_used = (nnsid==1) ? zthr_1sid : zthr_2sid ;
685
686 EXRETURN ;
687 }
688
689 /*---------------------------------------------------------------------------*/
690 /* Load image (realization) for thresholding and clustering into fim array */
691
generate_image(float * fim,int ind)692 void generate_image( float *fim , int ind )
693 {
694 load_from_Xdataset( xinset , ind , fim ) ;
695 return ;
696 }
697
698 /*---------------------------------------------------------------------------*/
699 /* Global cluster collection:
700 Xclustar_g[icase][ipthr][iter] = cluster array at iteration iter,
701 for threshold ipthr, for case icase. The basic array created in main().
702 *//*-------------------------------------------------------------------------*/
703
704 /* pointer to pointer to pointer to pointer to struct! Isn't C great? */
705
706 static Xcluster_array ****Xclustar_g ;
707
708 /* Xclust_tot[icase][ipthr][nn] =
709 consolidated clusters from all sims,
710 for nn=0 .. nclust_tot[icase][ipthr]-1
711 Usually nclust_tot[][] is less than niter, since not all
712 iterations will have clusters in a given case/thresh combo */
713
714 static int nclust_max ;
715 static int **nclust_tot ; /* [ncase][npthr] */
716 static Xcluster ****Xclust_tot ; /* [ncase][npthr][nclust_tot[icase][ipthr]] */
717
718 /*---------------------------------------------------------------------------*/
719 /* Get a NN{nn}_{ss}sided cluster array at a particular thresh/sub-case
720 (icase), at a particular iteration (iter), and save it into the global
721 cluster collection Xclustar_g. nn=1 or 2 or 3; ss=1 or 2.
722 *//*-------------------------------------------------------------------------*/
723
gather_clusters(int icase,int ipthr,float * fim,MRI_IMAGE * tfim,int nn,int ss,int iter)724 void gather_clusters( int icase, int ipthr,
725 float *fim, MRI_IMAGE *tfim, int nn,int ss, int iter )
726 {
727 register int ii ; register float thr ; float *tfar = MRI_FLOAT_PTR(tfim) ;
728
729 thr = fabsf(zthr_used[ipthr]) ;
730
731 #if 1 /** older code **/
732 if( ss == 1 ){
733 if( iter%2 == 0 ){
734 for( ii=0 ; ii < nxyz ; ii++ )
735 tfar[ii] = (fim[ii] >= thr) ? fim[ii] : 0.0f ;
736 } else {
737 for( ii=0 ; ii < nxyz ; ii++ )
738 tfar[ii] = (fim[ii] <= -thr) ? fim[ii] : 0.0f ;
739 }
740 } else {
741 for( ii=0 ; ii < nxyz ; ii++ )
742 tfar[ii] = (fabsf(fim[ii]) >= thr) ? fim[ii] : 0.0f ;
743 }
744 Xclustar_g[icase][ipthr][iter] = find_Xcluster_array( tfim,nn,iter, 0,NULL,NULL,NULL ) ;
745
746 #else /** new code: 1-sided uses both pos and neg [doesn't work] **/
747 if( ss == 1 ){
748 Xcluster_array *pcar , *ncar ;
749 for( ii=0 ; ii < nxyz ; ii++ )
750 tfar[ii] = (fim[ii] >= thr) ? fim[ii] : 0.0f ;
751 pcar = find_Xcluster_array( tfim,nn,iter, 0,NULL,NULL,NULL ) ;
752 for( ii=0 ; ii < nxyz ; ii++ )
753 tfar[ii] = (fim[ii] <= -thr) ? fim[ii] : 0.0f ;
754 ncar = find_Xcluster_array( tfim,nn,iter, 0,NULL,NULL,NULL ) ;
755 MERGE_Xcluster_arrays(pcar,ncar) ; /* ncar is deleted */
756 Xclustar_g[icase][ipthr][iter] = pcar ;
757 } else {
758 for( ii=0 ; ii < nxyz ; ii++ )
759 tfar[ii] = (fabsf(fim[ii]) >= thr) ? fim[ii] : 0.0f ;
760 Xclustar_g[icase][ipthr][iter] = find_Xcluster_array( tfim,nn,iter, 0,NULL,NULL,NULL ) ;
761 }
762 #endif
763
764 return ;
765 }
766
767 /*---------------------------------------------------------------------------*/
768 /* Dilate a cluster by 1 voxel [don't change FOM values] */
769 /* DALL is defined in mri_threshX.c [included way above] */
770
771 #define DILATE_point(pqr) \
772 do{ int i,j,k , npt=xc->npt; \
773 IJK_TO_THREE(pqr,i,j,k,nx,nxy) ; \
774 if( npt == xc->nall ){ \
775 xc->nall += DALL + xc->nall/4 ; \
776 xc->ip = (ind_t *)realloc(xc->ip,sizeof(ind_t)*xc->nall) ; \
777 xc->jp = (ind_t *)realloc(xc->jp,sizeof(ind_t)*xc->nall) ; \
778 xc->kp = (ind_t *)realloc(xc->kp,sizeof(ind_t)*xc->nall) ; \
779 xc->ijk= (int *) realloc(xc->ijk,sizeof(int) *xc->nall) ; \
780 } \
781 xc->ip[npt] = i; xc->jp[npt] = j; xc->kp[npt] = k; \
782 xc->ijk[npt]= pqr ; xc->npt++ ; \
783 } while(0)
784
785 /*---------------------------------------------------------------------------*/
786 /* per-thread work arrays for dilation;
787 the highest level arrays are created in main(),
788 and the sub-arrays are never free()-ed, to avoid OpenMP degradation */
789
790 static int **dilg_ijk=NULL ;
791 static int *ndilg=NULL ;
792
dilate_Xcluster(Xcluster * xc,int nnlev,int ithr)793 void dilate_Xcluster( Xcluster *xc , int nnlev , int ithr )
794 {
795 int npt,ntry,ii,jj , nx1=nx-1,ny1=ny-1,nz1=nz-1 ;
796 int *ijkar , xx,yy,zz , ijk ;
797
798 int pxpy=+1+nx , pxmy=+1-nx , mxpy=nx-1 , mxmy=-1-nx ; /* for NN2 */
799 int pxpz=+1+nxy, pxmz=+1-nxy, mxpz=nxy-1 , mxmz=-1-nxy ;
800 int pypz=nx+nxy, pymz=nx-nxy, mypz=nxy-nx, mymz=-nx-nxy;
801 int pxpypz = +1+nx+nxy , pxpymz = +1+nx-nxy , /* for NN3 */
802 pxmypz = +1-nx+nxy , pxmymz = +1-nx-nxy ,
803 mxpypz = -1+nx+nxy , mxpymz = -1+nx-nxy ,
804 mxmypz = -1-nx+nxy , mxmymz = -1-nx-nxy ;
805
806 if( xc == NULL || xc->npt == 0 ) return ; /* bad inputs */
807
808 npt = xc->npt ;
809
810 /* create the list of candidate points for dilation
811 (most will not survive the second step: already in cluster) */
812
813 ntry = npt * ( (nnlev<=1) ? 6 : (nnlev==2) ? 18 : 26 ) ;
814
815 if( ntry > ndilg[ithr] ){ /* make it bigger if need be */
816 dilg_ijk[ithr] = (int *)realloc(dilg_ijk[ithr],sizeof(int)*(ntry+64)) ;
817 ndilg [ithr] = ntry+64 ;
818 }
819
820 ijkar = dilg_ijk[ithr] ;
821
822 for( jj=ii=0 ; ii < npt ; ii++ ){ /* candidates = neighbors */
823 xx = xc->ip[ii] ; yy = xc->jp[ii] ; zz = xc->kp[ii] ; /* current pt */
824
825 /* no dilation from outer edges of 3D grid */
826 if( xx==0 || xx==nx1 || yy==0 || yy==ny1 || zz==0 || zz==nz1 ) continue ;
827
828 ijk = xx + nx*yy + nxy*zz ; /* 3D grid index */
829
830 ijkar[jj++] = ijk+1 ; ijkar[jj++] = ijk-1 ; /* NN1 dilation */
831 ijkar[jj++] = ijk+nx ; ijkar[jj++] = ijk-nx ; /* candidates */
832 ijkar[jj++] = ijk+nxy ; ijkar[jj++] = ijk-nxy ;
833 if( nnlev >= 2 ){
834 ijkar[jj++] = ijk+pxpy ; ijkar[jj++] = ijk+mxpy ; /* NN2 */
835 ijkar[jj++] = ijk+pxmy ; ijkar[jj++] = ijk+mxmy ;
836 ijkar[jj++] = ijk+pxpz ; ijkar[jj++] = ijk+mxpz ;
837 ijkar[jj++] = ijk+pxmz ; ijkar[jj++] = ijk+mxmz ;
838 ijkar[jj++] = ijk+pypz ; ijkar[jj++] = ijk+mypz ;
839 ijkar[jj++] = ijk+pymz ; ijkar[jj++] = ijk+mymz ;
840 if( nnlev >=3 ){
841 ijkar[jj++] = ijk+pxpypz ; ijkar[jj++] = ijk+pxpymz ; /* NN3 */
842 ijkar[jj++] = ijk+pxmypz ; ijkar[jj++] = ijk+pxmymz ;
843 ijkar[jj++] = ijk+mxpypz ; ijkar[jj++] = ijk+mxpymz ;
844 ijkar[jj++] = ijk+mxmypz ; ijkar[jj++] = ijk+mxmymz ;
845 }}
846 }
847
848 ntry = jj ;
849 qsort_int(ntry,ijkar) ; /* sorting makes it easy to skip duplicates */
850
851 /* for each new point:
852 compare to all other points in the cluster
853 if it isn't any of them, add it to the cluster */
854
855 for( jj=0 ; jj < ntry ; jj++ ){
856 ijk = ijkar[jj] ; if( !INMASK(ijk) ) continue ; /* not a useful pt */
857 if( jj > 0 && ijk==ijkar[jj-1] ) continue ; /* is duplicate */
858 for( ii=0 ; ii < xc->npt ; ii++ ){ /* check if already in cluster */
859 if( ijk == xc->ijk[ii] ) break ; /* oops, already present */
860 }
861 if( ii == xc->npt ) DILATE_point(ijk) ; /* add it in, Babee! */
862 } /* note: xc->npt will expand as the cluster grows */
863
864 return ; /* cluster is all dilated now */
865 }
866
867 #undef DILATE_point
868
869 /*---------------------------------------------------------------------------*/
870 /* Vector struct.
871 One for each point in the mask, to hold the FOM values found at that point.
872 *//*-------------------------------------------------------------------------*/
873
874 typedef struct {
875 ind_t ip,jp,kp ; int ijk ; /* where is this pt? */
876 int npt , nall ; /* array dimensions */
877 float *far ; /* array of FOM values */
878 float val ; /* summary FOM value */
879 } Xvector ;
880
881 #define CREATE_Xvector(xv,siz) \
882 do{ xv = (Xvector *)malloc(sizeof(Xvector)) ; \
883 xv->npt = 0 ; xv->nall = (siz) ; \
884 xv->far = (float *)malloc(sizeof(float)*(siz)) ; \
885 } while(0)
886
887 #define DESTROY_Xvector(xv) \
888 do{ free(xv->far); free(xv); } while(0)
889
890 #define ADDTO_Xvector(xv,val) \
891 do{ if( xv->npt >= xv->nall ){ \
892 xv->nall = xv->npt + DALL + xv->nall/4 ; \
893 xv->far = (float *)realloc(xv->far,sizeof(float)*xv->nall) ; \
894 } \
895 xv->far[xv->npt++] = (val) ; \
896 } while(0)
897
898 /* global FOM vector arrays -- created in main() */
899
900 /* these are workspaces (for hpow=0,1,2) */
901
902 static Xvector **fomvec0 = NULL ; /* [mask_ngood] */
903 static Xvector **fomvec1 = NULL ; /* [mask_ngood] */
904 static Xvector **fomvec2 = NULL ; /* [mask_ngood] */
905
906 /* these are sorted FOM vectors for each case, each pthr, each point */
907
908 static Xvector ****fomsort0 = NULL ; /* [ncase][npthr][mask_ngood] */
909 static Xvector ****fomsort1 = NULL ; /* [ncase][npthr][mask_ngood] */
910 static Xvector ****fomsort2 = NULL ; /* [ncase][npthr][mask_ngood] */
911
912 /*---------------------------------------------------------------------------*/
913 /* Process the clusters to FOM vectors for a range of 3D indexes,
914 for a given p-value thresh/subcase index icase.
915 This function is intended to be use in multiple threads, operating
916 over differing ijkbot..ijktop blocks:
917 A given fomvec will not be in thread conflict, since it is linked
918 to a given voxel in space (ijk), so there is no potential problem
919 when updating that fomvec, or the cluster's ijk entry.
920 *//*-------------------------------------------------------------------------*/
921
process_clusters_to_Xvectors(int ijkbot,int ijktop,int icase,int ipthr)922 void process_clusters_to_Xvectors( int ijkbot, int ijktop , int icase, int ipthr )
923 {
924 Xcluster **xcar = Xclust_tot[icase][ipthr] ; /* clusters we load from */
925 int ncar = nclust_tot[icase][ipthr] ; /* number of cluster arrays */
926 Xcluster *xc ;
927 int cc , pp,npt , ijk,vin ;
928 const int nx1=nx-1,ny1=ny-1,nz1=nz-1 ;
929
930 /* macro to add FOM values to the fomvecH vectors for voxel #iii */
931
932 #define PROCESS(iii) \
933 do{ if( iii >= ijkbot && ijk <= ijktop ){ \
934 vin = ijk_to_vec[iii] ; \
935 if( vin >= 0 ){ \
936 ADDTO_Xvector(fomvec0[vin],xc->fomh[0]) ; \
937 ADDTO_Xvector(fomvec1[vin],xc->fomh[1]) ; \
938 ADDTO_Xvector(fomvec2[vin],xc->fomh[2]) ; \
939 } \
940 }} while(0)
941
942 for( cc=0 ; cc < ncar ; cc++ ){ /* loop over clusters */
943 xc = xcar[cc] ; if( xc == NULL ) continue ;
944 if( xc->ijkmin > ijktop ) continue ;
945 if( xc->ijkmax < ijkbot ) continue ;
946 npt = xc->npt ; /* number of points in cluster */
947 for( pp=0 ; pp < npt ; pp++ ){ /* loop over pts inside cluster */
948 ijk = xc->ijk[pp] ; /* index of pt in grid */
949 PROCESS(ijk) ; /* if pt in region and is good, add FOM to this list */
950
951 #ifdef LARGER_FOOTPRINT /* not enabled -- dilation is used instead */
952 { int iqq ;
953 if( xc->ip[pp] < nx1 ){ iqq = ijk+1 ; PROCESS(iqq) ; }
954 if( xc->ip[pp] > 0 ){ iqq = ijk-1 ; PROCESS(iqq) ; }
955 if( xc->jp[pp] < ny1 ){ iqq = ijk+nx ; PROCESS(iqq) ; }
956 if( xc->jp[pp] > 0 ){ iqq = ijk-nx ; PROCESS(iqq) ; }
957 if( xc->kp[pp] < nz1 ){ iqq = ijk+nxy ; PROCESS(iqq) ; }
958 if( xc->kp[pp] > 0 ){ iqq = ijk-nxy ; PROCESS(iqq) ; }
959 }
960 #endif
961
962 } /* end of loop over pts inside cluster */
963 } /* end of loop over clusters */
964
965 #undef PROCESS
966
967 return ;
968 }
969
970 /*--------------------------------------------------------------------------*/
971 /* get median of number of counts in voxel fomvec's touched by this cluster */
972 /* that is, a measure of how often the voxels in this */
973 /* cluster are 'hit' by other clusters in the simulation */
974
975 static int *ncount_g ; /* workspaces initialized in main() */
976 static float **fcount_g ;
977
get_Xcluster_nbcount(Xcluster * xc,int ithr)978 int get_Xcluster_nbcount( Xcluster *xc , int ithr )
979 {
980 int ii ; float *fc ;
981
982 /* trivial cases */
983
984 if( xc == NULL || xc->npt == 0 ) return 0 ;
985 if( xc->npt == 1 ){
986 xc->nbcount = fomvec0[ijk_to_vec[xc->ijk[0]]]->npt ;
987 return xc->nbcount ;
988 }
989
990 /* size of this cluster is bigger than space we've set aside? */
991
992 if( ncount_g[ithr] < xc->npt ){
993 ncount_g[ithr] = xc->npt + DALL ;
994 fcount_g[ithr] = (float *)realloc(fcount_g[ithr],
995 sizeof(float)*ncount_g[ithr]) ;
996 }
997
998 /* assemble the count of FOM vector hits for each voxel in this cluster */
999
1000 fc = fcount_g[ithr] ; /* thread-specific temp array */
1001 for( ii=0 ; ii < xc->npt ; ii++ )
1002 fc[ii] = fomvec0[ijk_to_vec[xc->ijk[ii]]]->npt ; /* count */
1003
1004 xc->nbcount = (int)qmed_float(xc->npt,fc) ; /* find the median count */
1005 return xc->nbcount ;
1006 }
1007
1008 /*---------------------------------------------------------------------------*/
1009 /* find the intermediate value of x at alphat (t=target),
1010 given x at alpha0 and alpha1, where alpha is a tail (1-CDF) probability,
1011 assumed to be of the extreme value form
1012 alpha(x) = 1 - exp( -exp(-(x-p)/q) )
1013 log(-log(1-alpha(x))) = -x/q + p/q
1014 Define a(alpha) = log(-log(1-alpha)), so
1015 a(x) = -x/q + p/q or x(a) = p - q*a
1016 The alpha values are between 0 and 1; alphat is between alpha0 and alpha1.
1017 *//*-------------------------------------------------------------------------*/
1018
inverse_interp_extreme(float alpha0,float alpha1,float alphat,float x0,float x1)1019 static float inverse_interp_extreme( float alpha0, float alpha1, float alphat,
1020 float x0 , float x1 )
1021 {
1022 float a0,a1,at , xx ;
1023
1024 if( (alphat-alpha0)*(alphat-alpha1) < 0.0f ){ /* alphat is bracketed */
1025 a0 = logf(-logf(1.0f-alpha0)) ;
1026 a1 = logf(-logf(1.0f-alpha1)) ;
1027 at = logf(-logf(1.0f-alphat)) ;
1028 xx = (x0 + (x1-x0)/(a1-a0)*(at-a0)) ; /* linear interp in a-space */
1029 } else if( fabsf(alphat-alpha0) <= fabsf(alphat-alpha1) ){
1030 xx = x0 ; /* alphat is closer to alpha0 */
1031 } else {
1032 xx = x1 ; /* alphat is closer to alpha1 */
1033 }
1034 return xx ;
1035 }
1036
1037 /*===========================================================================*/
1038 /* this is too complicated, and should be function-ized :( */
1039 /* which will never happen :) */
1040 /*===========================================================================*/
1041
main(int argc,char * argv[])1042 int main( int argc , char *argv[] )
1043 {
1044 int qpthr , ii,xx,yy,zz,ijk , dijk , qcase ;
1045 int ndilstep , ndilated[4] , ndilsum , ndiltot ;
1046 int count_targ100 , count_targ80, count_targ60 ;
1047 THD_3dim_dataset *qset=NULL ;
1048 float *qar=NULL , ***gthresh0=NULL , ***gthresh1=NULL, ***gthresh2=NULL ;
1049 char qpr[1024] ;
1050 MRI_IMAGE *cim0 =NULL , *cim1 =NULL , *cim2 =NULL ;
1051 MRI_IMARR **cimar0=NULL , **cimar1=NULL , **cimar2=NULL ;
1052 float ***car0 =NULL , ***car1 =NULL , ***car2 =NULL , **carHP ;
1053 float *farar=NULL ;
1054 int nfomkeep , nfar , itrac , ithresh , hp,ibr, ithresh_list[MAXITE] ;
1055 float tfrac=0.0006f, farperc=0.0f,farcut=0.0f, tfracold,farpercold,ttemp, farlast ;
1056 int ifarp ;
1057 float *gthrout , *gthrH ;
1058 int ntfp=0 , ntfp_all=0 ; float *tfs=NULL , *fps=NULL ;
1059
1060 /*----- help me if you can (I'm feeling down) -----*/
1061
1062 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
1063 printf("\n"
1064 "This program takes as input random field simulations\n"
1065 "(e.g., from 3dttest++) and does the ETAC processing to\n"
1066 "find cluster figure of merit (FOM) thresholds that are\n"
1067 "equitable (AKA balanced) across\n"
1068 " * voxel-wise p-values (-pthr option)\n"
1069 " * blurring cases (-ncase option)\n"
1070 " * H power values (-hpow option) -- probably not useful\n"
1071 "as well as being balanced across space to produce\n"
1072 "a False Positive Rate (FPR) that is approximately the\n"
1073 "same for each location and for each sub-case listed\n"
1074 "above. The usual goal is a global FPR of 5%%.\n"
1075 "\n"
1076 "* This program can be slow and consume a LOT of memory!\n"
1077 " (And I mean a BIG LOT, not a small lot.)\n"
1078 "\n"
1079 "* The output is a set of multi-threshold (*.mthresh.*.nii)\n"
1080 " files -- one for each of the -ncase inputs.\n"
1081 "\n"
1082 "* These files can be used via program 3dMultiThresh\n"
1083 " to produce an 'activation' mask.\n"
1084 "\n"
1085 "* 3dXClustSim is intended to be used from 3dttest++\n"
1086 " (via its '-ETAC' option) or some other script.\n"
1087 "\n"
1088 "* It is not intended to be run directly by any but the most\n"
1089 " knowledgeable and astute users. Which is why this help is so terse.\n"
1090 ) ;
1091
1092 printf("\n"
1093 "--------\n"
1094 "OPTIONS:\n"
1095 "--------\n"
1096 "\n"
1097 " -inset mask sdata ... {MANDATORY} [from 3dtoXdataset or 3dttest++]\n"
1098 " -insdat Data files are in the '.sdat' format.\n"
1099 "\n"
1100 " -NN 1 or 2 or 3 [-NN1 or -NN2 or -NN3 will work; default = 2]\n"
1101 " -sid 1 or 2 [-1sid or -2sid will work; default = 2]\n"
1102 " -hpow 0 1 2 [or some subset of these; default = 2]\n"
1103 "\n"
1104 " -ncase N lab1 .. labN [multiple processing cases; e.g., blurs]\n"
1105 " [default = 1 A]\n"
1106 " [example = 4 b04 b06 b08 b10]\n"
1107 "\n"
1108 " -pthr list of values [default = 0.0100 0.0056 0.0031 0.0018 0.0010]\n"
1109 " [equiv z1= 2.326 2.536 2.731 2.911 3.090 ]\n"
1110 " [equiv z2= 2.576 2.770 2.958 3.121 3.291 ]\n"
1111 "\n"
1112 " -FPR ff set global FPR goal to ff%%, where ff is an integer\n"
1113 " from 2 to 9 (inclusive). Default value is 5.\n"
1114 "\n"
1115 " -multiFPR compute results for multiple FPR goals (2%%, 3%%, ... 9%%)\n"
1116 "\n"
1117 " -minclust M don't allow clusters smaller than M voxels [default M=5]\n"
1118 "\n"
1119 " -local do the 'local' (voxelwise) ETAC computations\n"
1120 " -global do the 'global' (volumewise) ETAC computations\n"
1121 " -nolocal don't do the 'local'\n"
1122 " -noglobal don't do the 'global'\n"
1123 "\n"
1124 " -splitfrac F split simulations into pieces ( 0.2 < F < 0.8 )\n"
1125 "\n"
1126 " -prefix something useful\n"
1127 " -verb be more verbose\n"
1128 " -quiet silentium est aureum\n"
1129
1130 #if 0 /* disabled */
1131 " -FOMcount turn on FOMcount output\n"
1132 " -FARvox turn on FARvox output\n"
1133 #endif
1134
1135 "\n"
1136 "**-----------------------------------------------------------\n"
1137 "** Authored by Lamont Cranston, also known as ... The Shadow.\n"
1138 "**-----------------------------------------------------------\n"
1139 "\n"
1140 ) ;
1141 exit(0) ;
1142 }
1143
1144 (void)COX_clock_time() ; /* initialize the timer */
1145
1146 /* startup paperwork */
1147
1148 mainENTRY("3dXClustSim") ; machdep() ;
1149 AFNI_logger("3dXClustSim",argc,argv);
1150 PRINT_VERSION("3dXClustSim"); AUTHOR("Lamont Cranston") ;
1151
1152 #ifdef DEBUG_MEM
1153 enable_mcw_malloc() ;
1154 #ifdef USE_OMP
1155 omp_set_num_threads(1) ; /* will be slow */
1156 #endif
1157 #endif
1158
1159 /*----- load command line options -----*/
1160
1161 get_options(argc,argv) ;
1162
1163 #ifdef DEBUG_MEM
1164 if( verb < 2 ) pause_mcw_malloc() ;
1165 #endif
1166
1167 /*----- get the number of threads -----------------------------------*/
1168
1169 nthr = 1 ;
1170 #ifdef USE_OMP
1171 #pragma omp parallel /*------------ start parallel section ----------*/
1172 {
1173 #pragma omp master
1174 { nthr = omp_get_num_threads() ; }
1175 } /*---------- end parallel section ----------*/
1176 #endif
1177
1178 if( verb ){
1179 if( nthr > 1 )
1180 INFO_message("3dXClustSim: Using %d OpenMP threads",nthr) ;
1181 else
1182 INFO_message("3dXClustSim: Using 1 thread -- this will be slow :-(") ;
1183 }
1184
1185 /* setup some arrays in mri_threshX.c */
1186
1187 mri_multi_threshold_setup() ;
1188
1189 /*--- initialize the global cluster arrays ---*/
1190
1191 /*--- Xclustar_g[qcase][qpthr][iter]
1192 = array of clusters for case qcase,
1193 at p-value pthr[qpthr],
1194 for realization iter ---*/
1195
1196 Xclustar_g = (Xcluster_array ****)malloc(sizeof(Xcluster_array ***)*ncase) ;
1197 for( qcase=0 ; qcase < ncase ; qcase++ ){
1198 Xclustar_g[qcase] = (Xcluster_array ***)malloc(sizeof(Xcluster_array **)*npthr) ;
1199 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
1200 Xclustar_g[qcase][qpthr] = (Xcluster_array **)malloc(sizeof(Xcluster_array *)*niter_clust) ;
1201 }
1202 }
1203
1204 MEMORY_CHECK("startup") ;
1205
1206 /*===================================================================*/
1207 /*--- STEP 1a: loop over realizations to load up Xclustar_g[][][] ---*/
1208
1209 if( verb )
1210 INFO_message("STEP 1a: start %d-sided clustering with NN=%d and %d p-thresholds",nnsid,nnlev,npthr) ;
1211
1212 AFNI_OMP_START ; /*------------ start parallel section ----------*/
1213 #pragma omp parallel
1214 { DECLARE_ithr ;
1215 int iter , ipthr , icase , isim , isub ;
1216 MRI_IMAGE *fim , *tfim ; float *far , *tfar ;
1217
1218 /* initialize thread-specific workspace images */
1219
1220 #pragma omp critical
1221 { fim = mri_new_conforming( imtemplate , MRI_float ) ;
1222 tfim = mri_new_conforming( fim , MRI_float ) ; }
1223
1224 far = MRI_FLOAT_PTR(fim) ;
1225 tfar = MRI_FLOAT_PTR(tfim) ;
1226
1227 /* clusterize each volume at each p-value threshold */
1228
1229 /* #pragma omp for schedule(dynamic,500) */
1230 #pragma omp for
1231 for( isim=0; isim < nsim ; isim++ ){ /* loop over realizations */
1232 icase = isim / niter_clust ; iter = isim % niter_clust ;
1233 generate_image( far , iter+icase*niter ) ; /* get image to cluster */
1234 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){ /* loop over thresholds */
1235 /* creates Xclustar_g[icase][ipthr][iter] */
1236 gather_clusters( icase, ipthr , far, tfim, nnlev,nnsid, iter ) ;
1237 }
1238 if( verb > 2 && nthr == 1 && isim%1000 == 666 ) fprintf(stderr,".") ;
1239 if( 0 && verb > 2 ){
1240 #pragma omp critical
1241 { char mmm[32] ; sprintf(mmm,"isim=%d",isim) ; MEMORY_CHECK(mmm) ; }
1242 }
1243 }
1244
1245 /* delete thread-specific workspace for this processing block */
1246
1247 #pragma omp critical
1248 { mri_free(fim) ; mri_free(tfim) ; }
1249 }
1250 AFNI_OMP_END ; /*---------- end parallel section ----------*/
1251
1252 MEMORY_CHECK("after clustering") ;
1253
1254 /* don't need xinset's data for the nonce (maybe) */
1255 if( !do_global_etac ){
1256 if( verb > 1 ) ININFO_message("un-mapping input datasets") ;
1257 unmap_Xdataset(xinset) ;
1258 }
1259
1260 /*=====================================================================*/
1261 /* STEP 1b: count num clusters total, and merge them into one big list */
1262
1263 { int qter, qq,pp ; Xcluster_array *xcar ;
1264
1265 /* create the _tot global arrays for the merged cluster lists */
1266
1267 nclust_tot = (int **)malloc(sizeof(int *)*ncase) ;
1268 for( qcase=0 ; qcase < ncase ; qcase++ )
1269 nclust_tot[qcase] = (int *)malloc(sizeof(int)*npthr) ;
1270
1271 Xclust_tot = (Xcluster ****)malloc(sizeof(Xcluster_array ***)*ncase) ;
1272 for( qcase=0 ; qcase < ncase ; qcase++ )
1273 Xclust_tot[qcase] = (Xcluster ***)malloc(sizeof(Xcluster **)*npthr) ;
1274
1275 nclust_max = 0 ; /* keep track of the largest situation */
1276
1277 if( verb )
1278 ININFO_message("STEP 1b: merge cluster lists (%.1fs elapsed)",COX_clock_time()) ;
1279
1280 for( qcase=0 ; qcase < ncase ; qcase++ ){
1281 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
1282 nclust_tot[qcase][qpthr] = 0 ;
1283 for( qter=0 ; qter < niter_clust ; qter++ ){ /* count total clusters */
1284 xcar = Xclustar_g[qcase][qpthr][qter] ;
1285 if( xcar != NULL ) nclust_tot[qcase][qpthr] += xcar->nclu ;
1286 }
1287 /* max number of clusters across all p-value thresh, all cases */
1288 if( nclust_tot[qcase][qpthr] > nclust_max ) nclust_max = nclust_tot[qcase][qpthr] ;
1289 if( verb > 1 )
1290 ININFO_message(" %d total clusters :: Case %s pthr=%.5f",
1291 nclust_tot[qcase][qpthr], lcase[qcase] , pthr[qpthr] );
1292
1293 /* assemble all clusters at this case and p-value thresh into one array */
1294
1295 Xclust_tot[qcase][qpthr] = (Xcluster **)malloc(sizeof(Xcluster *)*nclust_tot[qcase][qpthr]) ;
1296
1297 for( pp=qter=0 ; qter < niter_clust ; qter++ ){ /* loop over realizations */
1298 xcar = Xclustar_g[qcase][qpthr][qter] ;
1299 if( xcar == NULL ) continue ;
1300 for( qq=0 ; qq < xcar->nclu ; qq++ ){ /* over clusts frm realization */
1301 Xclust_tot[qcase][qpthr][pp++] = xcar->xclu[qq] ; /* add to output array */
1302 }
1303 free(xcar->xclu) ; free(xcar) ;
1304 }
1305 free(Xclustar_g[qcase][qpthr]) ;
1306 }
1307 free(Xclustar_g[qcase]) ;
1308 }
1309 free(Xclustar_g) ; Xclustar_g = NULL ; /* not needed no more */
1310 }
1311
1312 MEMORY_CHECK("after counting") ;
1313
1314 /*=========================================================================*/
1315 /*--- STEP 1c: find the global distributions and min thresholds -----------*/
1316
1317 #if 0
1318 # define GTHRESH_FAC 0.799999f /* factor for method 1 */
1319 # define GTHRESH_THA 0.008888f /* how far into clust table: method 1 (per %) */
1320 # define GTHRESH_THB 0.013579f /* how far into clust table: method 2 (per %) */
1321 #else
1322 # define GTHRESH_FAC 0.135799f /* factor for method 1 */
1323 # define GTHRESH_THA 0.011111f /* how far into clust table: method 1 (per %) */
1324 # define GTHRESH_THB 0.044444f /* how far into clust table: method 2 (per %) */
1325 #endif
1326
1327 /*------------------------------------------------------------------------*/
1328 if( do_local_etac ){ /*----- not needed for global ETAC -----*/
1329 int nfom,jj,nfff; Xcluster **xcc;
1330 float a0,a1,f0,f1,fta,ftb , fmax , fg ;
1331 float *fomg0, *fomg1, *fomg2 ;
1332
1333 fomg0 = calloc(sizeof(float),nclust_max) ; /* workspaces */
1334 fomg1 = calloc(sizeof(float),nclust_max) ; /* for this */
1335 fomg2 = calloc(sizeof(float),nclust_max) ; /* section */
1336
1337 if( verb )
1338 ININFO_message("STEP 1c: compute minimum thresholds (%.1fs)",COX_clock_time()) ;
1339
1340 /* arrays for global threshold:
1341 [ifarp=FPR goal][qcase=blurring][ipthr=p-value] */
1342
1343 gthresh0 = (float ***)malloc(sizeof(float **)*numfarp) ;
1344 gthresh1 = (float ***)malloc(sizeof(float **)*numfarp) ;
1345 gthresh2 = (float ***)malloc(sizeof(float **)*numfarp) ;
1346 for( ifarp=0 ; ifarp < numfarp ; ifarp++ ){
1347 gthresh0[ifarp] = (float **)malloc(sizeof(float *)*ncase) ;
1348 gthresh1[ifarp] = (float **)malloc(sizeof(float *)*ncase) ;
1349 gthresh2[ifarp] = (float **)malloc(sizeof(float *)*ncase) ;
1350 for( qcase=0 ; qcase < ncase ; qcase++ ){
1351 gthresh0[ifarp][qcase] = (float *)calloc(sizeof(float),npthr) ; /* saved */
1352 gthresh1[ifarp][qcase] = (float *)calloc(sizeof(float),npthr) ; /* global */
1353 gthresh2[ifarp][qcase] = (float *)calloc(sizeof(float),npthr) ; /* thresholds */
1354 }
1355 }
1356
1357 for( qcase=0 ; qcase < ncase ; qcase++ ){
1358 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
1359 /* load all FOM values from all clusters for this case/pthr combo */
1360 xcc = Xclust_tot[qcase][qpthr] ;
1361 for( nfom=ii=0 ; ii < nclust_tot[qcase][qpthr] ; ii++ ){
1362 if( xcc[ii] != NULL ){
1363 fomg0[nfom] = xcc[ii]->fomh[0] ;
1364 fomg1[nfom] = xcc[ii]->fomh[1] ;
1365 fomg2[nfom] = xcc[ii]->fomh[2] ; nfom++ ;
1366 }
1367 }
1368 if( nfom < 50 ) continue ; /* should never happen */
1369 nfff = niter_clust ;
1370 if( nfff > nfom ) nfff = nfom ; /* very very unlikely */
1371
1372 fmax = AFNI_numenv("AFNI_XCLUSTSIM_FMAX") ;
1373 if( fmax <= 0.01f || fmax > 1.0f ) fmax = 0.888f ;
1374
1375 /* global threshold computed from tail of FOM distribution */
1376
1377 qsort_float_rev( nfom, fomg0 ) ; /* hpow=0 */
1378 qsort_float_rev( nfom, fomg1 ) ; /* hpow=1 */
1379 qsort_float_rev( nfom, fomg2 ) ; /* hpow=2 */
1380
1381 for( ifarp=0 ; ifarp < numfarp ; ifarp++ ){
1382 fg = farplist[ifarp] ; /* FPR goal in % */
1383
1384 jj = (int)rintf(GTHRESH_THA*nfff*fg) ; /* method 1 */
1385 if( jj >= nfom ) jj = nfom-1 ; /* should be impossible */
1386 fta = GTHRESH_FAC*fomg0[jj] ;
1387 jj = (int)rintf(GTHRESH_THB*nfff*fg) ; /* method 2 */
1388 if( jj >= nfom ) jj = nfom-1 ; /* should be impossible */
1389 ftb = fomg0[jj] ;
1390 gthresh0[ifarp][qcase][qpthr] = (int)(fmax*MAX(fta,ftb)+(1.0f-fmax)*MIN(fta,ftb)) ;
1391 if( verb > 1 && do_hpow0 ){
1392 ININFO_message(" min threshold %.1f :: Case %s pthr=%.5f h=0 fgoal=%.1f%%",
1393 gthresh0[ifarp][qcase][qpthr],lcase[qcase],pthr[qpthr],fg) ;
1394 }
1395
1396 jj = (int)rintf(GTHRESH_THA*nfff*fg) ;
1397 if( jj >= nfom ) jj = nfom-1 ; /* should be impossible */
1398 fta = GTHRESH_FAC*fomg1[jj] ;
1399 jj = (int)rintf(GTHRESH_THB*nfff*fg) ;
1400 if( jj >= nfom ) jj = nfom-1 ; /* should be impossible */
1401 ftb = fomg1[jj] ;
1402 gthresh1[ifarp][qcase][qpthr] = (fmax*MAX(fta,ftb)+(1.0f-fmax)*MIN(fta,ftb)) ;
1403 if( verb > 1 && do_hpow1 ){
1404 ININFO_message(" min threshold %.1f :: Case %s pthr=%.5f h=1 fgoal=%.1f%%",
1405 gthresh1[ifarp][qcase][qpthr],lcase[qcase],pthr[qpthr],fg) ;
1406 }
1407
1408 jj = (int)rintf(GTHRESH_THA*nfff*fg) ;
1409 if( jj >= nfom ) jj = nfom-1 ; /* should be impossible */
1410 fta = GTHRESH_FAC*fomg2[jj] ;
1411 jj = (int)rintf(GTHRESH_THB*nfff*fg) ;
1412 if( jj >= nfom ) jj = nfom-1 ; /* should be impossible */
1413 ftb = fomg2[jj] ;
1414 gthresh2[ifarp][qcase][qpthr] = (fmax*MAX(fta,ftb)+(1.0f-fmax)*MIN(fta,ftb)) ;
1415 if( verb > 1 && do_hpow2 ){
1416 ININFO_message(" min threshold %.1f :: Case %s pthr=%.5f h=2 fgoal=%.1f%%",
1417 gthresh2[ifarp][qcase][qpthr],lcase[qcase],pthr[qpthr],fg) ;
1418 }
1419 } /* end of loop over FPR goals (farp) */
1420 } /* end loop over thresholds */
1421 } /* end of loop over cases */
1422
1423 free(fomg0) ; free(fomg1) ; free(fomg2) ;
1424 MEMORY_CHECK("after globalizing") ;
1425 }
1426 /*------------------------------------------------------------------------*/
1427
1428 /*==========================================================================*/
1429 /*--- STEP 1d: do ETAC globally [Sep 2018] ---------------------------------*/
1430 /*:::::::::::: much of this code is a rehash of the voxelwise STEP 4, infra */
1431 /*:::::::::::: (which, ironically, is no longer default or recommended) */
1432
1433 if( do_global_etac ){
1434 float ***fomglob0, ***fomglob1, ***fomglob2 ;
1435 Xcluster **xcc ; int nfom,nfom_indx , **nfomglob ;
1436 float **fthar0 , **fthar1 , **fthar2 ;
1437 float *fgar0=NULL,*fgar1=NULL,*fgar2=NULL ,
1438 *tpar0=NULL,*tpar1=NULL,*tpar2=NULL ; /* for use_indx [08 Mar 2019] */
1439 int *indx=NULL ;
1440
1441 if( AFNI_noenv("AFNI_XCLUSTSIM_SORTER") ) use_indx = 0 ; /* 10 Mar 2019 */
1442
1443 if( verb )
1444 ININFO_message("STEP 1d: compute global ETAC thresholds (%.1fs)",COX_clock_time()) ;
1445
1446 /*--- create vectors to hold all FOMs ---*/
1447
1448 nfomkeep = (int)(TOPFRAC*niter_clust) ; /* max number of FOMs to keep */
1449
1450 fomglob0 = (float ***)malloc(sizeof(float **)*ncase) ;
1451 fomglob1 = (float ***)malloc(sizeof(float **)*ncase) ;
1452 fomglob2 = (float ***)malloc(sizeof(float **)*ncase) ;
1453 nfomglob = (int **) malloc(sizeof(int *) *ncase) ;
1454 fthar0 = (float **) malloc(sizeof(float *) *ncase) ;
1455 fthar1 = (float **) malloc(sizeof(float *) *ncase) ;
1456 fthar2 = (float **) malloc(sizeof(float *) *ncase) ;
1457 if( use_indx ){ /* temp arrays for sorting */
1458 indx = (int *) malloc(sizeof(int) *nclust_max) ;
1459 tpar0 = (float *) malloc(sizeof(float) *nclust_max) ;
1460 tpar1 = (float *) malloc(sizeof(float) *nclust_max) ;
1461 tpar2 = (float *) malloc(sizeof(float) *nclust_max) ;
1462 }
1463
1464 /* loops over diffusion cases, then over p-thresholds */
1465 /* extract arrays of FOM values for each case, for each p */
1466
1467 for( qcase=0 ; qcase < ncase ; qcase++ ){
1468 fomglob0[qcase] = (float **)malloc(sizeof(float *)*npthr) ;
1469 fomglob1[qcase] = (float **)malloc(sizeof(float *)*npthr) ;
1470 fomglob2[qcase] = (float **)malloc(sizeof(float *)*npthr) ;
1471 nfomglob[qcase] = (int *) malloc(sizeof(int) *npthr) ;
1472 fthar0[qcase] = (float *) malloc(sizeof(float) *npthr) ;
1473 fthar1[qcase] = (float *) malloc(sizeof(float) *npthr) ;
1474 fthar2[qcase] = (float *) malloc(sizeof(float) *npthr) ;
1475 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
1476 fgar0 = fomglob0[qcase][qpthr] = (float *)malloc(sizeof(float)*nclust_max);
1477 fgar1 = fomglob1[qcase][qpthr] = (float *)malloc(sizeof(float)*nclust_max);
1478 fgar2 = fomglob2[qcase][qpthr] = (float *)malloc(sizeof(float)*nclust_max);
1479 xcc = Xclust_tot[qcase][qpthr] ;
1480 for( nfom=ii=0 ; ii < nclust_tot[qcase][qpthr] ; ii++ ){
1481 if( xcc[ii] != NULL ){
1482 if( use_indx ) indx[nfom] = xcc[ii]->indx ; /* for sorting below */
1483 fgar0[nfom] = xcc[ii]->fomh[0] ;
1484 fgar1[nfom] = xcc[ii]->fomh[1] ;
1485 fgar2[nfom] = xcc[ii]->fomh[2] ; nfom++ ;
1486 }
1487 }
1488 #if 0
1489 ININFO_message("-- found %d FOMs for qcase=%d qpthr=%d out of %d clusters",nfom,qcase,qpthr,nclust_tot[qcase][qpthr]) ;
1490 #endif
1491
1492 /*- Order the results, biggest first, for cluster-threshold finding -*/
1493
1494 if( use_indx ){ /* eliminate multiple results from same iteration index */
1495
1496 /* step A = sort on indx to bring all results from same
1497 iteration together, if they aren't already. */
1498
1499 AAmemcpy(tpar0,fgar0,sizeof(float)*nfom) ; /* temp copies of */
1500 AAmemcpy(tpar1,fgar1,sizeof(float)*nfom) ; /* FOM values */
1501 AAmemcpy(tpar2,fgar2,sizeof(float)*nfom) ;
1502
1503 qsort_IFFF( nfom , indx , tpar0,tpar1,tpar2 ) ; /* sort on indx */
1504
1505 /* step B = for each separate value in indx, find the largest FOM,
1506 save into the original arrays, return number vals found */
1507
1508 nfom_indx = ixmax_IFFF( nfom , indx , tpar0,tpar1,tpar2 ,
1509 fgar0,fgar1,fgar2 ) ;
1510 if( verb > 1 )
1511 INFO_message(" ETAC sort @ p=%f: nfom before = %d after = %d",
1512 pthr[qpthr],nfom,nfom_indx) ;
1513 nfom = nfom_indx ;
1514 }
1515
1516 /* actual sorting of extracted FOM arrays now */
1517
1518 qsort_float_rev( nfom , fgar0 ) ; /* sort, largest first */
1519 qsort_float_rev( nfom , fgar1 ) ;
1520 qsort_float_rev( nfom , fgar2 ) ;
1521
1522 if( nfom > nfomkeep ){ /* trim the storage for frugality */
1523 fomglob0[qcase][qpthr] = (float *)realloc( fomglob0[qcase][qpthr],
1524 sizeof(float)*nfomkeep ) ;
1525 fomglob1[qcase][qpthr] = (float *)realloc( fomglob1[qcase][qpthr],
1526 sizeof(float)*nfomkeep ) ;
1527 fomglob2[qcase][qpthr] = (float *)realloc( fomglob2[qcase][qpthr],
1528 sizeof(float)*nfomkeep ) ;
1529 }
1530 nfomglob[qcase][qpthr] = MIN(nfom,nfomkeep) ; /* how many we have */
1531 if( verb > 1 )
1532 ININFO_message(" kept %d FOMs for qcase=%d qpthr=%d",
1533 nfomglob[qcase][qpthr],qcase,qpthr) ;
1534 } /* end of loop over pthr values */
1535 } /* end of loop over blurring cases */
1536
1537 /*--- now run the multithreshold simulations for global ETAC ---*/
1538 /*--- that is to say, this is where the results are computed ---*/
1539
1540 nit33 = 1.0f/(niter_clust+0.333f) ;
1541
1542 ntfp = 0 ; ntfp_all = 128 ;
1543 tfs = (float *)malloc(sizeof(float)*ntfp_all) ;
1544 fps = (float *)malloc(sizeof(float)*ntfp_all) ;
1545
1546 farlast = farplist[0] ;
1547 for( ifarp=0 ; ifarp < numfarp ; ifarp++ ){ /* loop over FAR goals */
1548
1549 farp_goal = farplist[ifarp] ;
1550
1551 if( verb )
1552 INFO_message("STEP 1d.%s: adjusting per-voxel FOM thresholds to reach FPR=%.2f%% (%.1fs)",
1553 abcd[ifarp] , farp_goal , COX_clock_time() ) ;
1554
1555 /* tfrac = FOM count fractional threshold;
1556 will be adjusted to find the farp_goal FPR goal */
1557
1558 if( ifarp == 0 ){ /* first time thru */
1559 tfrac = (4.0f+farp_goal)*0.000777f ;
1560 } else if( ntfp < 2 ){ /* if only 1 earlier calculation */
1561 tfrac *= 1.0777f * farp_goal / farlast ; /* adjust previous result */
1562 } else {
1563 int jj, jd ; float adf, mdf ; /* later: use closest in history */
1564 mdf = fabsf(fps[0]-farp_goal) ; jd = 0 ; /* to adjust starting point */
1565 for( jj=1 ; jj < ntfp ; jj++ ){
1566 adf = fabsf(fps[jj]-farp_goal) ;
1567 if( adf < mdf ){ mdf = adf ; jd = jj ; }
1568 }
1569 tfrac = 1.0111f * tfs[jd] * farp_goal / fps[jd] ;
1570 }
1571
1572 itrac = 0 ;
1573 farpercold = farperc = 0.0f ; tfracold = tfrac ;
1574
1575 /*--- Loop back to here with adjusted tfrac ---*/
1576
1577 GARP_LOOPBACK:
1578 {
1579 float min_tfrac,a0,a1 ; int nedge,nmin,ntest,npt,jthresh ;
1580 min_tfrac = 6.0f / niter_clust ; if( min_tfrac > 0.0001f ) min_tfrac = 0.0001f ;
1581
1582 itrac++ ; /* number of iterations */
1583 nfar = 0 ; /* total FAR count */
1584 nedge = nmin = ntest = 0 ; /* number of edge cases */
1585
1586 /* we take ithresh-th largest FOM for each case as the threshold */
1587 #if 0
1588 ithresh = (int)(tfrac*niter_clust) ; /* FOM count threshold */
1589 #else
1590 ithresh = (int)rintf(tfrac*(niter_clust-0.666f)+0.333f) ;
1591 #endif
1592
1593 /* Check if trying to re-litigate previous case [Cinco de Mayo 2017] */
1594
1595 ithresh_list[itrac-1] = ithresh ;
1596 if( itrac > 3 &&
1597 (ithresh_list[itrac-2] == ithresh || ithresh_list[itrac-3] == ithresh) ){
1598 ININFO_message(" #%d: would re-iterate at %g ==> %d ; breaking out",
1599 itrac,tfrac,ithresh ) ;
1600 goto GARP_BREAKOUT ;
1601 }
1602
1603 if( verb )
1604 ININFO_message(" #%d: Testing global thresholds at %g ==> %d",itrac,tfrac,ithresh) ;
1605
1606 /* compute the global thresholds for each case */
1607
1608 for( qcase=0 ; qcase < ncase ; qcase++ ){
1609 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
1610 npt = nfomglob[qcase][qpthr] ; /* how many FOM values here */
1611 jthresh = ithresh ; /* default index of FOM thresh */
1612
1613 if( jthresh+1 >= npt ){ /* edge case: too few FOM values */
1614 jthresh = npt-2 ; nedge++ ; /* should not happen */
1615 }
1616 a0 = ((float)jthresh+0.666f)*nit33 ;
1617 a1 = a0 + nit33 ;
1618 fthar0[qcase][qpthr] = inverse_interp_extreme(
1619 a0,a1,tfrac,
1620 fomglob0[qcase][qpthr][jthresh],
1621 fomglob0[qcase][qpthr][jthresh+1] ) ;
1622 fthar1[qcase][qpthr] = inverse_interp_extreme(
1623 a0,a1,tfrac,
1624 fomglob1[qcase][qpthr][jthresh],
1625 fomglob1[qcase][qpthr][jthresh+1] ) ;
1626 fthar2[qcase][qpthr] = inverse_interp_extreme(
1627 a0,a1,tfrac,
1628 fomglob2[qcase][qpthr][jthresh],
1629 fomglob2[qcase][qpthr][jthresh+1] ) ;
1630 }}
1631
1632 /*-------------- now loop over realizations and threshold them --------*/
1633 AFNI_OMP_START ; /*------------ start parallel section ------------------*/
1634 #pragma omp parallel
1635 { DECLARE_ithr ;
1636 MRI_IMAGE *fim , *tfim ; float *far , *tfar ;
1637 int iter,npt,iv,jthresh , ipthr,icase , qq,hh ;
1638 float a0,a1,f0,f1,ft ;
1639
1640 /* workspace image array for each thread (from realizations) */
1641
1642 #pragma omp critical
1643 { fim = mri_new_conforming( imtemplate , MRI_float ) ;
1644 far = MRI_FLOAT_PTR(fim) ; }
1645
1646 /* use the thresholds computed above to multi-threshold
1647 each realization, and create count of total FA and in each voxel */
1648
1649 /* #pragma omp for schedule(dynamic,333) */
1650 #pragma omp for /* parallelized */
1651 for( iter=niter_tbot ; iter < niter_ttop ; iter++ ){ /* loop over iterations */
1652 for( icase=0 ; icase < ncase ; icase++ ){ /* over cases */
1653 generate_image( far , iter+icase*niter ) ; /* get image to test */
1654
1655 /* multi-threshold this image with global FOM thresholds */
1656
1657 tfim = mri_multi_threshold_Xcluster_fomth( fim, npthr,zthr_used,
1658 nnsid,nnlev,
1659 fthar0[icase] ,
1660 fthar1[icase] ,
1661 fthar2[icase] ,
1662 0, NULL , NULL ) ;
1663 if( tfim != NULL ){ /* nothing found ==> NULL was returned */
1664 #pragma omp critical
1665 { mri_free(tfim) ; nfar++ ;} /* nfar = count of total FA */
1666 break ; /* no need to continue thru cases (got a FA) */
1667 }
1668 } /* end of loop over cases */
1669 } /* end of loop over iterations (parallelized) */
1670
1671 #pragma omp critical
1672 { mri_free(fim) ; } /* tossing the trash for this thread */
1673 }
1674 AFNI_OMP_END ; /*------------ end parallel section ------------*/
1675
1676 /* compute the FAR percentage at this tfrac */
1677
1678 farpercold = farperc ; /* save last time result */
1679 farperc = (100.0f*nfar)/(float)niter_test ; /* what we got this time */
1680 farlast = farperc ; /* save for next FPR goal */
1681
1682 /* save results for later re-use */
1683 if( ntfp == ntfp_all ){
1684 ntfp_all += 128 ;
1685 tfs = (float *)realloc(tfs,sizeof(float)*ntfp_all) ;
1686 fps = (float *)realloc(tfs,sizeof(float)*ntfp_all) ;
1687 }
1688 tfs[ntfp] = tfrac ; fps[ntfp] = farperc ; ntfp++ ;
1689
1690 /* farcut = precision desired for our FPR goal */
1691 farcut = 0.222f ;
1692 if( itrac > 2 ) farcut += (itrac-2)*0.0321f ;
1693
1694 if( verb )
1695 ININFO_message(" global FPR=%.2f%% (%.1fs)", farperc,COX_clock_time() ) ;
1696 MEMORY_CHECK(" ") ;
1697
1698 /* if no substantial progress, quit */
1699
1700 if( itrac > 2 && fabsf(farperc-farpercold) < 0.0222f ){
1701 ININFO_message(" ((Progress too slow - breaking out **))") ;
1702 goto GARP_BREAKOUT ;
1703 }
1704
1705 /* try another tfrac to get closer to our goal? */
1706
1707 if( itrac < MAXITE && fabsf(farperc-FG_GOAL) > farcut ){
1708 float fff , dtt ;
1709 if( itrac == 1 || (farperc-FG_GOAL)*(farpercold-FG_GOAL) > 0.1f ){ /* scale */
1710 fff = FG_GOAL/farperc ;
1711 if( fff > 2.222f ) fff = 2.222f ; else if( fff < 0.450f ) fff = 0.450f ;
1712 dtt = (fff-1.0f)*tfrac ; /* tfrac step */
1713 if( itrac > 2 ) dtt *= (0.8666f+0.2222f*itrac) ; /* accelerate it */
1714 ttemp = tfrac ; tfrac += dtt ;
1715 } else { /* linear inverse interpolate */
1716 fff = (farperc-farpercold)/(tfrac-tfracold) ;
1717 ttemp = tfrac ; tfrac = tfracold + (FG_GOAL-farpercold)/fff ;
1718 }
1719 #define TFTOP (0.666f*TOPFRAC)
1720 tfracold = ttemp ;
1721 if( tfrac < min_tfrac ) tfrac = min_tfrac ;
1722 else if( tfrac > TFTOP ) tfrac = TFTOP ;
1723 goto GARP_LOOPBACK ;
1724 }
1725
1726 GARP_BREAKOUT: ; /*nada*/
1727 } /* end of iterations to find the ideal farperc */
1728
1729 /*---::: Write results out, then this farp_goal is done :::---*/
1730
1731 { char *qpref , *cpt , *qax[2] ;
1732 NI_element *qel ;
1733 float *qar ; int nqq , qdim[2] , iqq ; float gval ;
1734
1735 /* constants for the NIML element to be created for each case */
1736
1737 nqq = nhpow * npthr ; /* number of thresholds */
1738 qar = (float *)malloc(sizeof(float)*nqq) ; /* thresholds */
1739
1740 qdim[0] = npthr ; qdim[1] = nhpow; /* order of data */
1741 qax[0] = "zthr"; qax[1] = "FOM"; /* labels for fun */
1742
1743 qpref = strdup(prefix) ; /* mangle the prefix */
1744 cpt = strstr(qpref,".nii") ; if( cpt != NULL ) *cpt = '\0' ;
1745
1746 for( qcase=0 ; qcase < ncase ; qcase++ ){ /* loop over cases */
1747
1748 /* create a NIML element to hold the thresholds */
1749
1750 qel = NI_new_data_element( "global_ETAC_thresholds" , nqq ) ;
1751 NI_set_dimen( qel , 2 , qdim ) ;
1752 NI_set_axes ( qel , qax ) ;
1753
1754 /* store cluster-FOM thresholds we computed above */
1755
1756 for( iqq=qpthr=0 ; qpthr < npthr ; qpthr++ ){
1757 if( do_hpow0 ) qar[iqq++] = fthar0[qcase][qpthr] ;
1758 if( do_hpow1 ) qar[iqq++] = fthar1[qcase][qpthr] ;
1759 if( do_hpow2 ) qar[iqq++] = fthar2[qcase][qpthr] ;
1760
1761 if( do_local_etac ){ /* Possibly reduce local min thresholds [23 Jan 2019] */
1762 gval = 0.222f * fthar0[qcase][qpthr] ;
1763 if( do_hpow0 && gval < gthresh0[ifarp][qcase][qpthr] ){
1764 ININFO_message("deltamin threshold %.1f (%.1f):: Case %s pthr=%.5f h=0 fgoal=%.1f%%",
1765 gval,gthresh0[ifarp][qcase][qpthr],lcase[qcase],pthr[qpthr],farp_goal) ;
1766 gthresh0[ifarp][qcase][qpthr] = gval ;
1767 }
1768 gval = 0.222f * fthar1[qcase][qpthr] ;
1769 if( do_hpow1 && gval < gthresh1[ifarp][qcase][qpthr] ){
1770 ININFO_message("deltamin threshold %.1f (%.1f):: Case %s pthr=%.5f h=1 fgoal=%.1f%%",
1771 gval,gthresh1[ifarp][qcase][qpthr],lcase[qcase],pthr[qpthr],farp_goal) ;
1772 gthresh1[ifarp][qcase][qpthr] = gval ;
1773 }
1774 gval = 0.222f * fthar2[qcase][qpthr] ;
1775 if( do_hpow2 && gval < gthresh2[ifarp][qcase][qpthr] ){
1776 ININFO_message("deltamin threshold %.1f (%.1f):: Case %s pthr=%.5f h=2 fgoal=%.1f%%",
1777 gval,gthresh2[ifarp][qcase][qpthr],lcase[qcase],pthr[qpthr],farp_goal) ;
1778 gthresh2[ifarp][qcase][qpthr] = gval ;
1779 }
1780 }
1781
1782 }
1783
1784 /* store them into the element */
1785
1786 NI_add_column ( qel , NI_FLOAT , qar ) ;
1787 NI_set_column_label( qel , 0 , "thresh" ) ;
1788
1789 /* add some attributes for later use in 3dMultiThresh */
1790
1791 /* simple constants */
1792
1793 sprintf(qpr,"%d",nnlev) ; NI_set_attribute(qel,"NNlev",qpr) ;
1794 sprintf(qpr,"%d",nnsid) ; NI_set_attribute(qel,"NNsid",qpr) ;
1795 sprintf(qpr,"%d",min_clust) ; NI_set_attribute(qel,"min_clust",qpr) ;
1796 NI_set_attribute( qel , "case" , lcase[qcase] ) ;
1797
1798 /* list of FOMs used */
1799
1800 qpr[0] = '\0' ; iqq = 0 ;
1801 if( do_hpow0 ){ strcat(qpr,"hpow0") ; iqq++ ; }
1802 if( do_hpow1 ){
1803 if( iqq ) strcat(qpr,",") ;
1804 strcat(qpr,"hpow1") ; iqq++ ;
1805 }
1806 if( do_hpow2 ){
1807 if( iqq ) strcat(qpr,",") ;
1808 strcat(qpr,"hpow2") ; iqq++ ;
1809 }
1810 NI_set_attribute( qel , "FOM_list" , qpr ) ;
1811
1812 /* list of voxelwise thresholds used */
1813
1814 sprintf(qpr,"%.5f",zthr_used[0]) ;
1815 for( qpthr=1 ; qpthr < npthr ; qpthr++ )
1816 sprintf(qpr+strlen(qpr),",%.5f",zthr_used[qpthr]) ;
1817 NI_set_attribute( qel , "zthr_list" , qpr ) ;
1818
1819 /* FAR goal (not really needed) */
1820
1821 sprintf(qpr,"%.2f",farp_goal) ;
1822 NI_set_attribute( qel , "FAR_goal" , qpr ) ;
1823
1824 /* invent a filename for the output */
1825
1826 sprintf( qpr ,"globalETAC.mthresh.%s.%s.%02dperc.niml" ,
1827 qpref , lcase[qcase] , (int)rintf(farp_goal) ) ;
1828
1829 /* and write it out */
1830
1831 NI_write_element_tofile( qpr , qel , NI_TEXT_MODE ) ;
1832 NI_free_element(qel) ;
1833
1834 if( verb )
1835 ININFO_message("global ETAC thresholds written to %s",qpr) ;
1836
1837 } /* end of loop over cases */
1838
1839 free(qar) ; free(qpref) ;
1840
1841 } /* end of writing out NIML elements with global thresholds */
1842
1843 } /*--- end of loop over FAR goals ---*/
1844
1845 /*::: it's NOT the end of the world (yet) :::*/
1846
1847 /*--- delete vectors ---*/
1848
1849 if( use_indx ){
1850 free(indx) ; free(tpar0) ; free(tpar1) ; free(tpar2) ;
1851 }
1852
1853 for( qcase=0 ; qcase < ncase ; qcase++ ){
1854 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
1855 free(fomglob0[qcase][qpthr]) ;
1856 free(fomglob1[qcase][qpthr]) ;
1857 free(fomglob2[qcase][qpthr]) ;
1858 }
1859 free(fomglob0[qcase]); free(fomglob1[qcase]); free(fomglob2[qcase]);
1860 free(nfomglob[qcase]);
1861 free(fthar0[qcase]) ; free(fthar1[qcase]) ; free(fthar2[qcase]) ;
1862 }
1863 free(fomglob0); free(fomglob1); free(fomglob2); free(nfomglob);
1864 free(fthar0) ; free(fthar1) ; free(fthar2) ;
1865
1866 /*-------- can exit NOW NOW NOW unless local ETAC is on the menu --------*/
1867
1868 if( !do_local_etac ){
1869 INFO_message("=== 3dXClustSim ends: Elapsed time %.1fs",COX_clock_time()) ;
1870 exit(0) ;
1871 }
1872
1873 /* otherwise, can unmap simulations now, until needed later */
1874
1875 if( verb > 1 ) ININFO_message("un-mapping input datasets") ;
1876 unmap_Xdataset(xinset) ;
1877
1878 } /*--------------- end of global ETAC-ization ---------------*/
1879
1880 /*------------------------------------------------------------------------*/
1881 /*-------- At this time [March 2019] local ETAC should be regarded -------*/
1882 /*-------- as both highly experimental and uncertainly accurate :( -------*/
1883 /*------------------------------------------------------------------------*/
1884
1885 /***************************************************/
1886 /***** From here on is local (voxel-wise) ETAC *****/
1887 /***************************************************/
1888
1889 /*==========================================================================*/
1890 /*--- STEP 2: dilate the clusters ------------------------------------------*/
1891 /* The goal is to have at least so many "hits" at each voxel in the mask */
1892
1893 /*--- initialize the FOM vector array for each voxel ---*/
1894 /*--- this is the list of FOMs (e.g., cluster sizes)
1895 that were found at that voxel ---*/
1896
1897 fomvec0 = (Xvector **)malloc(sizeof(Xvector *)*mask_ngood) ;
1898 fomvec1 = (Xvector **)malloc(sizeof(Xvector *)*mask_ngood) ;
1899 fomvec2 = (Xvector **)malloc(sizeof(Xvector *)*mask_ngood) ;
1900 for( ii=0 ; ii < mask_ngood ; ii++ ){
1901 CREATE_Xvector(fomvec0[ii],DALL) ;
1902 fomvec0[ii]->ip = ipmask[ii] ; /* store location of this */
1903 fomvec0[ii]->jp = jpmask[ii] ; /* vector in the 3D grid */
1904 fomvec0[ii]->kp = kpmask[ii] ;
1905 fomvec0[ii]->ijk = ijkmask[ii] ;
1906 CREATE_Xvector(fomvec1[ii],DALL) ;
1907 fomvec1[ii]->ip = ipmask[ii] ;
1908 fomvec1[ii]->jp = jpmask[ii] ;
1909 fomvec1[ii]->kp = kpmask[ii] ;
1910 fomvec1[ii]->ijk = ijkmask[ii] ;
1911 CREATE_Xvector(fomvec2[ii],DALL) ;
1912 fomvec2[ii]->ip = ipmask[ii] ;
1913 fomvec2[ii]->jp = jpmask[ii] ;
1914 fomvec2[ii]->kp = kpmask[ii] ;
1915 fomvec2[ii]->ijk = ijkmask[ii] ;
1916 }
1917
1918 /* temp space to hold counts of FOM vec length for each voxel in a cluster */
1919 /* see get_Xcluster_nbcount() for usage */
1920
1921 ncount_g = (int *) malloc(sizeof(int) *nthr) ; /* dimens of fcount_g */
1922 fcount_g = (float **)malloc(sizeof(float *)*nthr) ; /* counts */
1923 for( ii=0 ; ii < nthr ; ii++ ){
1924 fcount_g[ii] = (float *)malloc(sizeof(float)*DALL) ;
1925 ncount_g[ii] = DALL ;
1926 }
1927
1928 /*--- thread-specific workspaces for clustering int dilate_Xcluster() ---*/
1929
1930 ndilg = (int * )malloc(sizeof(int )*nthr) ; /* size of dilg_ijk */
1931 dilg_ijk = (int **)malloc(sizeof(int *)*nthr) ; /* candidates for dilate */
1932 for( ii=0 ; ii < nthr ; ii++ ){
1933 ndilg [ii] = DALL ;
1934 dilg_ijk[ii] = (int *)malloc(sizeof(int)*DALL) ;
1935 }
1936
1937 /*----- loop over p-thresholds and dilate in voxel chunks -----*/
1938
1939 dijk = nxyz / nthr ; /* set chunk size */
1940 if( dijk > nxy ) dijk = nxy ;
1941 if( dijk > 4096 ) dijk = 4096 ;
1942
1943 /* target counts for voxel "hits" */
1944
1945 count_targ100 = (int)rintf(dilate_fac*niter_clust) ;
1946 if( count_targ100 > 333 ) /* why 333? It's reasonably sized */
1947 count_targ100 = (int)rintf(sqrtf(333.0f*count_targ100)) ;
1948 count_targ80 = (int)rintf(0.80f*count_targ100) ;
1949 count_targ60 = (int)rintf(0.60f*count_targ100) ;
1950
1951 if( verb )
1952 INFO_message("STEP 2: start cluster dilations (%.1fs)",COX_clock_time()) ;
1953
1954 for( qcase=0 ; qcase < ncase ; qcase++ ){ /* loop over cases */
1955 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){ /* loop over p-value thresh */
1956
1957 Xcluster **xcar = Xclust_tot[qcase][qpthr] ;
1958 int ncar = nclust_tot[qcase][qpthr] ;
1959
1960 #define NDILMAX 9 /* max number of dilation steps */
1961
1962 ndiltot = 0 ;
1963 for( ndilstep=0 ; ndilstep < NDILMAX ; ndilstep++ ){
1964 /* initialize counts of number of dilations for each NN type */
1965 ndilated[0] = ndilated[1] = ndilated[2] = ndilated[3] = 0 ;
1966
1967 #ifdef DEBUG_MEM
1968 if( verb > 1 ){
1969 char mmm[256] ;
1970 sprintf(mmm,"before dilation: qcase=%d qpthr=%d ndilstep=%d",qcase,qpthr,ndilstep) ;
1971 MEMORY_CHECK(mmm) ;
1972 }
1973 #endif
1974
1975 AFNI_OMP_START ; /*------------ start parallel section ----------*/
1976 #pragma omp parallel
1977 { DECLARE_ithr ;
1978 int iter, idil , ijkbot,ijktop , nclu=nclust_tot[qcase][qpthr] ;
1979
1980 /* reset the FOM vector to be empty, for each voxel */
1981
1982 #pragma omp for
1983 for( iter=0 ; iter < mask_ngood ; iter++ ){
1984 fomvec0[iter]->npt = 0 ;
1985 fomvec1[iter]->npt = 0 ;
1986 fomvec2[iter]->npt = 0 ;
1987 }
1988
1989 /* set the Xcluster index ranges */
1990
1991 #pragma omp for
1992 for( iter=0 ; iter < ncar ; iter++ ){
1993 set_Xcluster_ijkminmax( xcar[iter] ) ;
1994 }
1995
1996 /* load the FOM vectors (parallelized across 3D segments) */
1997
1998 /* #pragma omp for schedule(dynamic,1) */
1999 #pragma omp for
2000 for( ijkbot=0 ; ijkbot < nxyz ; ijkbot+=dijk ){
2001 ijktop = ijkbot + (dijk-1) ; if( ijktop >= nxyz ) ijktop = nxyz-1 ;
2002 process_clusters_to_Xvectors( ijkbot, ijktop , qcase,qpthr ) ;
2003 }
2004
2005 /* get median FOM counts in each cluster,
2006 then determine how to dilate (NN1, NN2, NN3)
2007 [parallelized across clusters = simulation iterations] */
2008
2009 /* #pragma omp for schedule(dynamic,666) */
2010 #pragma omp for
2011 for( iter=0 ; iter < nclu ; iter++ ){
2012 if( Xclust_tot[qcase][qpthr][iter] != NULL &&
2013 Xclust_tot[qcase][qpthr][iter]->nbcount < count_targ100 ){
2014 idil = get_Xcluster_nbcount( Xclust_tot[qcase][qpthr][iter] , ithr ) ;
2015 if( idil < count_targ100 ){ /* too few? */
2016 /* dilate: NN3 or NN2 or NN1 */
2017 idil = (idil < count_targ60) ? 3 :(idil < count_targ80) ? 2 : 1 ;
2018 dilate_Xcluster( Xclust_tot[qcase][qpthr][iter] , idil , ithr ) ;
2019 #pragma omp atomic
2020 ndilated[idil]++ ; /* count how many dilations of each type */
2021 }
2022 }
2023 } /* all clusters dilated now (or were OK) */
2024
2025 }
2026 AFNI_OMP_END ; /*---------- end parallel section ----------*/
2027
2028 /* number of significant dilations (NN2+NN3) */
2029
2030 ndilsum = ndilated[2] + ndilated[3] ;
2031 ndiltot += ndilated[1] + ndilated[2] + ndilated[3] ;
2032
2033 /* if not very many, then we are done with this p-value thresh */
2034
2035 if( ndilstep < NDILMAX-1 && ndilsum < niter_clust/50 ) break ;
2036 } /* end of loop over dilation steps */
2037 if( verb > 1 )
2038 ININFO_message(" %d dilation loops; %d total cluster dilations :: Case %s pthr=%.5f",
2039 ndilstep,ndiltot,lcase[qcase],pthr[qpthr]) ;
2040 } /* end of loop over p-value thresh cluster collection */
2041 } /* end of loop over cases */
2042
2043 /* free the dilation counting workspace for each thread */
2044
2045 for( ii=0 ; ii < nthr ; ii++ ) free(dilg_ijk[ii]) ;
2046 free(ndilg); free(dilg_ijk);
2047
2048 MEMORY_CHECK("after all dilations") ;
2049
2050 /*=======================================================*/
2051 /*--- STEP 3: create sorted and truncated FOM vectors ---*/
2052
2053 if( verb )
2054 INFO_message("STEP 3: re-loading & sorting FOM vectors after dilations (%.1fs)",COX_clock_time()) ;
2055
2056 /*--- initialize the final sorted FOM vector array for each voxel ---*/
2057
2058 fomsort0 = (Xvector ****)malloc(sizeof(Xvector ***)*ncase) ; /* hpow=0 */
2059 fomsort1 = (Xvector ****)malloc(sizeof(Xvector ***)*ncase) ; /* =1 */
2060 fomsort2 = (Xvector ****)malloc(sizeof(Xvector ***)*ncase) ; /* =2 */
2061 for( qcase=0 ; qcase < ncase ; qcase++ ){
2062 fomsort0[qcase] = (Xvector ***)malloc(sizeof(Xvector **)*npthr) ;
2063 fomsort1[qcase] = (Xvector ***)malloc(sizeof(Xvector **)*npthr) ;
2064 fomsort2[qcase] = (Xvector ***)malloc(sizeof(Xvector **)*npthr) ;
2065 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
2066 fomsort0[qcase][qpthr] = (Xvector **)malloc(sizeof(Xvector *)*mask_ngood) ;
2067 fomsort1[qcase][qpthr] = (Xvector **)malloc(sizeof(Xvector *)*mask_ngood) ;
2068 fomsort2[qcase][qpthr] = (Xvector **)malloc(sizeof(Xvector *)*mask_ngood) ;
2069 }
2070 }
2071
2072 nfomkeep = (int)(TOPFRAC*niter_clust) ; /* max number of FOMs to keep at 1 voxel */
2073
2074 for( qcase=0 ; qcase < ncase ; qcase++ ){ /* loop over cases */
2075 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){ /* loop over p-value thresh */
2076
2077 Xcluster **xcar = Xclust_tot[qcase][qpthr] ;
2078 int ncar = nclust_tot[qcase][qpthr] ;
2079
2080 #ifdef DEBUG_MEM
2081 if( verb > 1 ){
2082 char mmm[256] ;
2083 sprintf(mmm,"before FOMsort: qcase=%d qpthr=%d",qcase,qpthr) ;
2084 MEMORY_CHECK(mmm) ;
2085 }
2086 #else
2087 if( verb > 1 )
2088 ININFO_message(" start Case %s pthr=%.5f",lcase[qcase],pthr[qpthr]) ;
2089 #endif
2090
2091 AFNI_OMP_START ; /*------------ start parallel section ----------*/
2092 # pragma omp parallel
2093 { DECLARE_ithr ;
2094 int ijkbot,ijktop , iv,jj , npt , nclu=nclust_tot[qcase][qpthr] ;
2095 float a0,a1,f0,f1,ft ;
2096
2097 /* re-load FOM vectors for this threshold */
2098
2099 #pragma omp for
2100 for( iv=0 ; iv < mask_ngood ; iv++ ){
2101 fomvec0[iv]->npt = 0 ;
2102 fomvec1[iv]->npt = 0 ;
2103 fomvec2[iv]->npt = 0 ;
2104 }
2105
2106 /* set the Xcluster index ranges */
2107
2108 #pragma omp for
2109 for( iv=0 ; iv < ncar ; iv++ ){
2110 set_Xcluster_ijkminmax( xcar[iv] ) ;
2111 }
2112
2113 /* #pragma omp for schedule(dynamic,1) */
2114 #pragma omp for
2115 for( ijkbot=0 ; ijkbot < nxyz ; ijkbot+=dijk ){
2116 ijktop = ijkbot + (dijk-1) ;
2117 process_clusters_to_Xvectors( ijkbot, ijktop , qcase,qpthr ) ;
2118 }
2119
2120 /* delete Xclusters for this threshold,
2121 since we have the finalized FOM vectors */
2122
2123 #pragma omp master
2124 { for( jj=0 ; jj < nclu ; jj++ ){
2125 DESTROY_Xcluster(Xclust_tot[qcase][qpthr][jj]) ;
2126 }
2127 free(Xclust_tot[qcase][qpthr]) ; Xclust_tot[qcase][qpthr] = NULL ;
2128 }
2129 #pragma omp barrier
2130
2131 /* sort/truncate FOM vectors for this threshold */
2132
2133 /* vectors loaded for each pt in space;
2134 now sort them (biggest first) to determine equitable thresholds */
2135
2136 /* #pragma omp for schedule(dynamic,666) */
2137 #pragma omp for
2138 for( iv=0 ; iv < mask_ngood ; iv++ ){ /* loop over voxels in mask */
2139 if( fomvec0[iv]->npt < 3 ){ /* if it's way short */
2140 ADDTO_Xvector(fomvec0[iv],0.01f) ; /* put some padding in */
2141 ADDTO_Xvector(fomvec0[iv],0.01f) ;
2142 ADDTO_Xvector(fomvec0[iv],0.01f) ;
2143 ADDTO_Xvector(fomvec1[iv],0.01f) ;
2144 ADDTO_Xvector(fomvec1[iv],0.01f) ;
2145 ADDTO_Xvector(fomvec1[iv],0.01f) ;
2146 ADDTO_Xvector(fomvec2[iv],0.01f) ;
2147 ADDTO_Xvector(fomvec2[iv],0.01f) ;
2148 ADDTO_Xvector(fomvec2[iv],0.01f) ;
2149 }
2150
2151 /* sort into decreasing order (largest FOMs first) */
2152
2153 if( do_hpow0 ) qsort_float_rev( fomvec0[iv]->npt , fomvec0[iv]->far ) ;
2154 if( do_hpow1 ) qsort_float_rev( fomvec1[iv]->npt , fomvec1[iv]->far ) ;
2155 if( do_hpow2 ) qsort_float_rev( fomvec2[iv]->npt , fomvec2[iv]->far ) ;
2156
2157 /* create a new vector to get the keeper FOM counts */
2158
2159 if( do_hpow0 ){
2160 jj = MIN(fomvec0[iv]->npt,nfomkeep) ; /* how many to keep */
2161 CREATE_Xvector(fomsort0[qcase][qpthr][iv],jj) ;
2162 fomsort0[qcase][qpthr][iv]->ip = ipmask[iv] ; /* location */
2163 fomsort0[qcase][qpthr][iv]->jp = jpmask[iv] ;
2164 fomsort0[qcase][qpthr][iv]->kp = kpmask[iv] ;
2165 fomsort0[qcase][qpthr][iv]->ijk = ijkmask[iv] ;
2166 fomsort0[qcase][qpthr][iv]->npt = jj ; /* count */
2167 AAmemcpy(fomsort0[qcase][qpthr][iv]->far,fomvec0[iv]->far,sizeof(float)*jj) ;
2168 }
2169
2170 if( do_hpow1 ){
2171 jj = MIN(fomvec1[iv]->npt,nfomkeep) ; /* how many to keep */
2172 CREATE_Xvector(fomsort1[qcase][qpthr][iv],jj) ;
2173 fomsort1[qcase][qpthr][iv]->ip = ipmask[iv] ; /* location */
2174 fomsort1[qcase][qpthr][iv]->jp = jpmask[iv] ;
2175 fomsort1[qcase][qpthr][iv]->kp = kpmask[iv] ;
2176 fomsort1[qcase][qpthr][iv]->ijk = ijkmask[iv] ;
2177 fomsort1[qcase][qpthr][iv]->npt = jj ; /* count */
2178 AAmemcpy(fomsort1[qcase][qpthr][iv]->far,fomvec1[iv]->far,sizeof(float)*jj) ;
2179 }
2180
2181 if( do_hpow2 ){
2182 jj = MIN(fomvec2[iv]->npt,nfomkeep) ; /* how many to keep */
2183 CREATE_Xvector(fomsort2[qcase][qpthr][iv],jj) ;
2184 fomsort2[qcase][qpthr][iv]->ip = ipmask[iv] ; /* location */
2185 fomsort2[qcase][qpthr][iv]->jp = jpmask[iv] ;
2186 fomsort2[qcase][qpthr][iv]->kp = kpmask[iv] ;
2187 fomsort2[qcase][qpthr][iv]->ijk = ijkmask[iv] ;
2188 fomsort2[qcase][qpthr][iv]->npt = jj ; /* count */
2189 AAmemcpy(fomsort2[qcase][qpthr][iv]->far,fomvec2[iv]->far,sizeof(float)*jj) ;
2190 }
2191 } /* end of voxel loop */
2192 }
2193 AFNI_OMP_END ; /*---------- end parallel section ----------*/
2194
2195 /* save dataset of fomvec counts */
2196
2197 #ifdef ALLOW_EXTRAS
2198 if( do_FOMcount ){
2199 qset = EDIT_empty_copy(mask_dset) ;
2200 sprintf(qpr,".FOMc.%s.%d",lcase[qcase],qpthr) ;
2201 EDIT_dset_items( qset ,
2202 ADN_prefix , modify_afni_prefix(prefix,NULL,qpr) ,
2203 ADN_nvals , 1 ,
2204 ADN_none ) ;
2205 EDIT_substitute_brick( qset , 0 , MRI_float , NULL ) ;
2206 qar = DSET_ARRAY(qset,0) ;
2207 for( ii=0 ; ii < mask_ngood ; ii++ ){
2208 ijk = ijkmask[ii] ;
2209 qar[ijk] = (float)fomvec0[ii]->npt ;
2210 fomvec0[ii]->npt = 0 ;
2211 }
2212 EDIT_BRICK_LABEL(qset,0,"FOMcount") ;
2213 tross_Make_History( "3dXClustSim" , argc,argv , qset ) ;
2214 DSET_write(qset); WROTE_DSETI(qset);
2215 DSET_delete(qset); qset = NULL; qar = NULL;
2216 }
2217 #endif
2218
2219 } /* end of loop over qpthr */
2220 } /* end of loop over cases */
2221
2222 /* destroy the unsorted full FOM vectors */
2223
2224 for( ii=0 ; ii < mask_ngood ; ii++ ){ /* don't need no more */
2225 DESTROY_Xvector(fomvec0[ii]) ;
2226 DESTROY_Xvector(fomvec1[ii]) ;
2227 DESTROY_Xvector(fomvec2[ii]) ;
2228 }
2229 free(fomvec0) ; fomvec0 = NULL ;
2230 free(fomvec1) ; fomvec1 = NULL ;
2231 free(fomvec2) ; fomvec2 = NULL ;
2232
2233 MEMORY_CHECK("after all FOMsort") ;
2234
2235 /*=============================================================*/
2236 /*--- STEP 4: test FOM count thresholds to find desired FAR ---*/
2237
2238 /* get the noise simulations back */
2239
2240 if( verb > 1 ) ININFO_message("re-mapping input datasets") ;
2241 remap_Xdataset(xinset) ; /* de-nonce-ification */
2242
2243 /* create voxel-specific FOM threshold images (for each p-value thresh) */
2244
2245 if( do_hpow0 ){
2246 cimar0 = (MRI_IMARR **)malloc(sizeof(MRI_IMARR *)*ncase) ;
2247 car0 = (float ***)malloc(sizeof(float **)*ncase) ;
2248 for( qcase=0 ; qcase < ncase ; qcase++ ){
2249 INIT_IMARR(cimar0[qcase]) ;
2250 car0[qcase] = (float **)malloc(sizeof(float *)*npthr) ;
2251 for( qpthr=0; qpthr < npthr ; qpthr++ ){
2252 cim0 = mri_new_conforming( imtemplate , MRI_float ) ;
2253 ADDTO_IMARR(cimar0[qcase],cim0) ;
2254 car0[qcase][qpthr] = MRI_FLOAT_PTR(cim0) ;
2255 }
2256 }
2257 }
2258
2259 if( do_hpow1 ){
2260 cimar1 = (MRI_IMARR **)malloc(sizeof(MRI_IMARR *)*ncase) ;
2261 car1 = (float ***)malloc(sizeof(float **)*ncase) ;
2262 for( qcase=0 ; qcase < ncase ; qcase++ ){
2263 INIT_IMARR(cimar1[qcase]) ;
2264 car1[qcase] = (float **)malloc(sizeof(float *)*npthr) ;
2265 for( qpthr=0; qpthr < npthr ; qpthr++ ){
2266 cim1 = mri_new_conforming( imtemplate , MRI_float ) ;
2267 ADDTO_IMARR(cimar1[qcase],cim1) ;
2268 car1[qcase][qpthr] = MRI_FLOAT_PTR(cim1) ;
2269 }
2270 }
2271 }
2272
2273 if( do_hpow2 ){
2274 cimar2 = (MRI_IMARR **)malloc(sizeof(MRI_IMARR *)*ncase) ;
2275 car2 = (float ***)malloc(sizeof(float **)*ncase) ;
2276 for( qcase=0 ; qcase < ncase ; qcase++ ){
2277 INIT_IMARR(cimar2[qcase]) ;
2278 car2[qcase] = (float **)malloc(sizeof(float *)*npthr) ;
2279 for( qpthr=0; qpthr < npthr ; qpthr++ ){
2280 cim2 = mri_new_conforming( imtemplate , MRI_float ) ;
2281 ADDTO_IMARR(cimar2[qcase],cim2) ;
2282 car2[qcase][qpthr] = MRI_FLOAT_PTR(cim2) ;
2283 }
2284 }
2285 }
2286
2287 /* array to make map of false alarm count at each voxel */
2288
2289 #ifdef ALLOW_EXTRAS
2290 if( do_FARvox )
2291 farar = (float *)malloc(sizeof(float)*nxyz) ;
2292 #endif
2293
2294 /*--- loop over different FPR (farp) goals ---*/
2295
2296 nit33 = 1.0f/(niter_clust+0.333f) ;
2297
2298 if( ntfp_all == 0 ){ /* only needed if global ETAC was skipped earlier */
2299 ntfp_all = 128 ;
2300 tfs = (float *)malloc(sizeof(float)*ntfp_all) ;
2301 fps = (float *)malloc(sizeof(float)*ntfp_all) ;
2302 }
2303 ntfp = 0 ;
2304
2305 farlast = farplist[0] ;
2306 for( ifarp=0 ; ifarp < numfarp ; ifarp++ ){ /* 23 Aug 2017 */
2307
2308 farp_goal = farplist[ifarp] ;
2309
2310 if( verb )
2311 INFO_message("STEP 4.%s: adjusting per-voxel FOM thresholds to reach FPR=%.2f%% (%.1fs)",
2312 abcd[ifarp] , farp_goal , COX_clock_time() ) ;
2313
2314 /* tfrac = FOM count fractional threshold;
2315 will be adjusted to find the farp_goal FPR goal */
2316
2317 if( ifarp == 0 ){ /* first time thru */
2318 tfrac = (4.0f+farp_goal)*0.00005f ;
2319 } else if( ntfp < 2 ){ /* if only 1 earlier calculation */
2320 tfrac *= 1.0666f * farp_goal / farlast ; /* adjust previous result */
2321 } else {
2322 int jj, jd ; float adf, mdf ; /* later: use closest in history */
2323 mdf = fabsf(fps[0]-farp_goal) ; jd = 0 ; /* to adjust starting point */
2324 for( jj=1 ; jj < ntfp ; jj++ ){
2325 adf = fabsf(fps[jj]-farp_goal) ;
2326 if( adf < mdf ){ mdf = adf ; jd = jj ; }
2327 }
2328 tfrac = tfs[jd] * farp_goal / fps[jd] ;
2329 }
2330
2331 itrac = 0 ;
2332 farpercold = farperc = 0.0f ; tfracold = tfrac ;
2333
2334 /*--- Loop back to here with adjusted tfrac ---*/
2335
2336 FARP_LOOPBACK:
2337 {
2338 float min_tfrac ; int nedge,nmin,ntest ;
2339 min_tfrac = 6.0f / niter_clust ; if( min_tfrac > 0.0001f ) min_tfrac = 0.0001f ;
2340
2341 itrac++ ; /* number of iterations */
2342 nfar = 0 ; /* total FAR count */
2343 nedge = nmin = ntest = 0 ; /* number of edge cases */
2344
2345 /* we take ithresh-th largest FOM at each voxel as its FOM threshold */
2346 #if 0
2347 ithresh = (int)(tfrac*niter_clust) ; /* FOM count threshold */
2348 #else
2349 ithresh = (int)rintf(tfrac*(niter_clust-0.666f)+0.333f) ;
2350 #endif
2351
2352 /* Check if trying to re-litigate previous case [Cinco de Mayo 2017] */
2353
2354 ithresh_list[itrac-1] = ithresh ;
2355 if( itrac > 3 &&
2356 (ithresh_list[itrac-2] == ithresh || ithresh_list[itrac-3] == ithresh) ){
2357 ININFO_message(" #%d: would re-iterate at %g ==> %d ; breaking out",
2358 itrac,tfrac,ithresh ) ;
2359 goto FARP_BREAKOUT ;
2360 }
2361
2362 if( verb )
2363 ININFO_message(" #%d: Testing threshold images at %g ==> %d",itrac,tfrac,ithresh) ;
2364
2365 for( qcase=0 ; qcase < ncase ; qcase++ ){ /* initialize FOM thresholds to 0 */
2366 for( qpthr=0 ; qpthr < npthr ; qpthr++ ){
2367 if( do_hpow0 ){ AAmemset(car0[qcase][qpthr],0,sizeof(float)*nxyz) ; }
2368 if( do_hpow1 ){ AAmemset(car1[qcase][qpthr],0,sizeof(float)*nxyz) ; }
2369 if( do_hpow2 ){ AAmemset(car2[qcase][qpthr],0,sizeof(float)*nxyz) ; }
2370 }}
2371
2372 #ifdef ALLOW_EXTRAS
2373 if( do_FARvox ){
2374 AAmemset(farar,0,sizeof(float)*nxyz) ; /* zero out FARvox array */
2375 }
2376 #endif
2377
2378 AFNI_OMP_START ; /*------------ start parallel section ----------*/
2379 #pragma omp parallel
2380 { DECLARE_ithr ;
2381 MRI_IMAGE *fim , *tfim ; float *far , *tfar ;
2382 int iter,npt,iv,jthresh , ipthr,icase , qq,hh ;
2383 float a0,a1,f0,f1,ft ;
2384 Xvector ***fomsortH=NULL ; float **carH=NULL, *gthreshH ;
2385
2386 /* workspace image array for each thread */
2387
2388 #pragma omp critical
2389 { fim = mri_new_conforming( imtemplate , MRI_float ) ;
2390 far = MRI_FLOAT_PTR(fim) ; }
2391
2392 /* find the FOM threshold for each voxel, for each p-value thresh */
2393
2394 #pragma omp for
2395 for( iv=0 ; iv < mask_ngood ; iv++ ){ /* loop over voxels */
2396 for( icase=0 ; icase < ncase ; icase++ ){ /* over cases */
2397 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){ /* over p-value */
2398 for( hh=0 ; hh < 3 ; hh++ ){ /* over hpow values */
2399 if( hh==0 ){
2400 if( !do_hpow0 ) continue ;
2401 fomsortH = fomsort0[icase]; carH = car0[icase]; gthreshH = gthresh0[ifarp][icase];
2402 }
2403 if( hh==1 ){
2404 if( !do_hpow1 ) continue ;
2405 fomsortH = fomsort1[icase]; carH = car1[icase]; gthreshH = gthresh1[ifarp][icase];
2406 }
2407 if( hh==2 ){
2408 if( !do_hpow2 ) continue ;
2409 fomsortH = fomsort2[icase]; carH = car2[icase]; gthreshH = gthresh2[ifarp][icase];
2410 }
2411 npt = fomsortH[ipthr][iv]->npt ; /* how many FOM values here */
2412 jthresh = ithresh ; /* default index of FOM thresh */
2413
2414 if( jthresh > (int)(0.666f*npt) ){ /* edge case: */
2415 jthresh = (int)(0.666f*npt) ; /* not many FOM values here */
2416 #pragma omp atomic
2417 nedge++ ;
2418 }
2419 /* extract this FOM thresh by interpolation */
2420 #if 0
2421 a0 = ((float)jthresh)/((float)niter_clust) ;
2422 a1 = a0 + 1.0f/((float)niter_clust) ;
2423 #else
2424 a0 = ((float)jthresh+0.666f)*nit33 ;
2425 a1 = a0 + nit33 ;
2426 #endif
2427 f0 = fomsortH[ipthr][iv]->far[jthresh] ;
2428 f1 = fomsortH[ipthr][iv]->far[jthresh+1] ;
2429 ft = inverse_interp_extreme( a0,a1,tfrac , f0,f1 ) ;
2430
2431 #pragma omp atomic /* total number of tests */
2432 ntest++ ;
2433
2434 if( ft < gthreshH[ipthr] ){ /* min case: */
2435 ft = gthreshH[ipthr] ; /* threshold falls below global */
2436 #pragma omp atomic /* min computed for this situation */
2437 nmin++ ;
2438 }
2439 carH[ipthr][ijkmask[iv]] = ft ; /* = FOM threshold for this voxel */
2440 /* = the goal of this entire program! */
2441 } /* end loop over hh */
2442 } /* end loop over p-thresh */
2443 } /* end loop over cases */
2444 } /* end loop over voxels (the parallelized loop) */
2445
2446 /* now, use the thresholds just computed to multi-threshold
2447 each realization, and create count of total FA and in each voxel */
2448
2449 /* #pragma omp for schedule(dynamic,333) */
2450 #pragma omp for /* parallelized */
2451 for( iter=niter_tbot ; iter < niter_ttop ; iter++ ){ /* loop over iterations */
2452 for( icase=0 ; icase < ncase ; icase++ ){ /* over cases */
2453 generate_image( far , iter+icase*niter ) ; /* get image to test */
2454
2455 /* threshold this image */
2456
2457 tfim = mri_multi_threshold_Xcluster_cimar( fim, npthr,zthr_used,
2458 nnsid,nnlev,
2459 do_hpow0 ? cimar0[icase] : NULL ,
2460 do_hpow1 ? cimar1[icase] : NULL ,
2461 do_hpow2 ? cimar2[icase] : NULL ,
2462 0, NULL , NULL ) ;
2463 if( tfim != NULL ){ /* nothing found ==> NULL is returned */
2464 #ifdef ALLOW_EXTRAS
2465 tfar = MRI_FLOAT_PTR(tfim) ;
2466 for( qq=0 ; qq < nxyz ; qq++ ){
2467 if( do_FARvox && tfar[qq] != 0.0f )
2468 #pragma omp atomic
2469 farar[qq]++ ; /* count of FA at this voxel */
2470 }
2471 #endif
2472 #pragma omp critical
2473 { mri_free(tfim) ; nfar++ ;} /* nfar = count of total FA */
2474 break ; /* no need to continue thru cases (got a FA) */
2475 }
2476 } /* end of loop over cases */
2477 } /* end of loop over iterations (parallelized) */
2478
2479 #pragma omp critical
2480 { mri_free(fim) ; } /* tossing the trash for this thread */
2481 }
2482 AFNI_OMP_END ; /*------------ end parallel section ------------*/
2483
2484 /* compute the FAR percentage at this tfrac */
2485
2486 farpercold = farperc ; /* save what we got last time */
2487 farperc = (100.0f*nfar)/(float)niter_test ; /* what we got this time */
2488 farlast = farperc ; /* save for next FPR goal, if any */
2489
2490 /* save results for later re-use [22 Feb 2018] */
2491 if( ntfp == ntfp_all ){
2492 ntfp_all += 128 ;
2493 tfs = (float *)realloc(tfs,sizeof(float)*ntfp_all) ;
2494 fps = (float *)realloc(tfs,sizeof(float)*ntfp_all) ;
2495 }
2496 tfs[ntfp] = tfrac ; fps[ntfp] = farperc ; ntfp++ ;
2497
2498 /* farcut = precision desired for our FPR goal */
2499 farcut = 0.222f ;
2500 if( itrac > 2 ) farcut += (itrac-2)*0.0321f ;
2501
2502 if( verb )
2503 ININFO_message(" local FPR=%.2f%% (%.1fs)", farperc,COX_clock_time()) ;
2504 MEMORY_CHECK(" ") ;
2505
2506 /* if no substantial progress, quit */
2507
2508 if( itrac > 2 && fabsf(farperc-farpercold) < 0.0222f ){
2509 ININFO_message(" ((Progress too slow - breaking out **))") ;
2510 goto FARP_BREAKOUT ;
2511 }
2512
2513 /* try another tfrac to get closer to our goal? */
2514
2515 if( itrac < MAXITE && fabsf(farperc-FG_GOAL) > farcut ){
2516 float fff , dtt ;
2517 if( itrac == 1 || (farperc-FG_GOAL)*(farpercold-FG_GOAL) > 0.1f ){ /* scale */
2518 fff = FG_GOAL/farperc ;
2519 if( fff > 2.222f ) fff = 2.222f ; else if( fff < 0.450f ) fff = 0.450f ;
2520 dtt = (fff-1.0f)*tfrac ; /* tfrac step */
2521 dtt *= (0.8666f+0.2222f*itrac) ; /* accelerate it */
2522 ttemp = tfrac ; tfrac += dtt ;
2523 } else { /* linear inverse interpolate */
2524 fff = (farperc-farpercold)/(tfrac-tfracold) ;
2525 ttemp = tfrac ; tfrac = tfracold + (FG_GOAL-farpercold)/fff ;
2526 }
2527 #define TFTOP (0.666f*TOPFRAC)
2528 tfracold = ttemp ;
2529 if( tfrac < min_tfrac ) tfrac = min_tfrac ;
2530 else if( tfrac > TFTOP ) tfrac = TFTOP ;
2531 goto FARP_LOOPBACK ;
2532 }
2533
2534 FARP_BREAKOUT: ; /*nada*/
2535 } /* end of iterations to find the ideal farperc */
2536
2537 /*=====================================================*/
2538 /*--- Write stuff out, then this farp_goal is done ---*/
2539
2540 gthrout = (float *)malloc(sizeof(float)*npthr*nhpow) ;
2541
2542 for( qcase=0 ; qcase < ncase ; qcase++ ){ /* one dataset per case */
2543
2544 qset = EDIT_empty_copy(mask_dset) ;
2545
2546 sprintf(qpr,".mthresh.%s.%02dperc",lcase[qcase],(int)rintf(farp_goal)) ; /* prefix modifier */
2547
2548 EDIT_dset_items( qset ,
2549 ADN_prefix , modify_afni_prefix(prefix,NULL,qpr) ,
2550 ADN_nvals , npthr*nhpow ,
2551 ADN_datum_all , MRI_float ,
2552 ADN_none ) ;
2553
2554 /* attach bricks for each pthr, for each hpow */
2555
2556 for( ibr=0,qpthr=0 ; qpthr < npthr ; qpthr++ ){
2557 for( hp=0 ; hp < 3 ; hp++ ){
2558 if( hp==0 ){ if( !do_hpow0 ) continue;
2559 else { carHP=car0[qcase]; gthrH=gthresh0[ifarp][qcase]; }}
2560 if( hp==1 ){ if( !do_hpow1 ) continue;
2561 else { carHP=car1[qcase]; gthrH=gthresh1[ifarp][qcase]; }}
2562 if( hp==2 ){ if( !do_hpow2 ) continue;
2563 else { carHP=car2[qcase]; gthrH=gthresh2[ifarp][qcase]; }}
2564 EDIT_substitute_brick( qset , ibr , MRI_float , NULL ) ;
2565 qar = DSET_ARRAY(qset,ibr) ;
2566 AAmemcpy( qar , carHP[qpthr] , sizeof(float)*nxyz ) ;
2567 sprintf(qpr,"Mth:%.4f:h=%d",pthr[qpthr],hp) ;
2568 EDIT_BRICK_LABEL(qset,ibr,qpr) ;
2569 gthrout[ibr] = gthrH[qpthr] ;
2570 ibr++ ;
2571 }}
2572
2573 /* attach an attribute describing the multi-threshold setup */
2574
2575 { float *afl=malloc(sizeof(float)*(npthr+5)) ;
2576 afl[0] = (float)nnlev ; /* NN */
2577 afl[1] = (float)nnsid ; /* sidedness of t-test */
2578 afl[2] = (float)npthr ; /* number of z-thresholds */
2579 afl[3] = (float)(do_hpow0 + 2*do_hpow1 + 4*do_hpow2) ; /* hpow bits */
2580 for( qpthr=0 ; qpthr < npthr ; qpthr++ ) /* and the z-thresholds */
2581 afl[qpthr+4] = zthr_used[qpthr] ;
2582
2583 afl[npthr+4] = (float)min_clust ; /* min cluster size [21 Sep 2017] */
2584
2585 THD_set_float_atr( qset->dblk ,
2586 "MULTI_THRESHOLDS" , npthr+5 , afl ) ;
2587 free(afl) ;
2588
2589 /* attribute for the minimum thresholds per brick [20 Feb 2018] */
2590
2591 THD_set_float_atr( qset->dblk ,
2592 "MULTI_THRESHOLDS_MIN" , npthr*nhpow , gthrout ) ;
2593 }
2594
2595 /* output dataset */
2596
2597 tross_Make_History( "3dXClustSim" , argc,argv , qset ) ;
2598 DSET_write(qset); WROTE_DSETI(qset);
2599 DSET_delete(qset); qset = NULL; qar = NULL;
2600
2601 } /* end of loop over cases (blurs) */
2602
2603 #ifdef ALLOW_EXTRAS
2604 if( do_FARvox ){ /* output the per-voxel false alarm count */
2605 qset = EDIT_empty_copy(mask_dset) ;
2606 sprintf(qpr,".FARvox") ;
2607 EDIT_dset_items( qset ,
2608 ADN_prefix , modify_afni_prefix(prefix,NULL,qpr) ,
2609 ADN_nvals , 1 ,
2610 ADN_none ) ;
2611 EDIT_BRICK_LABEL(qset,0,"FARcount") ;
2612 EDIT_substitute_brick( qset , 0 , MRI_float , NULL ) ;
2613 qar = DSET_ARRAY(qset,0) ;
2614 AAmemcpy( qar , farar , sizeof(float)*nxyz ) ;
2615 tross_Make_History( "3dXClustSim" , argc,argv , qset ) ;
2616 DSET_write(qset); WROTE_DSETI(qset);
2617 DSET_delete(qset); qset = NULL; qar = NULL;
2618 }
2619 #endif
2620
2621 } /* end of loop over farp goals */
2622
2623 /* It's the end of the world, Calvin */
2624
2625 INFO_message("=== 3dXClustSim ends: Elapsed time = %.1fs",COX_clock_time()) ;
2626 MEMORY_CHECK("THE END") ;
2627 exit(0) ;
2628 }
2629
2630 /*============================================================================*/
2631 /*======= Functions to extract the largest FOM found for ===============*/
2632 /*======= each input iteration volume -- used for Global ETAC ===============*/
2633 /*============================================================================*/
2634
2635 /*----------------------------------------------------------------------------*/
2636 /*------- Sorting function for ETAC: int and 3 floats sorted together --------*/
2637 /*------- with the floats going along for the ride driven by the int. --------*/
2638
2639 #undef QS_CUTOFF
2640 #undef QS_SWAP
2641 #undef QS_ISWAP
2642 #define QS_CUTOFF 20 /* cutoff to switch from qsort to isort */
2643 #define QS_SWAP(x,y) (temp=(x), (x)=(y),(y)=temp)
2644 #define QS_ISWAP(x,y) (itemp=(x),(x)=(y),(y)=itemp)
2645 #ifndef QS_STACK
2646 # define QS_STACK 9999
2647 #endif
2648
2649 /*------------------------------------------------------------------------------*/
2650
isort_IFFF(int n,int * iar,float * xar,float * yar,float * zar)2651 void isort_IFFF( int n, int *iar, float *xar, float *yar, float *zar )
2652 {
2653 int j , p ; /* array indices */
2654 float xtemp,ytemp,ztemp ;
2655 int itemp ;
2656 float *xa = xar ;
2657 float *ya = yar ;
2658 float *za = zar ;
2659 int *ia = iar ;
2660
2661 if( n < 2 ) return ;
2662
2663 for( j=1 ; j < n ; j++ ){
2664
2665 if( ia[j] < ia[j-1] ){ /* out of order */
2666 p = j ;
2667 xtemp = xa[j] ; ytemp = ya[j] ; ztemp = za[j] ; itemp = ia[j] ;
2668 do{
2669 xa[p] = xa[p-1]; ya[p] = ya[p-1]; za[p] = za[p-1]; ia[p] = ia[p-1];
2670 p-- ;
2671 } while( p > 0 && itemp < ia[p-1] ) ;
2672 xa[p] = xtemp ; ya[p] = ytemp ; za[p] = ztemp ; ia[p] = itemp ;
2673 }
2674 }
2675 return ;
2676 }
2677
2678 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
2679
qsrec_IFFF(int n,int * iar,float * xar,float * yar,float * zar,int cutoff)2680 void qsrec_IFFF( int n, int *iar, float *xar, float *yar, float *zar, int cutoff )
2681 {
2682 int i , j ; /* scanning indices */
2683 float temp, xpivot, ypivot, zpivot ; /* holding places */
2684 int itemp, ipivot ;
2685 float *xa = xar ;
2686 float *ya = yar ;
2687 float *za = zar ;
2688 int *ia = iar ;
2689
2690 int left , right , mst ;
2691 int stack[QS_STACK] ; /* stack for qsort "recursion" */
2692
2693 /* return if too short (insertion sort will clean up) */
2694
2695 if( cutoff < 3 ) cutoff = 3 ;
2696 if( n < cutoff ) return ;
2697
2698 /* initialize stack to start with whole array */
2699
2700 stack[0] = 0 ;
2701 stack[1] = n-1 ;
2702 mst = 2 ;
2703
2704 /* loop while the stack is nonempty */
2705
2706 while( mst > 0 ){
2707 right = stack[--mst] ; /* work on subarray from left -> right */
2708 left = stack[--mst] ;
2709
2710 i = ( left + right ) / 2 ; /* middle of subarray */
2711
2712 /*----- sort the left, middle, and right a[]'s -----*/
2713
2714 if( ia[left] > ia[i] ){ QS_SWAP (xa[left] ,xa[i] ) ;
2715 QS_SWAP (ya[left] ,ya[i] ) ;
2716 QS_SWAP (za[left] ,za[i] ) ;
2717 QS_ISWAP(ia[left] ,ia[i] ) ; }
2718 if( ia[left] > ia[right] ){ QS_SWAP (xa[left] ,xa[right] ) ;
2719 QS_SWAP (ya[left] ,ya[right] ) ;
2720 QS_SWAP (za[left] ,za[right] ) ;
2721 QS_ISWAP(ia[left] ,ia[right] ) ; }
2722 if( ia[i] > ia[right] ){ QS_SWAP (xa[right] ,xa[i] ) ;
2723 QS_SWAP (ya[right] ,ya[i] ) ;
2724 QS_SWAP (za[right] ,za[i] ) ;
2725 QS_ISWAP(ia[right] ,ia[i] ) ; }
2726
2727 xpivot = xa[i] ; xa[i] = xa[right] ;
2728 ypivot = ya[i] ; ya[i] = ya[right] ;
2729 zpivot = za[i] ; za[i] = za[right] ;
2730 ipivot = ia[i] ; ia[i] = ia[right] ;
2731
2732 i = left ; j = right ; /* initialize scanning */
2733
2734 /*----- partition: move elements bigger than pivot up and elements
2735 smaller than pivot down, scanning in from ends -----*/
2736
2737 do{
2738 for( ; ia[++i] < ipivot ; ) ; /* scan i up, until ia[i] >= ipivot */
2739 for( ; ia[--j] > ipivot ; ) ; /* scan j down, until ia[j] <= ipivot */
2740
2741 if( j <= i ) break ; /* if j meets i, quit this loop */
2742
2743 QS_SWAP (xa[i] ,xa[j] ) ;
2744 QS_SWAP (ya[i] ,ya[j] ) ;
2745 QS_SWAP (za[i] ,za[j] ) ;
2746 QS_ISWAP(ia[i] ,ia[j] ) ;
2747 } while( 1 ) ;
2748
2749 /*----- at this point, the array is partitioned at the pivot -----*/
2750
2751 xa[right] = xa[i] ; xa[i] = xpivot ; /* restore the pivot */
2752 ya[right] = ya[i] ; ya[i] = ypivot ;
2753 za[right] = za[i] ; za[i] = zpivot ;
2754 ia[right] = ia[i] ; ia[i] = ipivot ;
2755
2756 /*----- signal the subarrays that need further work -----*/
2757
2758 if( (i-left) > cutoff ){ stack[mst++] = left ; stack[mst++] = i-1 ; }
2759 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right ; }
2760
2761 } /* end of while stack is non-empty */
2762 return ;
2763 }
2764
2765 /* quick_sort : sort an array partially recursively, and partially insertion */
2766
qsort_IFFF(int n,int * iar,float * xar,float * yar,float * zar)2767 void qsort_IFFF( int n, int *iar, float *xar, float *yar, float *zar )
2768 {
2769 if( n < 2 || iar == NULL || xar == NULL || yar == NULL || zar == NULL ) return ;
2770 qsrec_IFFF( n , iar , xar,yar,zar , QS_CUTOFF ) ;
2771 isort_IFFF( n , iar , xar,yar,zar ) ;
2772 return ;
2773 }
2774
2775 /*----------------------------------------------------------------------------*/
2776 /* Find the max FOM value that occured in any given iteration (iar).
2777 * inputs are the sorted-by-iar [xyz]ar arrays;
2778 * outputs [xyz]out are the largest [xyz]ar values for each
2779 distinct value of iar;
2780 * return value is the number of such distinct values found
2781 == number of values loaded into [xyz]out arrays.
2782 *//*--------------------------------------------------------------------------*/
2783
ixmax_IFFF(int n,int * iar,float * xar,float * yar,float * zar,float * xout,float * yout,float * zout)2784 int ixmax_IFFF( int n, int *iar, float *xar, float *yar, float *zar,
2785 float *xout, float *yout, float *zout )
2786 {
2787 int nfound , ival , jj,kk ;
2788 float xtop, ytop, ztop ;
2789
2790 for( nfound=jj=0 ; jj < n ; ){
2791 ival = iar[jj] ; xtop = xar[jj] ; ytop = yar[jj] ; ztop = zar[jj] ;
2792 for( ++jj ; jj < n && iar[jj] == ival ; jj++ ){
2793 if( xar[jj] > xtop ) xtop = xar[jj] ;
2794 if( yar[jj] > ytop ) ytop = yar[jj] ;
2795 if( zar[jj] > ztop ) ztop = zar[jj] ;
2796 }
2797 xout[nfound] = xtop; yout[nfound] = ytop; zout[nfound] = ztop; nfound++;
2798 }
2799
2800 return nfound ;
2801 }
2802 /*============================================================================*/
2803