1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 This program performs agglomerative hierarchical clustering of voxels
10 for user specified parameter sub-bricks.
11
12
13 File: 3dStatClust.c
14 Author: B. Douglas Ward
15 Date: 08 October 1999
16
17 Mod: Replaced C "pow" function, significantly improving execution speed.
18 Date: 11 October 1999
19
20 Mod: Replaced dataset interface code with call to THD_open_dataset.
21 Restructured code for initializing hierarchical clustering.
22 Date: 19 October 1999
23
24 Mod: At each output cluster agglomeration step, print to the screen
25 which clusters are to be combined, and their distance.
26 Date: 05 September 2000
27
28 Mod: Corrected error in sort_clusters routine.
29 Date: 30 April 2001
30
31 Mod: Added call to AFNI_logger.
32 Date: 15 August 2001
33
34 */
35 /*---------------------------------------------------------------------------*/
36
37 #define PROGRAM_NAME "3dStatClust" /* name of this program */
38 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
39 #define PROGRAM_INITIAL "08 October 1999" /* date of initial program release */
40 #define PROGRAM_LATEST "15 August 2001" /* date of last program revision */
41
42 #define MAX_PARAMETERS 100
43
44 /*---------------------------------------------------------------------------*/
45 /*
46 Include header files.
47 */
48
49 #include "mrilib.h"
50 #include "matrix.h"
51
52 /*---------------------------------------------------------------------------*/
53
54 /** macro to test a malloc-ed pointer for validity **/
55
56 #define MTEST(ptr) \
57 if((ptr)==NULL) \
58 ( printf ("Cannot allocate memory \n"), exit(1) )
59
60
61 /*---------------------------------------------------------------------------*/
62
63 /** macro to open a dataset and make it ready for processing **/
64
65 #define DOPEN(ds,name) \
66 do{ int pv ; (ds) = THD_open_dataset((name)) ; \
67 if( !ISVALID_3DIM_DATASET((ds)) ){ \
68 fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
69 DSET_load((ds)) ; \
70 pv = DSET_PRINCIPAL_VALUE((ds)) ; \
71 if( DSET_ARRAY((ds),pv) == NULL ){ \
72 fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \
73 if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \
74 fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1); } \
75 break ; } while (0)
76
77
78 /*---------------------------------------------------------------------------*/
79
80 /** macro to return pointer to correct location in brick for current processing **/
81
82 #define SUB_POINTER(ds,vv,ind,ptr) \
83 do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \
84 default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \
85 case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \
86 (ptr) = (void *)( fim + (ind) ) ; \
87 } break ; \
88 case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \
89 (ptr) = (void *)( fim + (ind) ) ; \
90 } break ; \
91 case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \
92 (ptr) = (void *)( fim + (ind) ) ; \
93 } break ; } break ; } while(0)
94
95
96 /*-------------------------- global data ------------------------------------*/
97
98 static int SC_nvox = -1; /* # voxels */
99 static int SC_verb = 0; /* verbose? */
100 static float SC_thr = -1.0; /* threshold */
101 static int SC_nclust = 10; /* max. output clusters */
102 static int SC_statdist = 0; /* dist. calc. method */
103 static int SC_dimension= 0; /* number of parameters */
104
105 static char SC_thr_filename[THD_MAX_NAME] = "";
106 static char SC_output_prefix[THD_MAX_PREFIX] = "SC" ;
107 static char SC_session[THD_MAX_NAME] = "./" ;
108
109 static int * SC_voxels = NULL; /* indices for voxels above thr. */
110 static float ** SC_parameters = NULL; /* parameters for voxels above thr. */
111
112 static char * commandline = NULL ; /* command line for history notes */
113
114
115 /*---------------------------------------------------------------------------*/
116 /*
117 Print error message and stop.
118 */
119
SC_error(char * message)120 void SC_error (char * message)
121 {
122 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
123 exit(1);
124 }
125
126
127 /*---------------------------------------------------------------------------*/
128 /*
129 Include source code files.
130 */
131
132 #include "StatClust.c"
133
134
135 /*---------------------------------------------------------------------------*/
136 /*
137 Display help file.
138 */
139
SC_Syntax(void)140 void SC_Syntax(void)
141 {
142 printf(
143 "Perform agglomerative hierarchical clustering for user specified \n"
144 "parameter sub-bricks, for all voxels whose threshold statistic \n"
145 "is above a user specified value.\n"
146 "\nUsage: 3dStatClust options datasets \n"
147 "where the options are:\n"
148 ) ;
149
150 printf(
151 "-prefix pname = Use 'pname' for the output dataset prefix name.\n"
152 " OR [default='SC']\n"
153 "-output pname\n"
154 "\n"
155 "-session dir = Use 'dir' for the output dataset session directory.\n"
156 " [default='./'=current working directory]\n"
157 "-verb = Print out verbose output as the program proceeds.\n"
158 "\n"
159 "Options for calculating distance between parameter vectors: \n"
160 " -dist_euc = Calculate Euclidean distance between parameters \n"
161 " -dist_ind = Statistical distance for independent parameters \n"
162 " -dist_cor = Statistical distance for correlated parameters \n"
163 "The default option is: Euclidean distance. \n"
164 "\n"
165 "-thresh t tname = Use threshold statistic from file tname. \n"
166 " Only voxels whose threshold statistic is greater \n"
167 " than t in abolute value will be considered. \n"
168 " [If file tname contains more than 1 sub-brick, \n"
169 " the threshold stat. sub-brick must be specified!]\n"
170 "-nclust n = This specifies the maximum number of clusters for \n"
171 " output (= number of sub-bricks in output dataset).\n"
172 "\n"
173 "Command line arguments after the above are taken as parameter datasets.\n"
174 "\n"
175 ) ;
176
177 printf("\n" MASTER_SHORTHELP_STRING ) ;
178
179 PRINT_COMPILE_DATE ; exit(0) ;
180
181 }
182
183
184 /*---------------------------------------------------------------------------*/
185 /*
186 Read the arguments, load the global variables
187
188 */
189
SC_read_opts(int argc,char * argv[])190 int SC_read_opts( int argc , char * argv[] )
191 {
192 int nopt = 1 , ii ;
193 char dname[THD_MAX_NAME] ;
194 char subv[THD_MAX_NAME] ;
195 char * cpt ;
196 THD_3dim_dataset * dset ;
197 int * svar ;
198 char * str;
199 int ok, ilen, nlen , max_nsub=0 ;
200
201 char message[80]; /* error message */
202
203
204 while( nopt < argc ){
205
206 /**** -prefix prefix ****/
207
208 if( strncmp(argv[nopt],"-prefix",6) == 0 ||
209 strncmp(argv[nopt],"-output",6) == 0 ){
210 nopt++ ;
211 if( nopt >= argc ){
212 SC_error (" need argument after -prefix!");
213 }
214 MCW_strncpy( SC_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
215 continue ;
216 }
217
218 /**** -session directory ****/
219
220 if( strncmp(argv[nopt],"-session",6) == 0 ){
221 nopt++ ;
222 if( nopt >= argc ){
223 SC_error (" need argument after -session!");
224 }
225 MCW_strncpy( SC_session , argv[nopt++] , THD_MAX_NAME ) ;
226 continue ;
227 }
228
229 /**** -verb ****/
230
231 if( strncmp(argv[nopt],"-verb",5) == 0 ){
232 SC_verb++ ;
233 nopt++ ; continue ;
234 }
235
236 /**** -dist_euc ****/
237
238 if( strncmp(argv[nopt],"-dist_euc",9) == 0 ){
239 SC_statdist = 0 ;
240 nopt++ ; continue ;
241 }
242
243 /**** -dist_ind ****/
244
245 if( strncmp(argv[nopt],"-dist_ind",9) == 0 ){
246 SC_statdist = 1 ;
247 nopt++ ; continue ;
248 }
249
250 /**** -dist_cor ****/
251
252 if( strncmp(argv[nopt],"-dist_cor",9) == 0 ){
253 SC_statdist = 2 ;
254 nopt++ ; continue ;
255 }
256
257 /**** -nclust n ****/
258
259 if( strncmp(argv[nopt],"-nclust",7) == 0 ){
260 int ival;
261 nopt++ ;
262 if( nopt >= argc ){
263 SC_error (" need argument after -nclust!");
264 }
265 sscanf (argv[nopt], "%d", &ival);
266 if ((ival < 1) || (ival > 255)){
267 SC_error (" Require 1 <= nclust <= 255 ");
268 }
269 SC_nclust = ival;
270 nopt++;
271 continue;
272 }
273
274
275 /**** -thresh thr fname ****/
276
277 if( strncmp(argv[nopt],"-thresh",7) == 0 ){
278 float fval;
279 nopt++ ;
280 if( nopt+1 >= argc ){
281 SC_error (" need 2 arguments after -thresh!");
282 }
283 sscanf (argv[nopt], "%f", &fval);
284 if (fval < 0.0){
285 SC_error (" Require thr >= 0.0 ");
286 }
287 SC_thr = fval;
288 nopt++;
289
290 strcpy (SC_thr_filename, argv[nopt]);
291 nopt++;
292 continue;
293 }
294
295
296 /*----- Invalid option -----*/
297 if( argv[nopt][0] == '-' ){
298 sprintf (message, " Unknown option: %s ", argv[nopt]);
299 SC_error (message);
300 }
301
302
303 /*----- Remaining inputs should be parameter datasets -----*/
304 break;
305
306
307 } /* end of loop over command line arguments */
308
309
310 return (nopt) ;
311 }
312
313
314 /*---------------------------------------------------------------------------*/
315 /*
316 Routine to initialize the program: get all operator inputs; get indices
317 for voxels above threshold; get parameter vectors for all voxels above
318 threshold.
319 */
320
initialize_program(int argc,char * argv[])321 THD_3dim_dataset * initialize_program ( int argc , char * argv[] )
322 {
323 const int MIN_NVOX = 10; /* minimum number of voxels above threshold */
324
325 THD_3dim_dataset * thr_dset=NULL; /* threshold dataset */
326 THD_3dim_dataset * param_dset=NULL; /* parameter dataset(s) */
327
328 int nx, ny, nz; /* dataset dimensions in voxels */
329 int iv; /* index number of sub-brick */
330 void * vfim = NULL; /* sub-brick data pointer */
331 float * ffim = NULL; /* sub-brick data in floating point format */
332 int ivox, nvox, icount; /* voxel indices */
333 int nopt; /* points to current input option */
334 int ibrick, nbricks; /* sub-brick indices */
335 char message[80]; /* error message */
336
337
338 /*----- Save command line for history notes -----*/
339 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
340
341
342 /*----- Does user request help menu? -----*/
343 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) SC_Syntax() ;
344
345 /*----- Add to program log -----*/
346 AFNI_logger (PROGRAM_NAME,argc,argv);
347
348
349 /*----- Read input options -----*/
350 nopt = SC_read_opts( argc , argv ) ;
351
352
353 /*----- Open the threshold dataset -----*/
354 if (SC_verb) printf ("Reading threshold dataset: %s \n", SC_thr_filename);
355 DOPEN (thr_dset, SC_thr_filename);
356
357 if (thr_dset == NULL)
358 {
359 sprintf (message, "Cannot open threshold dataset %s", SC_thr_filename);
360 SC_error (message);
361 }
362
363 if (DSET_NVALS(thr_dset) != 1)
364 SC_error ("Must specify single sub-brick for threshold data");
365
366
367 /*----- Save dimensions of threshold dataset for compatibility test -----*/
368 nx = DSET_NX(thr_dset); ny = DSET_NY(thr_dset); nz = DSET_NZ(thr_dset);
369
370
371 /*----- Allocate memory for float data -----*/
372 nvox = DSET_NVOX (thr_dset);
373 ffim = (float *) malloc (sizeof(float) * nvox); MTEST (ffim);
374
375
376 /*----- Convert threshold dataset sub-brick to floats (in ffim) -----*/
377 iv = DSET_PRINCIPAL_VALUE (thr_dset);
378 SUB_POINTER (thr_dset, iv, 0, vfim);
379 EDIT_coerce_scale_type (nvox, DSET_BRICK_FACTOR(thr_dset,iv),
380 DSET_BRICK_TYPE(thr_dset,iv), vfim, /* input */
381 MRI_float , ffim); /* output */
382
383 /*----- Delete threshold dataset -----*/
384 THD_delete_3dim_dataset (thr_dset, False); thr_dset = NULL ;
385
386
387 /*----- Count number of voxels above threshold -----*/
388 SC_nvox = 0;
389 for (ivox = 0; ivox < nvox; ivox++)
390 if (fabs(ffim[ivox]) > SC_thr) SC_nvox++;
391 if (SC_verb) printf ("Number of voxels above threshold = %d \n", SC_nvox);
392 if (SC_nvox < MIN_NVOX)
393 {
394 sprintf (message, "Only %d voxels above threshold. Cannot continue.",
395 SC_nvox);
396 SC_error (message);
397 }
398
399
400
401 /*----- Allocate memory for voxel index array -----*/
402 SC_voxels = (int *) malloc (sizeof(int) * SC_nvox);
403 MTEST (SC_voxels);
404
405
406 /*----- Save indices of voxels above threshold -----*/
407 icount = 0;
408 for (ivox = 0; ivox < nvox; ivox++)
409 if (fabs(ffim[ivox]) > SC_thr)
410 {
411 SC_voxels[icount] = ivox;
412 icount++;
413 }
414
415
416 /*----- Allocate memory for parameter array -----*/
417 SC_parameters = (float **) malloc (sizeof(float *) * MAX_PARAMETERS);
418 MTEST (SC_parameters);
419
420
421 /*----- Begin loop over parameter datasets -----*/
422 SC_dimension = 0;
423 while (nopt < argc)
424 {
425 /*----- Check if this is an input option -----*/
426 if (argv[nopt][0] == '-')
427 SC_error ("ALL input options must precede ALL parameter datasets");
428
429 /*----- Open the parameter dataset -----*/
430 if (SC_verb) printf ("Reading parameter dataset: %s \n", argv[nopt]);
431 DOPEN (param_dset, argv[nopt]);
432
433 if (param_dset == NULL)
434 {
435 sprintf (message, "Cannot open parameter dataset %s", argv[nopt]);
436 SC_error (message);
437 }
438
439 /*----- Test for dataset compatibility -----*/
440 if ((nx != DSET_NX(param_dset)) || (ny != DSET_NY(param_dset))
441 || (nz != DSET_NZ(param_dset)))
442 {
443 sprintf (message, "Parameter dataset %s has incompatible dimensions",
444 argv[nopt]);
445 SC_error (message);
446 }
447
448
449 /*----- Get number of parameters specified by this dataset -----*/
450 nbricks = DSET_NVALS(param_dset);
451
452
453 /*----- Loop over sub-bricks selected from parameter dataset -----*/
454 for (ibrick = 0; ibrick < nbricks; ibrick++)
455 {
456 if (SC_verb) printf ("Reading parameter #%2d \n", SC_dimension+1);
457
458 SUB_POINTER (param_dset, ibrick, 0, vfim);
459 EDIT_coerce_scale_type (nvox, DSET_BRICK_FACTOR(param_dset,ibrick),
460 DSET_BRICK_TYPE(param_dset,ibrick), vfim, /* input */
461 MRI_float , ffim); /* output */
462
463 /*----- Allocate memory for parameter data -----*/
464 SC_parameters[SC_dimension]
465 = (float *) malloc (sizeof(float) * SC_nvox);
466 MTEST (SC_parameters[SC_dimension]);
467
468 /*----- Save parameter data for all voxels above threshold -----*/
469 for (ivox = 0; ivox < SC_nvox; ivox++)
470 SC_parameters[SC_dimension][ivox] = ffim[SC_voxels[ivox]];
471
472 /*----- Increment count of parameters -----*/
473 SC_dimension++;
474
475 }
476
477 /*----- Delete parameter dataset -----*/
478 THD_delete_3dim_dataset (param_dset, False); param_dset = NULL ;
479
480 nopt++;
481 }
482
483
484 /*----- Delete floating point sub-brick -----*/
485 if (ffim != NULL) { free (ffim); ffim = NULL; }
486
487
488 if (SC_verb) printf ("Number of parameters = %d \n", SC_dimension);
489 if (SC_dimension < 1) SC_error ("No parameter data?");
490
491
492 return NULL ;
493 }
494
495
496 /*---------------------------------------------------------------------------*/
497 /*
498 Perform agglomerative hierarchical clustering.
499 */
500
form_clusters()501 THD_3dim_dataset * form_clusters ()
502
503 {
504 THD_3dim_dataset * new_dset = NULL; /* hierarchical clustering */
505 THD_3dim_dataset * thr_dset = NULL; /* threshold dataset */
506 int ivox, ixyz, nxyz; /* voxel indices */
507 int iclust; /* cluster index */
508 int ip, jp; /* parameter indices */
509 cluster * head_clust = NULL; /* last cluster */
510 float * parameters = NULL; /* parameters after normalization */
511 byte ** bar = NULL; /* array of cluster sub-bricks */
512 int nbricks; /* number of cluster sub-bricks */
513 int ibrick; /* cluster sub-brick index */
514 int ierror; /* number of errors in editing data */
515 int ok; /* Boolean for successful matrix calc. */
516
517 vector v, av; /* intermediate vector results */
518 matrix s; /* square root of covariance matrix */
519 matrix sinv; /* inverse of square root of covariance matrix */
520
521 char message[80]; /* error message */
522
523
524 /*----- Initialize vectors and matrices -----*/
525 vector_initialize (&v);
526 vector_initialize (&av);
527 matrix_initialize (&s);
528 matrix_initialize (&sinv);
529
530
531 /*----- Calculate covariance matrix for input parameters -----*/
532 if (SC_statdist) calc_covariance (&s, &sinv);
533 else
534 {
535 matrix_identity (SC_dimension, &s);
536 matrix_identity (SC_dimension, &sinv);
537 }
538
539
540 /*----- Set number of sub-bricks -----*/
541 if (SC_nvox < SC_nclust)
542 nbricks = SC_nvox;
543 else
544 nbricks = SC_nclust;
545 if (SC_verb) printf ("Output dataset will have %d sub-bricks\n\n", nbricks);
546
547
548 /*----- Open threshold dataset -----*/
549 thr_dset = THD_open_dataset (SC_thr_filename);
550 nxyz = DSET_NVOX (thr_dset);
551
552
553 /*-- Make an empty copy of threshold dataset, for eventual output --*/
554 new_dset = EDIT_empty_copy (thr_dset);
555
556
557 /*----- Record history of dataset -----*/
558 tross_Copy_History (thr_dset, new_dset);
559 if( commandline != NULL ) tross_Append_History( new_dset , commandline ) ;
560
561
562 /*----- Delete threshold dataset -----*/
563 THD_delete_3dim_dataset (thr_dset, False); thr_dset = NULL ;
564
565
566 /*----- Modify some structural properties. Note that the nbricks
567 just make empty sub-bricks, without any data attached. -----*/
568 ierror = EDIT_dset_items (new_dset,
569 ADN_prefix, SC_output_prefix,
570 ADN_directory_name, SC_session,
571 ADN_type, HEAD_FUNC_TYPE,
572 ADN_func_type, FUNC_BUCK_TYPE,
573 ADN_ntt, 0, /* no time */
574 ADN_nvals, nbricks,
575 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
576 ADN_none ) ;
577
578 if( ierror > 0 )
579 {
580 sprintf (message,
581 " %d errors in attempting to create bucket dataset! ",
582 ierror);
583 SC_error (message);
584 }
585
586 if ( THD_deathcon() && THD_is_file(DSET_HEADNAME(new_dset)))
587 {
588 sprintf (message,
589 " Output dataset file %s already exists--cannot continue! ",
590 DSET_HEADNAME(new_dset));
591 SC_error (message);
592 }
593
594
595 /*----- Allocate memory -----*/
596 bar = (byte **) malloc (sizeof(byte *) * nbricks);
597 MTEST (bar);
598
599
600 /*----- Build lowest level of cluster hierarchy -----*/
601 for (ivox = 0; ivox < SC_nvox; ivox++)
602 {
603 /*----- Allocate space for parameter vector -----*/
604 parameters = (float *) malloc (sizeof(float) * SC_dimension);
605 MTEST (parameters);
606
607 /*----- Copy the parameter array -----*/
608 for (ip = 0; ip < SC_dimension; ip++)
609 parameters[ip] = SC_parameters[ip][ivox];
610
611 /*----- If using stat. dist., transform the parameter vector -----*/
612 if (SC_statdist)
613 {
614 array_to_vector (SC_dimension, parameters, &v);
615 vector_multiply (sinv, v, &av);
616 vector_to_array (av, parameters);
617 }
618
619 /*----- Create new cluster containing single voxel -----*/
620 ixyz = SC_voxels[ivox];
621 head_clust = new_cluster (ixyz, parameters, head_clust);
622 }
623
624
625 /*----- Deallocate memory for parameter data -----*/
626 free (SC_voxels); SC_voxels = NULL;
627 for (ip = 0; ip < SC_dimension; ip++)
628 {
629 free (SC_parameters[ip]); SC_parameters[ip] = NULL;
630 }
631 free (SC_parameters); SC_parameters = NULL;
632
633
634 /*----- Agglomerate clusters, one-by-one -----*/
635 for (iclust = SC_nvox; iclust > 0; iclust--)
636 {
637 if (SC_verb && (iclust % 100 == 0))
638 printf ("# Clusters = %d \n", iclust);
639
640 if (iclust <= nbricks)
641 {
642 /*----- Sort clusters in order of size -----*/
643 head_clust = sort_clusters (head_clust);
644
645 /*----- Print cluster centroid parameters -----*/
646 if (SC_verb)
647 {
648 printf ("\n\n# Clusters = %d \n\n", iclust);
649 print_all_clusters (head_clust, s);
650 }
651
652 /*----- allocate memory for output sub-brick -----*/
653 ibrick = iclust-1;
654 bar[ibrick] = (byte *) malloc (sizeof(byte) * nxyz);
655 MTEST (bar[ibrick]);
656
657 /*----- Save clusters into output sub-brick -----*/
658 for (ixyz = 0; ixyz < nxyz; ixyz++)
659 bar[ibrick][ixyz] = 0;
660 save_all_clusters (head_clust, bar[ibrick]);
661
662 /*----- attach bar[ib] to be sub-brick #ibrick -----*/
663 EDIT_substitute_brick (new_dset, ibrick, MRI_byte, bar[ibrick]);
664
665 }
666
667 /*----- Agglomerate clusters -----*/
668 if (iclust > 1)
669 head_clust = agglomerate_clusters (head_clust,
670 SC_verb*(iclust <= nbricks));
671
672
673 }
674
675
676 /*----- Deallocate memory -----*/
677 vector_destroy (&v);
678 vector_destroy (&av);
679 matrix_destroy (&s);
680 matrix_destroy (&sinv);
681 destroy_cluster (head_clust);
682
683
684 /*----- Return hierarchical clustering -----*/
685 return (new_dset);
686
687 }
688
689
690 /*---------------------------------------------------------------------------*/
691
main(int argc,char * argv[])692 int main( int argc , char * argv[] )
693
694 {
695 THD_3dim_dataset * clust_dset = NULL; /* hierarchical clusters data set */
696 int ip; /* parameter index */
697
698
699 /*----- Identify software -----*/
700 #if 0
701 printf ("\n\n");
702 printf ("Program: %s \n", PROGRAM_NAME);
703 printf ("Author: %s \n", PROGRAM_AUTHOR);
704 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
705 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
706 printf ("\n");
707 #endif
708
709 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
710
711 PRINT_VERSION("3dStatClust") ; AUTHOR(PROGRAM_AUTHOR);
712 mainENTRY("3dStatClust main") ; machdep() ;
713
714 { int new_argc ; char ** new_argv ;
715 addto_args( argc , argv , &new_argc , &new_argv ) ;
716 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
717 }
718
719
720 /*----- Initialize program: get all operator inputs; get indices
721 for voxels above threshold; get parameter vectors for all voxels
722 above threshold -----*/
723 initialize_program (argc, argv);
724
725
726 /*----- Perform agglomerative hierarchical clustering -----*/
727 clust_dset = form_clusters ();
728
729
730 /*----- Output the hierarchical clustering dataset -----*/
731 if( SC_verb ) fprintf(stderr,"Computing sub-brick statistics\n") ;
732 THD_load_statistics( clust_dset ) ;
733
734 THD_write_3dim_dataset( NULL,NULL , clust_dset , True ) ;
735 if( SC_verb ) fprintf(stderr,"Wrote output to %s\n", DSET_BRIKNAME(clust_dset) );
736
737
738 /*----- Deallocate memory for cluster dataset -----*/
739 THD_delete_3dim_dataset( clust_dset , False ) ; clust_dset = NULL ;
740
741 exit(0) ;
742
743
744 }
745
746
747
748
749
750
751
752