00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This program performs agglomerative hierarchical clustering of voxels 
00010   for user specified parameter sub-bricks.
00011 
00012 
00013   File:    3dStatClust.c
00014   Author:  B. Douglas Ward
00015   Date:    08 October 1999
00016 
00017   Mod:     Replaced C "pow" function, significantly improving execution speed.
00018   Date:    11 October 1999
00019 
00020   Mod:     Replaced dataset interface code with call to THD_open_dataset.
00021            Restructured code for initializing hierarchical clustering.
00022   Date:    19 October 1999
00023 
00024   Mod:     At each output cluster agglomeration step, print to the screen
00025            which clusters are to be combined, and their distance.
00026   Date:    05 September 2000
00027 
00028   Mod:     Corrected error in sort_clusters routine.
00029   Date:    30 April 2001
00030 
00031   Mod:     Added call to AFNI_logger.
00032   Date:    15 August 2001
00033 
00034 */
00035 /*---------------------------------------------------------------------------*/
00036 
00037 #define PROGRAM_NAME "3dStatClust"                   /* name of this program */
00038 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00039 #define PROGRAM_INITIAL "08 October 1999" /* date of initial program release */
00040 #define PROGRAM_LATEST "15 August 2001"     /* date of last program revision */
00041 
00042 #define MAX_PARAMETERS 100
00043 
00044 /*---------------------------------------------------------------------------*/
00045 /*
00046   Include header files.
00047 */
00048 
00049 #include "mrilib.h"
00050 #include "matrix.h"
00051 
00052 
00053 /*---------------------------------------------------------------------------*/
00054 
00055 #ifndef myXtFree
00056 #   define myXtFree(xp) (XtFree((char *)(xp)) , (xp)=NULL)
00057 #endif
00058 
00059 /*---------------------------------------------------------------------------*/
00060 
00061 /** macro to test a malloc-ed pointer for validity **/
00062 
00063 #define MTEST(ptr) \
00064 if((ptr)==NULL) \
00065 ( printf ("Cannot allocate memory \n"),  exit(1) )
00066      
00067 
00068 /*---------------------------------------------------------------------------*/
00069 
00070 /** macro to open a dataset and make it ready for processing **/
00071 
00072 #define DOPEN(ds,name)                                                        \
00073 do{ int pv ; (ds) = THD_open_dataset((name)) ;                                \
00074        if( !ISVALID_3DIM_DATASET((ds)) ){                                     \
00075           fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00076        THD_load_datablock( (ds)->dblk ) ;                                     \
00077        pv = DSET_PRINCIPAL_VALUE((ds)) ;                                      \
00078        if( DSET_ARRAY((ds),pv) == NULL ){                                     \
00079           fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }   \
00080        if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                         \
00081           fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);        }     \
00082        break ; } while (0)
00083 
00084 
00085 /*---------------------------------------------------------------------------*/
00086 
00087 /** macro to return pointer to correct location in brick for current processing **/
00088 
00089 #define SUB_POINTER(ds,vv,ind,ptr)                                            \
00090    do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
00091          default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
00092             case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
00093                             (ptr) = (void *)( fim + (ind) ) ;                 \
00094             } break ;                                                         \
00095             case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
00096                             (ptr) = (void *)( fim + (ind) ) ;                 \
00097             } break ;                                                         \
00098             case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
00099                              (ptr) = (void *)( fim + (ind) ) ;                \
00100             } break ; } break ; } while(0)
00101 
00102 
00103 /*-------------------------- global data ------------------------------------*/
00104 
00105 static int                      SC_nvox     = -1;   /* # voxels */
00106 static int                      SC_verb     = 0;    /* verbose? */
00107 static float                    SC_thr      = -1.0; /* threshold */
00108 static int                      SC_nclust   = 10;   /* max. output clusters */
00109 static int                      SC_statdist = 0;    /* dist. calc. method */
00110 static int                      SC_dimension= 0;    /* number of parameters */
00111 
00112 static char SC_thr_filename[THD_MAX_NAME]    = "";
00113 static char SC_output_prefix[THD_MAX_PREFIX] = "SC" ;
00114 static char SC_session[THD_MAX_NAME]         = "./"   ;
00115 
00116 static int * SC_voxels = NULL;           /* indices for voxels above thr. */
00117 static float ** SC_parameters = NULL;    /* parameters for voxels above thr. */
00118 
00119 static char * commandline = NULL ;       /* command line for history notes */
00120 
00121 
00122 /*---------------------------------------------------------------------------*/
00123 /*
00124    Print error message and stop.
00125 */
00126 
00127 void SC_error (char * message)
00128 {
00129   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00130   exit(1);
00131 }
00132 
00133 
00134 /*---------------------------------------------------------------------------*/
00135 /*
00136   Include source code files.
00137 */
00138 
00139 #include "matrix.c"
00140 #include "StatClust.c"
00141 
00142 
00143 /*---------------------------------------------------------------------------*/
00144 /*
00145   Display help file.
00146 */
00147 
00148 void SC_Syntax(void)
00149 {
00150    printf(
00151     "Perform agglomerative hierarchical clustering for user specified \n"
00152     "parameter sub-bricks, for all voxels whose threshold statistic   \n"
00153     "is above a user specified value.\n"
00154     "\nUsage: 3dStatClust options datasets \n"
00155     "where the options are:\n"
00156    ) ;
00157 
00158    printf(
00159     "-prefix pname    = Use 'pname' for the output dataset prefix name.\n"
00160     "  OR                 [default='SC']\n"
00161     "-output pname\n"
00162     "\n"
00163     "-session dir     = Use 'dir' for the output dataset session directory.\n"
00164     "                     [default='./'=current working directory]\n"
00165     "-verb            = Print out verbose output as the program proceeds.\n"
00166     "\n"
00167     "Options for calculating distance between parameter vectors: \n"
00168     "   -dist_euc        = Calculate Euclidean distance between parameters \n"
00169     "   -dist_ind        = Statistical distance for independent parameters \n"
00170     "   -dist_cor        = Statistical distance for correlated parameters \n"
00171     "The default option is:  Euclidean distance. \n"
00172     "\n"
00173     "-thresh t tname  = Use threshold statistic from file tname. \n"
00174     "                   Only voxels whose threshold statistic is greater \n"
00175     "                   than t in abolute value will be considered. \n"
00176     "                     [If file tname contains more than 1 sub-brick, \n"
00177     "                     the threshold stat. sub-brick must be specified!]\n"
00178     "-nclust n        = This specifies the maximum number of clusters for \n"
00179     "                   output (= number of sub-bricks in output dataset).\n"
00180     "\n"
00181     "Command line arguments after the above are taken as parameter datasets.\n"
00182     "\n"
00183    ) ;
00184 
00185    printf("\n" MASTER_SHORTHELP_STRING ) ;
00186 
00187    exit(0) ;
00188 
00189 }
00190 
00191 
00192 /*---------------------------------------------------------------------------*/
00193 /*
00194    Read the arguments, load the global variables
00195 
00196 */
00197 
00198 int SC_read_opts( int argc , char * argv[] )
00199 {
00200    int nopt = 1 , ii ;
00201    char dname[THD_MAX_NAME] ;
00202    char subv[THD_MAX_NAME] ;
00203    char * cpt ;
00204    THD_3dim_dataset * dset ;
00205    int * svar ;
00206    char * str;
00207    int ok, ilen, nlen , max_nsub=0 ;
00208 
00209    char message[80];        /* error message */
00210 
00211 
00212    while( nopt < argc ){
00213 
00214       /**** -prefix prefix ****/
00215 
00216       if( strncmp(argv[nopt],"-prefix",6) == 0 ||
00217           strncmp(argv[nopt],"-output",6) == 0   ){
00218          nopt++ ;
00219          if( nopt >= argc ){
00220             SC_error (" need argument after -prefix!");
00221          }
00222          MCW_strncpy( SC_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00223          continue ;
00224       }
00225 
00226       /**** -session directory ****/
00227 
00228       if( strncmp(argv[nopt],"-session",6) == 0 ){
00229          nopt++ ;
00230          if( nopt >= argc ){
00231             SC_error (" need argument after -session!"); 
00232          }
00233          MCW_strncpy( SC_session , argv[nopt++] , THD_MAX_NAME ) ;
00234          continue ;
00235       }
00236 
00237       /**** -verb ****/
00238 
00239       if( strncmp(argv[nopt],"-verb",5) == 0 ){
00240          SC_verb++ ;
00241          nopt++ ; continue ;
00242       }
00243 
00244       /**** -dist_euc ****/
00245 
00246       if( strncmp(argv[nopt],"-dist_euc",9) == 0 ){
00247          SC_statdist = 0 ;
00248          nopt++ ; continue ;
00249       }
00250 
00251       /**** -dist_ind ****/
00252 
00253       if( strncmp(argv[nopt],"-dist_ind",9) == 0 ){
00254          SC_statdist = 1 ;
00255          nopt++ ; continue ;
00256       }
00257 
00258       /**** -dist_cor ****/
00259 
00260       if( strncmp(argv[nopt],"-dist_cor",9) == 0 ){
00261          SC_statdist = 2 ;
00262          nopt++ ; continue ;
00263       }
00264 
00265       /**** -nclust n ****/
00266 
00267       if( strncmp(argv[nopt],"-nclust",7) == 0 ){
00268          int ival;
00269          nopt++ ;
00270          if( nopt >= argc ){
00271             SC_error (" need argument after -nclust!");
00272          }
00273          sscanf (argv[nopt], "%d", &ival); 
00274          if ((ival < 1) || (ival > 255)){
00275             SC_error (" Require 1 <= nclust <= 255 ");
00276          }
00277          SC_nclust = ival;
00278          nopt++;
00279          continue;
00280       }
00281 
00282 
00283       /**** -thresh thr fname ****/
00284 
00285       if( strncmp(argv[nopt],"-thresh",7) == 0 ){
00286          float fval;
00287          nopt++ ;
00288          if( nopt+1 >= argc ){
00289             SC_error (" need 2 arguments after -thresh!"); 
00290          }
00291          sscanf (argv[nopt], "%f", &fval); 
00292          if (fval < 0.0){
00293             SC_error (" Require thr >= 0.0 ");
00294          }
00295          SC_thr = fval;
00296          nopt++;
00297 
00298          strcpy (SC_thr_filename, argv[nopt]);
00299          nopt++;
00300          continue;
00301       }
00302 
00303 
00304       /*----- Invalid option -----*/
00305       if( argv[nopt][0] == '-' ){
00306         sprintf (message, " Unknown option: %s ", argv[nopt]);
00307         SC_error (message);
00308       }
00309 
00310 
00311       /*----- Remaining inputs should be parameter datasets -----*/
00312       break;
00313 
00314 
00315    }  /* end of loop over command line arguments */
00316 
00317 
00318    return (nopt) ;
00319 }
00320 
00321 
00322 /*---------------------------------------------------------------------------*/
00323 /*
00324   Routine to initialize the program: get all operator inputs; get indices
00325   for voxels above threshold; get parameter vectors for all voxels above
00326   threshold.
00327 */
00328 
00329 THD_3dim_dataset * initialize_program ( int argc , char * argv[] )
00330 {
00331   const int MIN_NVOX = 10;    /* minimum number of voxels above threshold */
00332 
00333   THD_3dim_dataset * thr_dset=NULL;     /* threshold dataset */
00334   THD_3dim_dataset * param_dset=NULL;   /* parameter dataset(s) */
00335 
00336   int nx, ny, nz;          /* dataset dimensions in voxels */
00337   int iv;                  /* index number of sub-brick */
00338   void * vfim = NULL;      /* sub-brick data pointer */
00339   float * ffim = NULL;     /* sub-brick data in floating point format */
00340   int ivox, nvox, icount;  /* voxel indices */
00341   int nopt;                /* points to current input option */
00342   int ibrick, nbricks;     /* sub-brick indices */
00343   char message[80];        /* error message */
00344 
00345 
00346   /*----- Save command line for history notes -----*/
00347   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00348 
00349 
00350   /*----- Does user request help menu? -----*/
00351   if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) SC_Syntax() ;
00352 
00353   /*----- Add to program log -----*/
00354   AFNI_logger (PROGRAM_NAME,argc,argv); 
00355 
00356 
00357   /*----- Read input options -----*/
00358   nopt = SC_read_opts( argc , argv ) ;
00359 
00360 
00361   /*----- Open the threshold dataset -----*/
00362   if (SC_verb)  printf ("Reading threshold dataset: %s \n", SC_thr_filename);
00363   DOPEN (thr_dset, SC_thr_filename);
00364 
00365   if (thr_dset == NULL)
00366     {
00367       sprintf (message, "Cannot open threshold dataset %s", SC_thr_filename); 
00368       SC_error (message);
00369     }
00370 
00371   if (DSET_NVALS(thr_dset) != 1)
00372     SC_error ("Must specify single sub-brick for threshold data");
00373 
00374 
00375   /*----- Save dimensions of threshold dataset for compatibility test -----*/
00376   nx = DSET_NX(thr_dset);   ny = DSET_NY(thr_dset);   nz = DSET_NZ(thr_dset);
00377 
00378 
00379   /*----- Allocate memory for float data -----*/
00380   nvox = DSET_NVOX (thr_dset);
00381   ffim = (float *) malloc (sizeof(float) * nvox);   MTEST (ffim);
00382 
00383 
00384   /*----- Convert threshold dataset sub-brick to floats (in ffim) -----*/
00385   iv = DSET_PRINCIPAL_VALUE (thr_dset);
00386   SUB_POINTER (thr_dset, iv, 0, vfim);
00387   EDIT_coerce_scale_type (nvox, DSET_BRICK_FACTOR(thr_dset,iv),
00388                           DSET_BRICK_TYPE(thr_dset,iv), vfim,     /* input  */
00389                           MRI_float                   , ffim);    /* output */
00390   
00391   /*----- Delete threshold dataset -----*/
00392   THD_delete_3dim_dataset (thr_dset, False);  thr_dset = NULL ;
00393 
00394 
00395   /*----- Count number of voxels above threshold -----*/
00396   SC_nvox = 0;
00397   for (ivox = 0;  ivox < nvox;  ivox++)
00398     if (fabs(ffim[ivox]) > SC_thr)  SC_nvox++;
00399   if (SC_verb)  printf ("Number of voxels above threshold = %d \n", SC_nvox);
00400   if (SC_nvox < MIN_NVOX)  
00401     {
00402       sprintf (message, "Only %d voxels above threshold.  Cannot continue.",
00403                SC_nvox);
00404       SC_error (message);
00405     }
00406 
00407 
00408 
00409   /*----- Allocate memory for voxel index array -----*/
00410   SC_voxels = (int *) malloc (sizeof(int) * SC_nvox);
00411   MTEST (SC_voxels);
00412 
00413 
00414   /*----- Save indices of voxels above threshold -----*/
00415   icount = 0;
00416   for (ivox = 0;  ivox < nvox;  ivox++)
00417     if (fabs(ffim[ivox]) > SC_thr)
00418       {
00419         SC_voxels[icount] = ivox;
00420         icount++;
00421       }
00422 
00423   
00424   /*----- Allocate memory for parameter array -----*/
00425   SC_parameters = (float **) malloc (sizeof(float *) * MAX_PARAMETERS);
00426   MTEST (SC_parameters);
00427 
00428 
00429   /*----- Begin loop over parameter datasets -----*/
00430   SC_dimension = 0;
00431   while (nopt < argc)
00432     {
00433       /*----- Check if this is an input option -----*/
00434       if (argv[nopt][0] == '-')
00435         SC_error ("ALL input options must precede ALL parameter datasets");
00436 
00437       /*----- Open the parameter dataset -----*/
00438       if (SC_verb)  printf ("Reading parameter dataset: %s \n", argv[nopt]);
00439       DOPEN (param_dset, argv[nopt]);
00440 
00441       if (param_dset == NULL)
00442         {
00443           sprintf (message, "Cannot open parameter dataset %s", argv[nopt]); 
00444           SC_error (message);
00445         }
00446 
00447       /*----- Test for dataset compatibility -----*/
00448       if ((nx != DSET_NX(param_dset)) || (ny != DSET_NY(param_dset))
00449           || (nz != DSET_NZ(param_dset)))
00450         {
00451           sprintf (message, "Parameter dataset %s has incompatible dimensions",
00452                    argv[nopt]); 
00453           SC_error (message);
00454         }
00455         
00456      
00457       /*----- Get number of parameters specified by this dataset -----*/
00458       nbricks = DSET_NVALS(param_dset);
00459 
00460 
00461       /*----- Loop over sub-bricks selected from parameter dataset -----*/
00462       for (ibrick = 0;  ibrick < nbricks;  ibrick++)
00463         {
00464           if (SC_verb)  printf ("Reading parameter #%2d \n", SC_dimension+1);
00465 
00466           SUB_POINTER (param_dset, ibrick, 0, vfim);
00467           EDIT_coerce_scale_type (nvox, DSET_BRICK_FACTOR(param_dset,ibrick),
00468                      DSET_BRICK_TYPE(param_dset,ibrick), vfim,   /* input  */
00469                      MRI_float                         , ffim);  /* output */
00470   
00471           /*----- Allocate memory for parameter data -----*/
00472           SC_parameters[SC_dimension] 
00473             = (float *) malloc (sizeof(float) * SC_nvox);
00474           MTEST (SC_parameters[SC_dimension]);
00475           
00476           /*----- Save parameter data for all voxels above threshold -----*/
00477           for (ivox = 0;  ivox < SC_nvox;  ivox++)
00478             SC_parameters[SC_dimension][ivox] = ffim[SC_voxels[ivox]];
00479           
00480           /*----- Increment count of parameters -----*/
00481           SC_dimension++;
00482 
00483         }
00484 
00485       /*----- Delete parameter dataset -----*/
00486       THD_delete_3dim_dataset (param_dset, False);  param_dset = NULL ;
00487      
00488       nopt++;
00489     }
00490 
00491 
00492   /*----- Delete floating point sub-brick -----*/
00493   if (ffim != NULL) { free (ffim);   ffim = NULL; }
00494 
00495 
00496   if (SC_verb)  printf ("Number of parameters = %d \n", SC_dimension);
00497   if (SC_dimension < 1)  SC_error ("No parameter data?");
00498 
00499 
00500 }
00501 
00502 
00503 /*---------------------------------------------------------------------------*/
00504 /*
00505   Perform agglomerative hierarchical clustering.
00506 */
00507 
00508 THD_3dim_dataset * form_clusters ()
00509 
00510 {
00511   THD_3dim_dataset * new_dset = NULL;   /* hierarchical clustering */
00512   THD_3dim_dataset * thr_dset = NULL;   /* threshold dataset */
00513   int ivox, ixyz, nxyz;            /* voxel indices */
00514   int iclust;                      /* cluster index */
00515   int ip, jp;                      /* parameter indices */
00516   cluster * head_clust = NULL;     /* last cluster */
00517   float * parameters = NULL;       /* parameters after normalization */
00518   byte ** bar = NULL;              /* array of cluster sub-bricks */
00519   int nbricks;                     /* number of cluster sub-bricks */
00520   int ibrick;                      /* cluster sub-brick index */
00521   int ierror;                      /* number of errors in editing data */
00522   int ok;                          /* Boolean for successful matrix calc. */
00523 
00524   vector v, av;               /* intermediate vector results */
00525   matrix s;                   /* square root of covariance matrix */
00526   matrix sinv;                /* inverse of square root of covariance matrix */
00527 
00528   char message[80];           /* error message */
00529 
00530 
00531   /*----- Initialize vectors and matrices -----*/
00532   vector_initialize (&v);
00533   vector_initialize (&av);
00534   matrix_initialize (&s);
00535   matrix_initialize (&sinv);
00536 
00537 
00538   /*----- Calculate covariance matrix for input parameters -----*/
00539   if (SC_statdist)  calc_covariance (&s, &sinv);
00540   else
00541     {
00542       matrix_identity (SC_dimension, &s);
00543       matrix_identity (SC_dimension, &sinv);
00544     }
00545   
00546 
00547   /*----- Set number of sub-bricks -----*/
00548   if (SC_nvox < SC_nclust)
00549     nbricks = SC_nvox;
00550   else
00551     nbricks = SC_nclust;
00552   if (SC_verb) printf ("Output dataset will have %d sub-bricks\n\n", nbricks);
00553 
00554 
00555   /*----- Open threshold dataset -----*/
00556   thr_dset = THD_open_dataset (SC_thr_filename);
00557   nxyz = DSET_NVOX (thr_dset);
00558 
00559 
00560   /*-- Make an empty copy of threshold dataset, for eventual output --*/
00561   new_dset = EDIT_empty_copy (thr_dset);
00562 
00563 
00564   /*----- Record history of dataset -----*/
00565   tross_Copy_History (thr_dset, new_dset);
00566   if( commandline != NULL ) tross_Append_History( new_dset , commandline ) ;
00567 
00568   
00569   /*----- Delete threshold dataset -----*/
00570   THD_delete_3dim_dataset (thr_dset, False);  thr_dset = NULL ;
00571 
00572 
00573   /*----- Modify some structural properties.  Note that the nbricks
00574           just make empty sub-bricks, without any data attached. -----*/
00575   ierror = EDIT_dset_items (new_dset,
00576                             ADN_prefix,          SC_output_prefix,
00577                             ADN_directory_name,  SC_session,
00578                             ADN_type,            HEAD_FUNC_TYPE,
00579                             ADN_func_type,       FUNC_BUCK_TYPE,
00580                             ADN_ntt,             0,               /* no time */
00581                             ADN_nvals,           nbricks,
00582                             ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,  
00583                             ADN_none ) ;
00584   
00585   if( ierror > 0 )
00586     {
00587       sprintf (message, 
00588               " %d errors in attempting to create bucket dataset! ", 
00589               ierror);
00590       SC_error (message);
00591     }
00592   
00593   if (THD_is_file(DSET_HEADNAME(new_dset))) 
00594     {
00595       sprintf (message,
00596               " Output dataset file %s already exists--cannot continue! ",
00597               DSET_HEADNAME(new_dset));
00598       SC_error (message);
00599     }
00600 
00601 
00602   /*----- Allocate memory -----*/
00603   bar  = (byte **) malloc (sizeof(byte *) * nbricks);
00604   MTEST (bar);
00605   
00606 
00607   /*----- Build lowest level of cluster hierarchy -----*/
00608   for (ivox = 0;  ivox < SC_nvox;  ivox++)
00609     {
00610       /*----- Allocate space for parameter vector -----*/
00611       parameters = (float *) malloc (sizeof(float) * SC_dimension);
00612       MTEST (parameters);
00613 
00614       /*----- Copy the parameter array -----*/
00615       for (ip = 0;  ip < SC_dimension;  ip++)
00616         parameters[ip] = SC_parameters[ip][ivox];
00617       
00618       /*----- If using stat. dist., transform the parameter vector -----*/
00619       if (SC_statdist)
00620         {
00621           array_to_vector (SC_dimension, parameters, &v);
00622           vector_multiply (sinv, v, &av);
00623           vector_to_array (av, parameters);
00624         }
00625 
00626       /*----- Create new cluster containing single voxel -----*/
00627       ixyz = SC_voxels[ivox];
00628       head_clust = new_cluster (ixyz, parameters, head_clust);
00629     }
00630 
00631 
00632   /*----- Deallocate memory for parameter data -----*/
00633   free (SC_voxels);   SC_voxels = NULL;
00634   for (ip = 0;  ip < SC_dimension;  ip++)
00635     {
00636       free (SC_parameters[ip]);   SC_parameters[ip] = NULL;
00637     }
00638   free (SC_parameters);   SC_parameters = NULL;
00639 
00640 
00641   /*----- Agglomerate clusters, one-by-one -----*/
00642   for (iclust = SC_nvox;  iclust > 0;  iclust--)
00643     {
00644       if (SC_verb && (iclust % 100 == 0))
00645         printf ("# Clusters = %d \n", iclust);
00646 
00647       if (iclust <= nbricks)
00648         {
00649           /*----- Sort clusters in order of size -----*/
00650           head_clust = sort_clusters (head_clust);
00651 
00652           /*----- Print cluster centroid parameters -----*/
00653           if (SC_verb)
00654             {
00655               printf ("\n\n# Clusters = %d \n\n", iclust);
00656               print_all_clusters (head_clust, s);
00657             }
00658      
00659           /*----- allocate memory for output sub-brick -----*/
00660           ibrick = iclust-1;
00661           bar[ibrick]  = (byte *) malloc (sizeof(byte) * nxyz);
00662           MTEST (bar[ibrick]);
00663 
00664           /*----- Save clusters into output sub-brick -----*/
00665           for (ixyz = 0;  ixyz < nxyz;  ixyz++)    
00666             bar[ibrick][ixyz] = 0;
00667           save_all_clusters (head_clust, bar[ibrick]); 
00668 
00669           /*----- attach bar[ib] to be sub-brick #ibrick -----*/
00670           EDIT_substitute_brick (new_dset, ibrick, MRI_byte, bar[ibrick]);
00671 
00672         }
00673 
00674       /*----- Agglomerate clusters -----*/
00675       if (iclust > 1)
00676         head_clust = agglomerate_clusters (head_clust, 
00677                                            SC_verb*(iclust <= nbricks));
00678           
00679 
00680     }
00681 
00682 
00683   /*----- Deallocate memory -----*/
00684   vector_destroy (&v);
00685   vector_destroy (&av);
00686   matrix_destroy (&s);
00687   matrix_destroy (&sinv);
00688   destroy_cluster (head_clust);
00689 
00690 
00691   /*----- Return hierarchical clustering -----*/
00692   return (new_dset);
00693   
00694 }
00695 
00696 
00697 /*---------------------------------------------------------------------------*/
00698 
00699 int main( int argc , char * argv[] )
00700 
00701 {
00702   THD_3dim_dataset * clust_dset = NULL;    /* hierarchical clusters data set */
00703   int ip;                                  /* parameter index */
00704 
00705   
00706   /*----- Identify software -----*/
00707   printf ("\n\n");
00708   printf ("Program:          %s \n", PROGRAM_NAME);
00709   printf ("Author:           %s \n", PROGRAM_AUTHOR); 
00710   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00711   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00712   printf ("\n");
00713 
00714 
00715    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00716 
00717    machdep() ; 
00718    { int new_argc ; char ** new_argv ;
00719      addto_args( argc , argv , &new_argc , &new_argv ) ;
00720      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00721    }
00722 
00723 
00724   /*----- Initialize program:  get all operator inputs; get indices
00725     for voxels above threshold; get parameter vectors for all voxels 
00726     above threshold  -----*/
00727   initialize_program (argc, argv);
00728 
00729 
00730   /*----- Perform agglomerative hierarchical clustering -----*/
00731   clust_dset = form_clusters ();
00732 
00733   
00734   /*----- Output the hierarchical clustering dataset -----*/
00735   if( SC_verb ) fprintf(stderr,"Computing sub-brick statistics\n") ;
00736   THD_load_statistics( clust_dset ) ;
00737 
00738   THD_write_3dim_dataset( NULL,NULL , clust_dset , True ) ;
00739   if( SC_verb ) fprintf(stderr,"Wrote output to %s\n", DSET_BRIKNAME(clust_dset) );
00740   
00741 
00742   /*----- Deallocate memory for cluster dataset -----*/   
00743   THD_delete_3dim_dataset( clust_dset , False ) ; clust_dset = NULL ;
00744   
00745   exit(0) ;
00746   
00747 
00748 }
00749 
00750 
00751 
00752 
00753 
00754 
00755