1 /*** Adapted from 3dClustSim.c ***/
2
3 #include "mrilib.h"
4
5 #include <time.h>
6 #include <sys/types.h>
7 #include <unistd.h>
8
9 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
10 #define MAX_CLUSTER_SIZE 99999 /* max. size of cluster for freq. table */
11
12 /*---------------------------------------------------------------------------*/
13 /* Global data */
14
15 static int max_cluster_size = MAX_CLUSTER_SIZE ;
16
17 static int nx ;
18 static int ny ;
19 static int nz ;
20 static int nxy ;
21 static int nxyz ;
22
23 static int do_niml = 1 ;
24 static int do_1D = 1 ;
25 static int do_2sid = 1 ;
26
27 #define PMAX 0.2
28
29 static int npthr = 29 ;
30 static double pthr[] = { 0.10, 0.09, 0.08, 0.07, 0.06,
31 0.05, 0.04, 0.03, 0.02, 0.015, 0.01,
32 0.007, 0.005, 0.003, 0.002, 0.0015, 0.001,
33 0.0007, 0.0005, 0.0003, 0.0002, 0.00015, 0.0001,
34 0.00007, 0.00005, 0.00003, 0.00002, 0.000015, 0.00001 } ;
35
36 static int nathr = 10 ;
37 static double athr[] = { 0.10, 0.09, .08, .07, .06, .05, .04, .03, .02, .01 } ;
38
39 static int verb = 1 ;
40 static int nthr = 1 ;
41
42 static char *prefix = NULL ;
43 static int do_final = 0 ;
44
45 static int ndset = 0 ;
46 static char **fdset = NULL ;
47
48 #define do_lohi 0
49 #define nodec 0
50
51 /*---------------------------------------------------------------------------*/
52
display_help_menu()53 void display_help_menu()
54 {
55 printf(
56 "Usage: 3dClustCount [options] dataset1 ... \n"
57 "\n"
58 "This program takes as input 1 or more datasets, thresholds them at various\n"
59 "levels, and counts up the number of clusters of various sizes. It is\n"
60 "adapted from 3dClustSim, but only does the cluster counting functions --\n"
61 "where the datasets come from is the user's business. It is intended for\n"
62 "use in a simulation script.\n"
63 #if 0
64 "\n"
65 "The input datasets can be in one of two forms.\n"
66 "(1) A single sub-brick, which should be coded as a statistic so that\n"
67 " a voxel-wise (AKA uncorrected) 2-sided p-value threshold can be\n"
68 " applied.\n"
69 "(2) Two sub-bricks, the first of which is some voxel-wise measurement\n"
70 " of the effect (e.g., percent signal change), and the second of which\n"
71 " is the statistical threshold sub-brick.\n"
72 "In case (1), only statistics of the cluster sizes are kept.\n"
73 "In case (2), statistics of the effect sub-brick are also kept for the\n"
74 "clusters of each size.\n"
75 #endif
76 "\n"
77 "-------\n"
78 "OPTIONS\n"
79 "-------\n"
80 " -prefix sss = Use string 'sss' as the prefix of the filename into which\n"
81 " results will be summed. The actual filename will be\n"
82 " 'sss.clustcount.niml'. If this file already exists, then\n"
83 " the results from the current run will be summed into the\n"
84 " existing results, and the file then re-written.\n"
85 #if 0
86 "\n"
87 " -addin aaa = Also add in the results from file 'aaa.clustcount.niml'\n"
88 " to the cumulating counts. This option lets you run\n"
89 " multiple 3dClustCount jobs in parallel scripts, then\n"
90 " merge the results.\n"
91 #endif
92 "\n"
93 " -final = If this option is given, then the results will be output\n"
94 " in a format like that used from 3dClustSim -- as 1D and\n"
95 " NIML formatted files with probabilities of various\n"
96 " cluster sizes.\n"
97 " ++ You can use '-final' without any input datasets if\n"
98 " you want to create the final output files from the\n"
99 " saved '.clustcount.niml' output file from earlier runs.\n"
100 "\n"
101 " -quiet = Don't print out the progress reports, etc.\n"
102 " ++ Put this option first to quiet most informational messages.\n"
103 "\n"
104 "--------\n"
105 "EXAMPLE:\n"
106 "-------\n"
107 "The steps here are\n"
108 " (a) Create a set of 250 3dGroupInCorr results from a set of 190 subjects,\n"
109 " using 250 randomly located seed locations. Note the use of '-sendall'\n"
110 " to get the individual subject results -- these are used in the next\n"
111 " step, and are in sub-bricks 2..191 -- the collective 3dGroupInCorr\n"
112 " results (in sub-bricks 0..1) are not actually used here.\n"
113 " (b) For each of these 250 output datasets, create 80 random splittings\n"
114 " into 2 groups of 95 subjects each, and carry out a 2-sample t-test\n"
115 " between these groups.\n"
116 " ++ Note the use of program 2perm to create the random splittings into\n"
117 " files QQ_A and QQ_B, drawn from sub-bricks 2..191 of the ${fred}\n"
118 " datasets.\n"
119 " ++ Note the use of the '[1dcat filename]' construction to specify\n"
120 " which sub-bricks of the ${fred} dataset are used for input to\n"
121 " the '-setX' options of 3dttest++.\n"
122 " (c) Count clusters from the '[1]' sub-brick of the 80 t-test outputs --\n"
123 " the t-statistic sub-brick.\n"
124 " ++ Note the use of a wildcard filename with a sub-brick selector:\n"
125 " 'QQ*.HEAD[1]' -- 3dClustCount will do the wildcard expansion\n"
126 " internally, then add the sub-brick selector '[1]' to each expanded\n"
127 " dataset filename.\n"
128 " (d) Produce the final report files for empirical cluster-size thresholds\n"
129 " for 3dGroupInCorr analyses -- rather than rely on 3dClustSim's assumption\n"
130 " of Gaussian-shaped spatial correlation structure.\n"
131 "The syntax is C-shell (tcsh), naturally.\n"
132 "\n"
133 " \\rm -f ABscat*\n"
134 " 3dGroupInCorr -setA A.errts.grpincorr.niml \\\n"
135 " -setB B.errts.grpincorr.niml \\\n"
136 " -labelA A -labelB B -seedrad 5 -nosix -sendall \\\n"
137 " -batchRAND 250 ABscat\n"
138 " foreach fred ( ABscat*.HEAD )\n"
139 " foreach nnn ( `count -dig 2 0 79` )\n"
140 " 2perm -prefix QQ 2 191\n"
141 " 3dttest++ -setA ${fred}'[1dcat QQ_A]' \\\n"
142 " -setB ${fred}'[1dcat QQ_B]' \\\n"
143 " -no1sam -prefix QQ${nnn}\n"
144 " end\n"
145 " 3dClustCount -prefix ABcount 'QQ*.HEAD[1]'\n"
146 " \\rm -f QQ*\n"
147 " end\n"
148 " 3dClustCount -final -prefix ABcount\n"
149 " \\rm -f ABscat*\n"
150 "\n"
151 "--------------------------------\n"
152 "---- RW Cox -- August 2012 -----\n"
153 "--------------------------------\n"
154 ) ;
155
156 PRINT_COMPILE_DATE ;
157 exit(0);
158 }
159
160 /*---------------------------------------------------------------------------*/
161 /* Routine to initialize the input options (values are in global variables). */
162
get_options(int argc,char ** argv)163 void get_options( int argc , char **argv )
164 {
165 int nopt=1 , ii ;
166 int nexp,iex,didex ; char **fexp ;
167
168 PUTENV("AFNI_GLOB_SELECTORS","YES") ;
169
170 while( nopt < argc && argv[nopt][0] == '-' ){
171
172 /*----- -input -----*/
173
174 if( strcasecmp(argv[nopt],"-input") == 0 ){
175
176 if( ++nopt >= argc ) ERROR_exit("need argument after '%s'",argv[nopt-1]) ;
177
178 for( ; nopt < argc && argv[nopt][0] != '-' ; nopt++ ){
179 if( HAS_WILDCARD(argv[nopt]) ){
180 MCW_file_expand( 1,argv+nopt, &nexp,&fexp ) ; didex = 1 ;
181 if( nexp < 1 ){
182 WARNING_message("wildcard name '%s' doesn't match any files",argv[nopt]) ;
183 continue ;
184 }
185 } else {
186 nexp = 1 ; fexp = argv+nopt ; didex = 0 ;
187 }
188 fdset = (char **)realloc(fdset,sizeof(char *)*(ndset+nexp)) ;
189 for( iex=0 ; iex < nexp ; iex++ ) fdset[ndset++] = strdup(fexp[iex]) ;
190 if( didex ) MCW_free_expand(nexp,fexp) ;
191 }
192 continue ;
193 }
194
195 #if 0
196 /*----- -2sided ------*/
197
198 if( strcasecmp(argv[nopt],"-2sided") == 0 ){
199 do_2sid = 1 ; nopt++ ; continue ;
200 }
201 if( strcasecmp(argv[nopt],"-1sided") == 0 ){
202 do_2sid = 0 ; nopt++ ; continue ;
203 }
204 #endif
205
206 /*----- -prefix -----*/
207
208 if( strcmp(argv[nopt],"-prefix") == 0 ){
209 nopt++ ; if( nopt >= argc ) ERROR_exit("need argument after -prefix!") ;
210 prefix = strdup(argv[nopt]) ;
211 if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix option!") ;
212 nopt++ ; continue ;
213 }
214
215 /*---- -quiet ----*/
216
217 if( strcasecmp(argv[nopt],"-quiet") == 0 ){
218 verb = 0 ; nopt++ ; continue ;
219 }
220
221 /*----- -final -----*/
222
223 if( strcasecmp(argv[nopt],"-final") == 0 ){
224 do_final = 1 ; nopt++ ; continue ;
225 }
226
227 /*----- unknown option -----*/
228
229 ERROR_exit("3dClustCount -- unknown option '%s'",argv[nopt]) ;
230 }
231
232 /*----- anything left is a dataset filename -----*/
233
234 for( ; nopt < argc ; nopt++ ){
235 if( HAS_WILDCARD(argv[nopt]) ){
236 MCW_file_expand( 1,argv+nopt, &nexp,&fexp ) ; didex = 1 ;
237 if( nexp < 1 ){
238 WARNING_message("wildcard name '%s' doesn't match any files",argv[nopt]) ;
239 continue ;
240 }
241 } else {
242 nexp = 1 ; fexp = argv+nopt ; didex = 0 ;
243 }
244 fdset = (char **)realloc(fdset,sizeof(char *)*(ndset+nexp)) ;
245 for( iex=0 ; iex < nexp ; iex++ ) fdset[ndset++] = strdup(fexp[iex]) ;
246 if( didex ) MCW_free_expand(nexp,fexp) ;
247 }
248
249 if( ndset == 0 && !do_final ) ERROR_exit("No input datasets?!") ;
250
251 if( prefix == NULL ){
252 prefix = "ClustCount" ;
253 INFO_message("No -prefix option ==> using prefix = '%s'",prefix) ;
254 }
255
256 return ;
257 }
258
259 /*---------------------------------------------------------------------------*/
260
261 #define DALL MAX_CLUSTER_SIZE
262
263 /*! Put (i,j,k) into the current cluster, if it is nonzero. */
264
265 #define CPUT(i,j,k) \
266 do{ ijk = THREE_TO_IJK(i,j,k,nx,nxy) ; \
267 if( mmm[ijk] ){ \
268 if( nnow == nall ){ /* increase array lengths */ \
269 nall += DALL + nall/2 ; \
270 inow = (short *)realloc(inow,sizeof(short)*nall) ; \
271 jnow = (short *)realloc(jnow,sizeof(short)*nall) ; \
272 know = (short *)realloc(know,sizeof(short)*nall) ; \
273 } \
274 inow[nnow] = i ; jnow[nnow] = j ; know[nnow] = k ; \
275 nnow++ ; mmm[ijk] = 0 ; \
276 } } while(0)
277
278 #define USE_MEMCHR /* faster, but trickier to understand */
279
280 static int nall = 0 ; /* arrays for clusterizing */
281 static short *inow = NULL ;
282 static short *jnow = NULL ;
283 static short *know = NULL ;
284
285 /*----------------------------------------------------------------------------*/
286
find_largest_cluster_NN1(byte * mmm)287 int find_largest_cluster_NN1( byte *mmm )
288 {
289 int max_size=0 ;
290 int ii,jj,kk, icl , ijk , ijk_last ;
291 int ip,jp,kp , im,jm,km , nnow ;
292 #ifdef USE_MEMCHR
293 byte *mch ;
294 #endif
295
296 ijk_last = 0 ; /* start scanning at the start of voxel-land */
297
298 while(1) {
299 /* find next nonzero point in mmm array */
300
301 #ifndef USE_MEMCHR
302 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
303 #else
304 mch = memchr( mmm+ijk_last , 1 , nxyz-ijk_last ) ; /* quicker search */
305 if( mch == NULL ) ijk = nxyz ;
306 else ijk = mch - mmm ;
307 #endif
308 if( ijk >= nxyz ) break ; /* didn't find any! */
309 ijk_last = ijk+1 ; /* start here next time */
310
311 mmm[ijk] = 0 ; /* clear found point */
312
313 nnow = 1 ; /* # pts in cluster */
314 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
315
316 /* loop over points in cluster, checking their neighbors,
317 growing the cluster if we find any that belong therein */
318
319 for( icl=0 ; icl < nnow ; icl++ ){
320 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
321 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
322 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
323
324 if( im >= 0 ) CPUT(im,jj,kk) ; /* NN1 */
325 if( ip < nx ) CPUT(ip,jj,kk) ;
326 if( jm >= 0 ) CPUT(ii,jm,kk) ;
327 if( jp < ny ) CPUT(ii,jp,kk) ;
328 if( km >= 0 ) CPUT(ii,jj,km) ;
329 if( kp < nz ) CPUT(ii,jj,kp) ;
330 }
331
332 if( nnow > max_size ) max_size = nnow ;
333 }
334
335 return max_size ;
336 }
337
338 /*----------------------------------------------------------------------------*/
339
find_largest_cluster_NN2(byte * mmm)340 int find_largest_cluster_NN2( byte *mmm )
341 {
342 int max_size=0 ;
343 int ii,jj,kk, icl , ijk , ijk_last ;
344 int ip,jp,kp , im,jm,km , nnow ;
345 #ifdef USE_MEMCHR
346 byte *mch ;
347 #endif
348
349 ijk_last = 0 ; /* start scanning at the start */
350
351 while(1) {
352 /* find next nonzero point in mmm array */
353
354 #ifndef USE_MEMCHR
355 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
356 #else
357 mch = memchr( mmm+ijk_last , 1 , nxyz-ijk_last ) ; /* quicker search */
358 if( mch == NULL ) ijk = nxyz ;
359 else ijk = mch - mmm ;
360 #endif
361 if( ijk == nxyz ) break ; /* didn't find any! */
362 ijk_last = ijk+1 ; /* start here next time */
363
364 mmm[ijk] = 0 ; /* clear found point */
365
366 nnow = 1 ; /* # pts in cluster */
367 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
368
369 /* loop over points in cluster, checking their neighbors,
370 growing the cluster if we find any that belong therein */
371
372 for( icl=0 ; icl < nnow ; icl++ ){
373 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
374 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
375 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
376
377 if( im >= 0 ){ CPUT(im,jj,kk) ; /* NN1 */
378 if( jm >= 0 ) CPUT(im,jm,kk) ; /* NN2 */
379 if( jp < ny ) CPUT(im,jp,kk) ; /* NN2 */
380 if( km >= 0 ) CPUT(im,jj,km) ; /* NN2 */
381 if( kp < nz ) CPUT(im,jj,kp) ; /* NN2 */
382 }
383 if( ip < nx ){ CPUT(ip,jj,kk) ; /* NN1 */
384 if( jm >= 0 ) CPUT(ip,jm,kk) ; /* NN2 */
385 if( jp < ny ) CPUT(ip,jp,kk) ; /* NN2 */
386 if( km >= 0 ) CPUT(ip,jj,km) ; /* NN2 */
387 if( kp < nz ) CPUT(ip,jj,kp) ; /* NN2 */
388 }
389 if( jm >= 0 ){ CPUT(ii,jm,kk) ; /* NN1 */
390 if( km >= 0 ) CPUT(ii,jm,km) ; /* NN2 */
391 if( kp < nz ) CPUT(ii,jm,kp) ; /* NN2 */
392 }
393 if( jp < ny ){ CPUT(ii,jp,kk) ; /* NN1 */
394 if( km >= 0 ) CPUT(ii,jp,km) ; /* NN2 */
395 if( kp < nz ) CPUT(ii,jp,kp) ; /* NN2 */
396 }
397 if( km >= 0 ) CPUT(ii,jj,km) ; /* NN1 */
398 if( kp < nz ) CPUT(ii,jj,kp) ; /* NN1 */
399 }
400
401 if( nnow > max_size ) max_size = nnow ;
402 }
403
404 return max_size ;
405 }
406
407 /*----------------------------------------------------------------------------*/
408
find_largest_cluster_NN3(byte * mmm)409 int find_largest_cluster_NN3( byte *mmm )
410 {
411 int max_size=0 ;
412 int ii,jj,kk, icl , ijk , ijk_last ;
413 int ip,jp,kp , im,jm,km , nnow ;
414 #ifdef USE_MEMCHR
415 byte *mch ;
416 #endif
417
418 ijk_last = 0 ; /* start scanning at the start */
419
420 while(1) {
421 /* find next nonzero point in mmm array */
422
423 #ifndef USE_MEMCHR
424 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( mmm[ijk] ) break ;
425 #else
426 mch = memchr( mmm+ijk_last , 1 , nxyz-ijk_last ) ; /* quicker search */
427 if( mch == NULL ) ijk = nxyz ;
428 else ijk = mch - mmm ;
429 #endif
430 if( ijk == nxyz ) break ; /* didn't find any! */
431 ijk_last = ijk+1 ; /* start here next time */
432
433 mmm[ijk] = 0 ; /* clear found point */
434
435 nnow = 1 ; /* # pts in cluster */
436 IJK_TO_THREE(ijk, inow[0],jnow[0],know[0] , nx,nxy) ;
437
438 /* loop over points in cluster, checking their neighbors,
439 growing the cluster if we find any that belong therein */
440
441 for( icl=0 ; icl < nnow ; icl++ ){
442 ii = inow[icl] ; jj = jnow[icl] ; kk = know[icl] ;
443 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
444 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
445
446 if( im >= 0 ){ CPUT(im,jj,kk) ; /* NN1 */
447 if( jm >= 0 ) CPUT(im,jm,kk) ; /* NN2 */
448 if( jp < ny ) CPUT(im,jp,kk) ; /* NN2 */
449 if( km >= 0 ) CPUT(im,jj,km) ; /* NN2 */
450 if( kp < nz ) CPUT(im,jj,kp) ; /* NN2 */
451 if( jm >= 0 && km >= 0 ) CPUT(im,jm,km) ; /* NN3 */
452 if( jm >= 0 && kp < nz ) CPUT(im,jm,kp) ; /* NN3 */
453 if( jp < ny && km >= 0 ) CPUT(im,jp,km) ; /* NN3 */
454 if( jp < ny && kp < nz ) CPUT(im,jp,kp) ; /* NN3 */
455 }
456 if( ip < nx ){ CPUT(ip,jj,kk) ; /* NN1 */
457 if( jm >= 0 ) CPUT(ip,jm,kk) ; /* NN2 */
458 if( jp < ny ) CPUT(ip,jp,kk) ; /* NN2 */
459 if( km >= 0 ) CPUT(ip,jj,km) ; /* NN2 */
460 if( kp < nz ) CPUT(ip,jj,kp) ; /* NN2 */
461 if( jm >= 0 && km >= 0 ) CPUT(ip,jm,km) ; /* NN3 */
462 if( jm >= 0 && kp < nz ) CPUT(ip,jm,kp) ; /* NN3 */
463 if( jp < ny && km >= 0 ) CPUT(ip,jp,km) ; /* NN3 */
464 if( jp < ny && kp < nz ) CPUT(ip,jp,kp) ; /* NN3 */
465 }
466 if( jm >= 0 ){ CPUT(ii,jm,kk) ; /* NN1 */
467 if( km >= 0 ) CPUT(ii,jm,km) ; /* NN2 */
468 if( kp < nz ) CPUT(ii,jm,kp) ; /* NN2 */
469 }
470 if( jp < ny ){ CPUT(ii,jp,kk) ; /* NN1 */
471 if( km >= 0 ) CPUT(ii,jp,km) ; /* NN2 */
472 if( kp < nz ) CPUT(ii,jp,kp) ; /* NN2 */
473 }
474 if( km >= 0 ) CPUT(ii,jj,km) ; /* NN1 */
475 if( kp < nz ) CPUT(ii,jj,kp) ; /* NN1 */
476 }
477
478 if( nnow > max_size ) max_size = nnow ;
479 }
480
481 return max_size ;
482 }
483
484 /*---------------------------------------------------------------------------*/
485 /* Find clusters, save some info, re-populate array? */
486
gather_stats_NN1(float thr,float * fim,byte * bfim,int * mtab)487 void gather_stats_NN1( float thr , float *fim , byte *bfim , int *mtab )
488 {
489 register int ii ; int siz ;
490
491 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] >= thr) ;
492 siz = find_largest_cluster_NN1(bfim) ;
493 if( siz > max_cluster_size ) siz = max_cluster_size ;
494 if( siz > 0 ) mtab[siz]++ ;
495
496 return ;
497 }
498
499 /*---------------------------------------------------------------------------*/
500 /* Find clusters, save some info, re-populate array? */
501
gather_stats_NN2(float thr,float * fim,byte * bfim,int * mtab)502 void gather_stats_NN2( float thr , float *fim , byte *bfim , int *mtab )
503 {
504 register int ii ; int siz ;
505
506 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] >= thr) ;
507 siz = find_largest_cluster_NN2(bfim) ;
508 if( siz > max_cluster_size ) siz = max_cluster_size ;
509 if( siz > 0 ) mtab[siz]++ ;
510
511 return ;
512 }
513
514 /*---------------------------------------------------------------------------*/
515 /* Find clusters, save some info, re-populate array? */
516
gather_stats_NN3(float thr,float * fim,byte * bfim,int * mtab)517 void gather_stats_NN3( float thr , float *fim , byte *bfim , int *mtab )
518 {
519 register int ii ; int siz ;
520
521 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (fim[ii] >= thr) ;
522 siz = find_largest_cluster_NN3(bfim) ;
523 if( siz > max_cluster_size ) siz = max_cluster_size ;
524 if( siz > 0 ) mtab[siz]++ ;
525
526 return ;
527 }
528
529 /*===========================================================================*/
530
main(int argc,char ** argv)531 int main( int argc , char **argv )
532 {
533 int **max_table[4] ; int nnn , ipthr , first_mask=1 , niter ;
534 char *fnam ;
535 NI_element *neldat=NULL ;
536
537 /*----- does user request help menu? -----*/
538
539 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) display_help_menu() ;
540
541 /*----- get the list of things to do -----*/
542
543 mainENTRY("3dClustCount"); machdep();
544 AFNI_logger("3dClustCount",argc,argv);
545 PRINT_VERSION("3dClustCount"); AUTHOR("Zhark the Enumerator");
546
547 get_options( argc , argv ) ;
548
549 /*----- create some space for the results -----*/
550
551 for( nnn=1 ; nnn <= 3 ; nnn++ ){
552 max_table[nnn] = (int **)malloc(sizeof(int *)*npthr) ; /* array of tables */
553 for( ipthr=0 ; ipthr < npthr ; ipthr++ ) /* create tables */
554 max_table[nnn][ipthr] = (int *)calloc(sizeof(int),(max_cluster_size+1)) ;
555 }
556
557 /* in 3dClustSim, the corresponding {...} code is inside OpenMP */
558 {
559 int ipthr, **mt[4] , nnn , ids,nds,ival , statcode , ii ;
560 MRI_IMAGE *fim = NULL; float *far, *stataux , thr ; byte *bfar ;
561 THD_3dim_dataset *dset ;
562
563 nall = DALL ;
564 inow = (short *)malloc(sizeof(short)*DALL) ;
565 jnow = (short *)malloc(sizeof(short)*DALL) ;
566 know = (short *)malloc(sizeof(short)*DALL) ;
567 for( nnn=1 ; nnn <= 3 ; nnn++ ) mt[nnn] = max_table[nnn] ;
568
569 /*--- loop over datasets ---*/
570
571 for( niter=ids=0 ; ids < ndset ; ids++ ){
572
573 /* get the next dataset */
574
575 dset = THD_open_dataset( fdset[ids] ) ;
576 if( !ISVALID_DSET(dset) ){
577 ERROR_message("Can't open dataset %s -- skipping it :-(",fdset[ids]) ;
578 continue ;
579 }
580 DSET_load(dset) ;
581 if( !DSET_LOADED(dset) ){
582 ERROR_message("Can't load dataset %s -- skipping it :-(",fdset[ids]) ;
583 DSET_delete(dset) ;
584 }
585 nx = DSET_NX(dset) ;
586 ny = DSET_NY(dset) ;
587 nz = DSET_NZ(dset) ; nxy = nx*ny ; nxyz = nx*ny*nz ;
588 bfar = (byte * )malloc(sizeof(byte)*nxyz) ;
589
590 if( verb )
591 fprintf(stderr,"++ Dataset %s",fdset[ids]) ;
592
593 /* process each statistical sub-brick */
594
595 for( nds=ival=0 ; ival < DSET_NVALS(dset) ; ival++ ){
596 statcode = DSET_BRICK_STATCODE(dset,ival) ; /* type of data */
597 stataux = DSET_BRICK_STATAUX (dset,ival) ;
598 if( !FUNC_IS_STAT(statcode) ) continue ; /* bad */
599 fim = THD_extract_float_brick(ival,dset) ;
600 if( fim == NULL ) continue ; /* bad */
601 far = MRI_FLOAT_PTR(fim) ;
602 for( ii=0 ; ii < nxyz ; ii++ ) far[ii] = fabsf(far[ii]) ;
603 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
604 thr = THD_pval_to_stat( pthr[ipthr] , statcode , stataux ) ;
605 gather_stats_NN1( thr , far , bfar , mt[1][ipthr] ) ;
606 gather_stats_NN2( thr , far , bfar , mt[2][ipthr] ) ;
607 gather_stats_NN3( thr , far , bfar , mt[3][ipthr] ) ;
608 }
609 niter++ ; nds++ ; /* one more iteration in the bag */
610 if( verb && nds == 1 ) fprintf(stderr," -- ") ;
611 } /* end of loop over sub-bricks */
612
613 if( verb && nds > 0 )
614 fprintf(stderr,"%d / %d sub-bricks processed -- %d total now\n",
615 nds,DSET_NVALS(dset),niter) ;
616 else if( verb && nds == 0 )
617 fprintf(stderr," -- 0 sub-bricks processed\n") ;
618
619 mri_free(fim) ; free(bfar) ; DSET_delete(dset) ;
620
621 } /* end loop over datasets */
622
623 free(know) ; free(jnow) ; free(inow) ;
624
625 } /* end of processing the inputs */
626
627 if( niter == 0 && !do_final )
628 ERROR_exit("No data was processed :-(") ;
629
630 if( verb && ndset > 1 && niter > 1 )
631 INFO_message("Processed a total of %d sub-bricks from %d datasets",niter,ndset) ;
632
633 /*---------- now read the existing file, if any, and merge it ----------*/
634
635 fnam = (char *)malloc(sizeof(char)*(strlen(prefix)+128)) ;
636 strcpy(fnam,prefix) ;
637 if( strstr(fnam,".clustcount.niml") == NULL ) strcat(fnam,".clustcount.niml") ;
638 neldat = NI_read_element_fromfile(fnam) ;
639 if( neldat != NULL ){
640
641 #define NELERR(sf) \
642 do{ ERROR_message("File %s has bad data [%s]",fnam,sf) ; \
643 NI_free_element(neldat) ; neldat = NULL ; goto NelDone ; } while(0)
644
645 char *atr ;
646 int ii,jj, nelite, nelmax, nelnum, **vel[4]={NULL,NULL,NULL,NULL} ;
647
648 if( neldat->type != NI_ELEMENT_TYPE ) NELERR("type") ;
649
650 nelmax = neldat->vec_len ;
651 if( nelmax > max_cluster_size ) nelmax = max_cluster_size+1 ;
652 if( nelmax < 2 ) NELERR("length") ;
653
654 nelnum = neldat->vec_num ;
655 if( nelnum < 3*npthr ) NELERR("number") ;
656 for( nnn=1 ; nnn <= 3 ; nnn++ ){
657 vel[nnn] = (int **)malloc(sizeof(int *)*npthr) ;
658 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
659 jj = (nnn-1)*npthr + ipthr ;
660 if( neldat->vec_typ[jj] != NI_INT ) NELERR("datum") ;
661 vel[nnn][ipthr] = (int *)neldat->vec[jj] ;
662 }
663 }
664
665 atr = NI_get_attribute( neldat , "niter" ) ;
666 if( atr == NULL ) NELERR("niter") ;
667
668 nelite = (int)strtod(atr,NULL) ;
669 if( nelite <= 0 ) NELERR("niter") ;
670
671 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
672 for( ii=0 ; ii < nelmax ; ii++ ){
673 max_table[1][ipthr][ii] += vel[1][ipthr][ii] ;
674 max_table[2][ipthr][ii] += vel[2][ipthr][ii] ;
675 max_table[3][ipthr][ii] += vel[3][ipthr][ii] ;
676 }
677 }
678
679 niter += nelite ;
680 if( verb ) INFO_message("merged %d sub-brick counts from file %s",nelite,fnam) ;
681
682 #undef NELERR
683 NelDone:
684 if( vel[1] != NULL ){ free(vel[1]); free(vel[2]); free(vel[3]); }
685 NI_free_element(neldat) ; neldat = NULL ;
686 }
687
688 /*---------- save counting data to file ----------*/
689
690 if( niter > 0 ){
691 int itop=0 , ii ; char buf[128] ;
692
693 for( ii=max_cluster_size ; ii > 0 && itop == 0 ; ii-- ){
694 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
695 if( max_table[1][ipthr][ii] > 0 ||
696 max_table[2][ipthr][ii] > 0 ||
697 max_table[3][ipthr][ii] > 0 ){ itop=ii; break; }
698 }
699 }
700 if( itop == 0 ) ERROR_exit("No clusters found at all!") ;
701
702 neldat = NI_new_data_element( "ClustCounts" , itop+1 ) ;
703 for( nnn=1 ; nnn <= 3 ; nnn++ ){
704 for( ii=0 ; ii < npthr ; ii++ ){
705 NI_add_column( neldat , NI_INT , max_table[nnn][ipthr] ) ;
706 }}
707
708 sprintf(buf,"%d",niter) ;
709 NI_set_attribute( neldat , "niter" , buf ) ;
710
711 NI_write_element_tofile( fnam , neldat , NI_TEXT_MODE ) ;
712
713 NI_free_element(neldat) ; neldat = NULL ;
714
715 } /* end of save */
716
717 /*---------- compute and print the output table ----------*/
718
719 if( do_final ){
720 double *alpha , aval , ahi,alo ;
721 float **clust_thresh , cmax=0.0f ;
722 int ii , itop , iathr ;
723 char *commandline = tross_commandline("3dClustCount",argc,argv) ;
724 char fname[THD_MAX_NAME] , pname[THD_MAX_NAME] ;
725 char *amesg = NULL ; /* 20 Jan 2011 */
726
727 alpha = (double *)malloc(sizeof(double)*(max_cluster_size+1)) ;
728 clust_thresh = (float **)malloc(sizeof(float *)*npthr) ;
729 for( ipthr=0 ; ipthr < npthr ; ipthr++ )
730 clust_thresh[ipthr] = (float *)malloc(sizeof(float)*nathr) ;
731
732 for( nnn=1 ; nnn <= 3 ; nnn++ ){
733 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
734 for( itop=ii=1 ; ii <= max_cluster_size ; ii++ ){
735 alpha[ii] = max_table[nnn][ipthr][ii] / (double)niter ;
736 if( alpha[ii] > 0.0 ) itop = ii ;
737 }
738 for( ii=itop-1 ; ii >= 1 ; ii-- ) alpha[ii] += alpha[ii+1] ;
739 #if 0
740 INFO_message("pthr[%d]=%g itop=%d",ipthr,pthr[ipthr],itop) ;
741 for( ii=1 ; ii <= itop ; ii++ )
742 fprintf(stderr," %d=%g",ii,alpha[ii]) ;
743 fprintf(stderr,"\n") ;
744 #endif
745 for( iathr=0 ; iathr < nathr ; iathr++ ){
746 aval = athr[iathr] ;
747 if( aval > alpha[1] ){ /* unpleasant situation */
748 ii = 1 ;
749 amesg = THD_zzprintf(
750 amesg ,
751 " NN=%d pthr=%9.6f alpha=%6.3f [max simulated alpha=%6.3f]\n" ,
752 nnn , pthr[ipthr] , aval , alpha[1] ) ;
753 } else {
754 for( ii=1 ; ii < itop ; ii++ ){
755 if( alpha[ii] >= aval && alpha[ii+1] <= aval ) break ;
756 }
757 }
758
759 alo=alpha[ii] ; ahi=alpha[ii+1] ;
760 if( do_lohi ){
761 aval = ii ; /* for debugging */
762 } else {
763 if( alo >= 1.0 ) alo = 1.0 - 0.1/niter ;
764 if( ahi <= 0.0 ) ahi = 0.1/niter ;
765 if( ahi >= alo ) ahi = 0.1*alo ;
766 aval = log(-log(1.0-aval)) ;
767 alo = log(-log(1.0-alo)) ;
768 ahi = log(-log(1.0-ahi)) ;
769 aval = ii + (aval-alo)/(ahi-alo) ;
770 if( aval < 1.0 ) aval = 1.0 ;
771 else if( nodec ) aval = (int)(aval+0.951) ;
772 }
773 clust_thresh[ipthr][iathr] = aval ;
774
775 if( clust_thresh[ipthr][iathr] > cmax ) cmax = clust_thresh[ipthr][iathr] ;
776 }
777 }
778
779 if( do_lohi == 0 ){
780 /* edit each column to increase as pthr increases [shouldn't be needed] */
781
782 for( iathr=0 ; iathr < nathr ; iathr++ ){
783 for( ipthr=npthr-2 ; ipthr >= 0 ; ipthr-- ){
784 if( clust_thresh[ipthr][iathr] < clust_thresh[ipthr+1][iathr] )
785 clust_thresh[ipthr][iathr] = clust_thresh[ipthr+1][iathr] ;
786 }
787 }
788
789 /* edit each row to increase as athr decreases [shouldn't be needed] */
790
791 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
792 for( iathr=1 ; iathr < nathr ; iathr++ ){
793 if( clust_thresh[ipthr][iathr] < clust_thresh[ipthr][iathr-1] )
794 clust_thresh[ipthr][iathr] = clust_thresh[ipthr][iathr-1] ;
795 }
796 }
797 }
798
799 #if 0
800 if( !nodec && !do_niml && cmax > 9999.9f ){ /* if largest is way big, */
801 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){ /* then truncate to ints. */
802 for( iathr=0 ; iathr < nathr ; iathr++ ){
803 aval = clust_thresh[ipthr][iathr] ;
804 aval = (int)(aval+0.951) ;
805 clust_thresh[ipthr][iathr] = aval ;
806 }
807 }
808 nodec = 1 ;
809 }
810 #endif
811
812 MPROBE ;
813
814 if( prefix != NULL ){
815 sprintf(pname,"%s.NN%d.",prefix,nnn) ;
816 } else {
817 fflush(stderr) ; fflush(stdout) ;
818 }
819
820 if( do_1D ){ /* output in 1D format */
821 FILE *fp = stdout ;
822 if( prefix != NULL ){
823 strcpy(fname,pname) ; strcat(fname,"1D") ; fp = fopen(fname,"w") ;
824 if( fp == NULL ){
825 ERROR_message("Can't open file %s -- using stdout",fname) ;
826 fp = stdout ;
827 }
828 }
829 fprintf(fp,
830 "# %s\n"
831 "#\n"
832 "# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels\n"
833 "# -NN %d | alpha = Prob(Cluster >= given size)\n"
834 "# pthr |" ,
835 commandline , nnn ) ;
836 for( iathr=0 ; iathr < nathr ; iathr++ ) fprintf(fp," %6.3f",athr[iathr]) ;
837 fprintf(fp,"\n"
838 "# ------ |" ) ;
839 for( iathr=0 ; iathr < nathr ; iathr++ ) fprintf(fp," ------") ;
840 fprintf(fp,"\n") ;
841 for( ipthr=0 ; ipthr < npthr ; ipthr++ ){
842 fprintf(fp,"%9.6f ",pthr[ipthr]) ;
843 for( iathr=0 ; iathr < nathr ; iathr++ ){
844 if( nodec )
845 fprintf(fp,"%7d" ,(int)clust_thresh[ipthr][iathr]) ;
846 else if( clust_thresh[ipthr][iathr] <= 9999.9f )
847 fprintf(fp,"%7.1f", clust_thresh[ipthr][iathr]) ;
848 else
849 fprintf(fp,"%7.0f", clust_thresh[ipthr][iathr]) ;
850 }
851 fprintf(fp,"\n") ;
852 }
853 }
854
855 if( do_niml ){ /* output in NIML format */
856 NI_element *nel ; float *vec ; char buf[1024] , *bbb ; NI_float_array nfar ;
857 sprintf(buf,"3dClustCount_NN%d",nnn) ;
858 nel = NI_new_data_element( buf , npthr ) ;
859 vec = (float *)malloc(sizeof(float)*MAX(npthr,nathr)) ;
860 for( iathr=0 ; iathr < nathr ; iathr++ ){
861 for( ipthr=0 ; ipthr < npthr ; ipthr++ )
862 vec[ipthr] = clust_thresh[ipthr][iathr] ;
863 NI_add_column( nel , NI_FLOAT , vec ) ;
864 }
865 NI_set_attribute( nel , "commandline" , commandline ) ;
866 sprintf(buf,"%d,%d,%d",nx,ny,nz) ;
867 NI_set_attribute(nel,"nxyz",buf) ;
868 #if 0
869 sprintf(buf,"%.3f,%.3f,%.3f",dx,dy,dz) ;
870 NI_set_attribute(nel,"dxyz",buf) ;
871 sprintf(buf,"%.2f,%.2f,%.2f",fwhm_x,fwhm_y,fwhm_z) ;
872 NI_set_attribute(nel,"fwhmxyz",buf) ;
873 #endif
874 sprintf(buf,"%d",niter) ;
875 NI_set_attribute(nel,"iter",buf) ;
876 for( ipthr=0 ; ipthr < npthr ; ipthr++ ) vec[ipthr] = pthr[ipthr] ;
877 nfar.num = npthr ; nfar.ar = vec ; bbb = NI_encode_float_list(&nfar,",") ;
878 NI_set_attribute(nel,"pthr",bbb) ; NI_free(bbb) ;
879 for( iathr=0 ; iathr < nathr ; iathr++ ) vec[iathr] = athr[iathr] ;
880 nfar.num = nathr ; nfar.ar = vec ; bbb = NI_encode_float_list(&nfar,",") ;
881 NI_set_attribute(nel,"athr",bbb) ; NI_free(bbb) ;
882 #if 0
883 if( mask_dset != NULL ){
884 NI_set_attribute(nel,"mask_dset_idcode",DSET_IDCODE_STR(mask_dset)) ;
885 NI_set_attribute(nel,"mask_dset_name" ,DSET_HEADNAME(mask_dset)) ;
886 sprintf(buf,"%d",mask_ngood) ;
887 NI_set_attribute(nel,"mask_count",buf) ;
888 }
889 #endif
890 if( prefix != NULL ){ strcpy(fname,pname) ; strcat(fname,"niml") ; }
891 else strcpy(fname,"stdout:") ;
892 NI_write_element_tofile( fname , nel , NI_TEXT_MODE ) ;
893 NI_free_element( nel ) ;
894 } /* end of NIML output */
895
896 } /* end of loop over nnn = NN degree */
897
898 if( amesg != NULL ){
899 WARNING_message("Simulation not effective for these cases:\n\n"
900 "%s\n"
901 "*+ This means that not enough clusters, of any size, +*\n"
902 " of voxels at or below each pthr threshold, were +*\n"
903 " found to estimate at each alpha level. +*\n"
904 "*+ In other words, the probability that noise-only +*\n"
905 " data (of the given smoothness) will cause +*\n"
906 " above-threshold (at the given pthr) clusters is +*\n"
907 " smaller than the desired alpha levels. +*\n"
908 "*+ This problem can arise when the masked region +*\n"
909 " being simulated is small and at the same time +*\n"
910 " the smoothness (FWHM) is large. +*\n"
911 "*+ Read the 'CAUTION and CAVEAT' section at the end +*\n"
912 " of the '-help' output for a longer explanation. +*\n\n"
913 , amesg ) ;
914 free(amesg) ;
915 }
916
917 } /* end of outputizationing */
918
919 /*-------- run away screaming into the night ----- AAUUGGGHHH!!! --------*/
920
921 exit(0);
922 }
923