src/avovk/3dkmeans.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "../mrilib.h"
#include "cluster_floatNOMASK.h"
#include "thd_segtools_fNM.h"


static void display_help(int detail)
{ printf ("3d+t Clustering segmentation, command-line version.\n");
  printf ("    Based on The C clustering library.\n");
  printf ("    Copyright (C) 2002 Michiel Jan Laurens de Hoon.\n");
  printf ("USAGE: 3dkmeans [options]\n");
  printf ("options:\n");
  printf ("  -v, --version Version information\n");
  printf ("  -f filename: Input data to be clustered.   \n");
  printf ("                You can specify multiple filenames in sequence\n"
          "                and they will be catenated internally.\n"
          "         e.g: -f F1+orig F2+orig F3+orig ...\n"
          "           or -f F1+orig -f F2+orig -f F3+orig ...\n" );
  printf ("  -input filename: Same as -f\n");
  printf(
 " -mask mset   Means to use the dataset 'mset' as a mask:\n"
 "                 Only voxels with nonzero values in 'mset'\n"
 "                 will be printed from 'dataset'.  Note\n"
 "                 that the mask dataset and the input dataset\n"
 "                 must have the same number of voxels.\n"
 " -mrange a b  Means to further restrict the voxels from\n"
 "                 'mset' so that only those mask values\n"
 "                 between 'a' and 'b' (inclusive) will\n"
 "                 be used.  If this option is not given,\n"
 "                 all nonzero values from 'mset' are used.\n"
 "                 Note that if a voxel is zero in 'mset', then\n"
 "                 it won't be included, even if a < 0 < b.\n"
 " -cmask 'opts' Means to execute the options enclosed in single\n"
 "                  quotes as a 3dcalc-like program, and produce\n"
 "                  produce a mask from the resulting 3D brick.\n"
 "       Examples:\n"
 "        -cmask '-a fred+orig[7] -b zork+orig[3] -expr step(a-b)'\n"
 "                  produces a mask that is nonzero only where\n"
 "                  the 7th sub-brick of fred+orig is larger than\n"
 "                  the 3rd sub-brick of zork+orig.\n"
 "        -cmask '-a fred+orig -expr 1-bool(k-7)'\n"
 "                  produces a mask that is nonzero only in the\n"
 "                  7th slice (k=7); combined with -mask, you\n"
 "                  could use this to extract just selected voxels\n"
 "                  from particular slice(s).\n"
 "       Notes: * You can use both -mask and -cmask in the same\n"
 "                  run - in this case, only voxels present in\n"
 "                  both masks will be dumped.\n"
 "              * Only single sub-brick calculations can be\n"
 "                  used in the 3dcalc-like calculations -\n"
 "                  if you input a multi-brick dataset here,\n"
 "                  without using a sub-brick index, then only\n"
 "                  its 0th sub-brick will be used.\n"
 "              * Do not use quotes inside the 'opts' string!\n"
 "\n");
#if 0
   /* Not implemented */
  printf ("  -cg a|m       Specifies whether to center each row\n"
          "                in the data\n"
          "                a: Subtract the mean of each row\n"
          "                m: Subtract the median of each row\n"
          "                (default is no centering)\n");
  printf ("  -ng           Specifies to normalize each row in the data\n"
          "                (default is no normalization)\n");
  printf ("  -ca a|m       Specifies whether to center each column \n"
          "                in the data\n"
          "                a: Subtract the mean of each column\n"
          "                m: Subtract the median of each column\n"
          "                (default is no centering)\n");
  printf ("  -na           Specifies to normalize each column in the data\n"
          "                (default is no normalization)\n");
#endif
  printf ("  -u jobname    Allows you to specify a different name for the \n"
          "                output files.\n"
          "                (default is derived from the input file name)\n");
  printf ("  -prefix PREFIX Allows you to specify a prefix for the output \n"
          "                 volumes. Default is the same as jobname\n"
          "                 There are two output volumes, one for the cluster\n"
          "                 membership and one with distance measures.\n"
          "                 The distance dataset, mostly for debugging puposes\n"
          "                 is formatted as follows:\n"
          "                 Sub-brick 0: Dc = 100*(1-Ci)+100*Di/(Dmax)\n"
          "                 with Ci the cluster number for voxel i, Di the \n"
          "                 distance of voxel i to the centroid of its \n"
          "                 assigned cluster, Dmax is the maximum distance in\n"
          "                 cluster Ci.\n"
          "                 Sub-bricks 1..k: Dc0k contains the distance of a\n"
          "                 voxel's data to the centroid of cluster k.\n"
          "                 Sub-brick k+1: Dc_norm = (1.0-Di/Ei)*100.0, where \n"
          "                 Ei is the smallest distance of voxel i to \n"
          "                 the remaining clusters that is larger than Di.\n");
  printf ("  -g [0..8]     Specifies distance measure for clustering\n" );
  printf ("                Note: Weight is a vector as long as the signatures\n"
          "                and used when computing distances. However for the\n"
          "                moment, all weights are set to 1\n"
          "                0: No clustering\n"
          "                1: Uncentered correlation distance\n"
          "                    Same as Pearson distance, except\n"
          "                    the means of v and s are not removed\n"
          "                    when computing correlation.\n"
          "                2: Pearson distance\n"
          "                    = (1-Weighted_Pearson_Correlation(v,s))\n"
          "                3: Uncentered correlation distance, absolute value\n"
          "                    Same as abs(Pearson distance), except\n"
          "                    the means of v and s are not removed\n"
          "                    when computing correlation.\n"
          "                4: Pearson distance, absolute value\n"
          "                    = (1-abs(Weighted_Pearson_Correlation(v,s)))\n"
          "                5: Spearman's rank distance\n"
          "                    = (1-Spearman_Rank_Correlation(v,s))\n"
          "                   No weighting is used\n"
          "                6: Kendall's distance\n"
          "                    = (1-Kendall_Tau(v,s))\n"
          "                   No weighting is used\n"
          "                7: Euclidean distance between v and s\n"
          "                    = 1/sum(weight) * sum(weight[i]*(v[i]-s[i])^2)\n"
          "                8: City-block distance\n"
          "                    = 1/sum(weight) * sum(weight[i]*abs(v[i]-s[i]))\n"
          "\n"
          "       (default for -g is 1, 7 if input has one value per voxel)\n"
          "\n");
#if 0
  printf ("  -k number     Specifies whether to run k-means clustering\n"
          "                instead of hierarchical clustering, and the number\n"
          "                of clusters k to use. \n"
          "                Default is kmeans with k = 3 clusters\n");
#endif
  printf ("  -k number     Specify number of clusters\n");
  printf ("  -remap  METH  Reassign clusters numbers based on METH:\n"
          "                   NONE: No remapping (default)\n"
          "                   COUNT: based on cluster size ascending\n"
          "                  iCOUNT: COUNT, descending\n"
          "                   MAG:  based on ascending magnitude of centroid\n"
          "                  iMAG: MAG, descending\n");
  printf ("  -labeltable LTFILE: Attach labeltable LTFILE to clustering\n"
          "                      output. This labeltable will overwrite\n"
          "                      a table that is taken from CLUST_INIT\n"
          "                      should you use -clust_init option.\n");
  printf ("  -clabels LAB1 LAB2 ...: Provide a label for each cluster.\n"
          "                          Labels cannot start with '-'.\n");
  printf ("  -clust_init CLUST_INIT: Specify a dataset to initialize \n"
          "                          clustering. This option sets -r 0 .\n"
          "                          If CLUST_INIT has a labeltable and \n"
          "                          you do not specify one then CLUST_INIT's\n"
          "                          table is used for the output\n");
#if 0
  printf ("  -c number     Force the program to do hierarchical clsutering\n"
          "                and specifies the number of clusters for tree\n"
          "                cutting after hierarchical clustering.\n"
          "       Options -c and -k are mutually exclusive\n");
#endif
  printf ("  -r number     For k-means clustering, the number of times the\n"
          "                k-means clustering algorithm is run\n"
          "                (default: 0 with -clust_init, 1 otherwise)\n");
/*  printf ("  -m [msca]     Specifies which hierarchical clustering method to\n"
          "                use:\n"
          "                m: Pairwise complete-linkage\n"
          "                s: Pairwise single-linkage\n"
          "                c: Pairwise centroid-linkage\n"
          "                a: Pairwise average-linkage\n"
          "                (default: m)\n"); */
  printf ("  -rsigs SIGS   Calculate distances from each voxel's signature\n"
          "                to the signatures in SIGS. \n"
          "                SIGS is a multi-column 1D file with each column\n"
          "                being a signature.\n"
          "                The output is a dset the same size as the input\n"
          "                with as many sub-bricks as there are columns in \n"
          "                SIGS.\n"
          "                With this option, no clustering is done.\n");
  printf ("  -verb         verbose \n");
  printf ("  -write_dists  Output text files containing various measures.\n"
          "                FILE.kgg.1D : Cluster assignments \n"
          "                FILE.dis.1D : Distance between clusters\n"
          "                FILE.cen.1D : Cluster centroids\n"
          "                FILE.info1.1D: Within cluster sum of distances\n"
          "                FILE.info2.1D: Maximum distance within each cluster\n"
          "                FILE.vcd.1D: Distance from voxel to its centroid\n");
  printf ("  -voxdbg I J K Output debugging info for voxel I J K\n");
  printf ("  -seed SEED    Seed for the random number generator.\n"
          "                Default is 1234567\n");
  EXRETURN;
}


/* ========================================================================= */

int main(int argc, char **argv)
{
   int ii=0, ncol=0, nrow=0, nl=0, nc=0;
   int i = 1;
   char* filename[256];
   int l = 0;
   char method = 'm';
   char cg = '\0';
   char ca = '\0';
   int ng = 0;
   int na = 0;
   char *prefix = NULL;
   char *signame=NULL;
   THD_3dim_dataset *in_set=NULL, *clust_set=NULL, *clust_init=NULL;
   THD_3dim_dataset *mask_dset=NULL, *dist_set=NULL;
   byte *cmask=NULL ; int ncmask=0 ;
   byte *mask=NULL;
   int nmask=-1, mnx=-1, mny=-1, mnz=-1, iset=0, N_iset=0, mnxyz=-1;
   float mask_bot=666.0 , mask_top=-666.0 ;
   OPT_KMEANS oc;
   float *dvec=NULL, **D=NULL;
   int n = 0, Ncoltot=0, nc0=0, nx=0, ny=0, nz=0;
   char *prefixvcd = NULL, *clust_init_name=NULL;


   mainENTRY("3dkmeans"); machdep();/* Used to be called 3dAclustering_fNM */
   PRINT_VERSION("3dkmeans"); AUTHOR("avovk") ;

   oc = new_kmeans_oc();

   oc.r = -1;
   oc.k = 0;
   oc.kh = 0;
   oc.jobname = NULL;
   oc.distmetric = '*';
   oc.verb = 0;
   oc.rand_seed = 1234567;
   oc.remap = NONE;
   oc.user_labeltable=NULL;
   oc.nclabels=0;

   for (i=0; i<4; ++i) oc.voxdebug[i] = -1;
   N_iset = 0;
   filename[N_iset] = NULL;


   i = 1;
   while (i < argc)
   { const char* const argument = argv[i];
    i++;
    if (strlen(argument)<2)
    { printf("ERROR: missing argument\n");
      RETURN(1);
    }
    if (argument[0]!='-')
    { printf("ERROR: unknown argument %s\n", argument);
      RETURN(1);
    }
    if(!strcmp(argument,"--version") || !strcmp(argument,"-v"))
    { clusterlib_display_version();
      RETURN(0);
    }
    if(     !strcmp(argument,"--help")
         || !strcmp(argument,"-h")
         || !strcmp(argument,"-help") )
    { display_help(strlen(argument) > 3 ? 2:1);
      RETURN(0);
    }
    if(     !strcmp(argument,"--verb")
         || !strcmp(argument,"-verb") )
    { oc.verb +=1;
      continue;
    }
    if(     !strcmp(argument,"--write_dists")
         || !strcmp(argument,"-write_dists") )
    { oc.writedists=1;
      continue;
    }
    if(     !strcmp(argument,"--remap")
         || !strcmp(argument,"-remap") )
    {
      if (i<argc) {
       if (!strcmp("COUNT",argv[i])) {
         oc.remap = COUNT;
       } else if (!strcmp("iCOUNT",argv[i])) {
         oc.remap = iCOUNT;
       } else if (!strcmp("MAG",argv[i])) {
         oc.remap = MAG;
       } else if (!strcmp("iMAG",argv[i])) {
         oc.remap = iMAG;
       } else if (!strcmp("NONE",argv[i])) {
         oc.remap = NONE;
       } else { printf ("Error reading command line argument for -remap\n");
        RETURN(1);
       }
      } else {
         printf ("Need parameter after -remap\n");
         RETURN(1);
      }
      i++;
      continue;
    }
    if(!strcmp(argument,"-cg"))
    { if (i==argc || strlen(argv[i])>1 || !strchr("am",argv[i][0]))
      { printf ("Error reading command line argument cg\n");
        RETURN(1);
      }
      cg = argv[i][0];
      i++;
      continue;
    }
    if(!strcmp(argument,"-ca"))
    { if (i==argc || strlen(argv[i])>1 || !strchr("am",argv[i][0]))
      { printf ("Error reading command line argument ca\n");
        RETURN(1);
      }
      ca = argv[i][0];
      i++;
      continue;
    }
    if(!strcmp(argument,"-labeltable"))
    { if (i==argc)
      { printf ("Error: Need filename after -labeltable\n");
        RETURN(1);
      }
      oc.user_labeltable = argv[i];
      i++;
      continue;
    }
    if(!strcmp(argument,"-prefix"))
    { if (i==argc)
      { printf ("Error: Need name after -prefix\n");
        RETURN(1);
      }
      prefix = strdup(argv[i]);
      i++;
      continue;
    }
    if(!strcmp(argument,"-clust_init"))
    { if (i==argc)
      { printf ("Error: Need dset after -clust_init\n");
        RETURN(1);
      }
      clust_init_name = argv[i];
      i++;
      continue;
    }
    if(!strcmp(argument,"-voxdbg"))
    { if (i+2==argc)
      { printf ("Error: Need 3 integers after -voxedbg\n");
        RETURN(1);
      }
      oc.voxdebug[0] = atoi(argv[i]);i++;
      oc.voxdebug[1] = atoi(argv[i]);i++;
      oc.voxdebug[2] = atoi(argv[i]);i++;
      continue;
    }

    if(!strcmp(argument,"-rsigs"))
    { if (i==argc)
      { printf ("Error: Need name after -rsigs\n");
        RETURN(1);
      }
      signame = argv[i];
      i++;
      continue;
    }

    if(!strcmp(argument,"-seed"))
    { if (i==argc)
      { printf ("Error: Need a +ve integer after -seed\n");
        RETURN(1);
      }
      oc.rand_seed = atoi(argv[i]);
      if (oc.rand_seed <=0) {
        printf ("Error: seed must be > 0\n");
        RETURN(1);
      }
      i++;
      continue;
    }

    if( strncmp(argument,"-mask",5) == 0 ){
       if( mask_dset != NULL )
         ERROR_exit("Cannot have two -mask options!\n") ;
       if( i >= argc )
         ERROR_exit("-mask option requires a following argument!\n");
       mask_dset = THD_open_dataset( argv[i] ) ;
       if( mask_dset == NULL )
         ERROR_exit("Cannot open mask dataset!\n") ;
       if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex )
         ERROR_exit("Cannot deal with complex-valued mask dataset!\n");
       i++ ; continue ;
    }

    if( strncmp(argument,"-mrange",5) == 0 ){
      if( i+1 >= argc )
        ERROR_exit("-mrange option requires 2 following arguments!\n");
       mask_bot = strtod( argv[  i] , NULL ) ;
       mask_top = strtod( argv[++i] , NULL ) ;
       if( mask_top < mask_top )
         ERROR_exit("-mrange inputs are illegal!\n") ;
       i++ ; continue ;
    }

    if( strcmp(argument,"-cmask") == 0 ){  /* 16 Mar 2000 */
       if( i >= argc )
          ERROR_exit("-cmask option requires a following argument!\n");
       cmask = EDT_calcmask( argv[i] , &ncmask, 0 ) ;
       if( cmask == NULL ) ERROR_exit("Can't compute -cmask!\n");
       i++ ; continue ;
    }

    if(!strcmp(argument,"-ng"))
    { ng = 1;
      continue;
    }
    if(!strcmp(argument,"-na"))
    { na = 1;
      continue;
    }

    if(!strcmp(argument,"-f") || !strcmp(argument,"-input")) {
         if (i==argc)
        { printf ("Error reading command line argument -f (or -input): "
                  "no file name specified\n");
          RETURN(1);
        }
        do {
         filename[N_iset] = argv[i];
         if (N_iset > 100) {
            printf ("Error: Too many input files!\n");
            RETURN(1);
         }
         ++N_iset; filename[N_iset] = NULL;
         i++;
        } while (i< argc && argv[i][0] != '-');
        continue;
    }
    if(!strcmp(argument,"-clabels")) {
         if (i==argc)
        { printf ("Error reading command line argument -clabels: "
                  "no labels specified\n");
          RETURN(1);
        }
        do {
         oc.clabels[oc.nclabels] = argv[i];
         if (oc.nclabels > 400) {
            printf ("Error: Too many labels!\n");
            RETURN(1);
         }
         ++oc.nclabels;
         i++;
        } while (i< argc && argv[i][0] != '-');
        continue;
    }

    switch (argument[1])
    { case 'l': l=1; break;
      case 'u':
      { if (i==argc)
        { printf ("Error reading command line argument u: "
                  "no job name specified\n");
          RETURN(1);
        }
        oc.jobname = clusterlib_setjobname(argv[i],0);
        i++;
        break;
      }
      #if 0
      case 'f':
      { if (i==argc)
        { printf ("Error reading command line argument f: "
                  "no file name specified\n");
          RETURN(1);
        }
        do {
         filename[N_iset] = argv[i];
         if (N_iset > 100) {
            printf ("Error: Too many input files!\n");
            RETURN(1);
         }
         ++N_iset; filename[N_iset] = NULL;
         i++;
        } while (i< argc && argv[i][0] != '-');
        break;
      }
      #endif
      case 'g':
      { int g;
        if (i==argc)
        { printf ("Error reading command line argument g: parameter missing\n");
          RETURN(1);
        }
        g = clusterlib_readnumber(argv[i]);
        if (g < 0 || g > 9)
        { printf ("Error reading command line argument g: "
                  "should be between 0 and 9 inclusive\n");
          RETURN(1);
        }
        i++;
        oc.distmetric = clusterlib_getmetric(g);
        break;
      }

      case 'k':
      { if (i==argc)
        { printf ("Error reading command line argument k: "
                  "parameter missing\n");
          RETURN(1);
        }
        if (oc.kh > 0) {
            ERROR_message("-k and -c options are mutually exclusive\n");
            RETURN(1);
        }
        oc.k = clusterlib_readnumber(argv[i]);
        if (oc.k < 1)
        { printf ("Error reading command line argument k: "
                  "a positive integer is required\n");
          RETURN(1);
        }
        i++;
        break;
      }
    case 'c':
      { if (i==argc)
        { printf ("Error reading command line argument c: parameter missing\n");
          RETURN(1);
        }
        if (oc.k > 0) {
            ERROR_message("-k and -c options are mutually exclusive\n");
            RETURN(1);
        }
        oc.kh = clusterlib_readnumber(argv[i]);
        if (oc.kh < 1)
        { printf ("Error reading command line argument c: "
                  "a positive integer is required\n");
          RETURN(1);
        }
        i++;
        break;
      }
   case 'r':
      { if (i==argc)
        { printf ("Error reading command line argument r: parameter missing\n");
          RETURN(1);
        }
        oc.r = clusterlib_readnumber(argv[i]);
        if (oc.r < 0)
        { printf ("Error reading command line argument r: "
                  "a >= 0 integer is required\n");
          RETURN(1);
        }
        i++;
        break;
      }
   case 'm':
      { if (i==argc || strlen(argv[i])>1 || !strchr("msca",argv[i][0]))
	  { printf ("Error reading command line argument m: "
               "should be 'm', 's', 'c', or 'a'\n");
	    RETURN(1);
	  }
        method = argv[i][0];
        i++;
        break;
      }
      default:
         printf ("Unknown option %s\n", argv[i-1]);
         suggest_best_prog_option(argv[0], argv[i-1]);
         RETURN(1);
    }

   }
   if (argc < 2) {
      printf ("Too few options!");
      display_help(0);
      RETURN(1);
   }

   if (oc.k <= 0 && oc.kh <= 0) oc.k = 3;

   if (clust_init_name) {
      if (oc.r > 0) {
         ERROR_message("Can't use -clust_init and -r > 0");
         RETURN(1);
      } else {
         oc.r = 0;
      }
   }

   if (oc.r == -1) oc.r = 1;


   if(oc.jobname == NULL) oc.jobname = clusterlib_setjobname(filename[0],1);

   if (oc.nclabels && oc.nclabels != oc.k && oc.nclabels != oc.kh) {
      ERROR_message("Have %d labels, but %d clusters\n",
               oc.nclabels, oc.k > 0 ? oc.k : oc.kh);
      RETURN(1);
   }

   /* load dsets and prepare array data for sending to clustering functions */

   if (!prefix) {
      prefix = (char *)calloc(strlen(oc.jobname)+15,sizeof(char));
      sprintf(prefix,"%s.k%d",oc.jobname, oc.k);
      THD_force_ok_overwrite(1) ;   /* don't worry about overwriting */
   }

   /* ------------- Mask business -----------------*/
   if( mask_dset == NULL ){
      mask = NULL ;
      if( oc.verb )
         INFO_message("Using all voxels in the entire dataset (no mask)\n") ;
   } else {
      mnx = DSET_NX(mask_dset);
      mny = DSET_NY(mask_dset);
      mnz = DSET_NZ(mask_dset);
      mnxyz = mnx*mny*mnz;
      mask = THD_makemask( mask_dset , 0 , mask_bot, mask_top ) ;
      if( mask == NULL ) ERROR_exit("Can't make mask") ;
      nmask = THD_countmask( mnx*mny*mnz  , mask ) ;
      if( oc.verb )
         INFO_message("%d voxels in the [%dx%dx%d] mask",nmask, mnx, mny, mnz) ;
      if( nmask <= 0 ) ERROR_exit("No voxels in the mask!\n") ;
      DSET_delete(mask_dset) ;
   }

   if( cmask != NULL ){
      if( mask != NULL ){
         if (mnxyz != ncmask) ERROR_exit("Mask and cmask dimension mismatch") ;
         for( ii=0 ; ii < mnxyz ; ii++ )
            mask[ii] = (mask[ii] && cmask[ii]) ;
         free(cmask) ;
         nmask = THD_countmask( mnxyz , mask ) ;
         if( nmask <= 0 ) ERROR_exit("No voxels in the mask+cmask!\n") ;
         if( oc.verb ) INFO_message("%d voxels in the mask+cmask\n",nmask) ;
      } else {
         mnx = -1; mny = -1; mnz = -1; /* unknown */
         mnxyz = ncmask;
         mask = cmask ;
         nmask = THD_countmask( mnxyz , mask ) ;
         if( nmask <= 0 ) ERROR_exit("No voxels in the cmask!\n") ;
         if( oc.verb ) INFO_message("%d voxels in the cmask\n",nmask) ;
      }
   }

   if (signame) {
      MRI_IMAGE *im = NULL;
      float *far = NULL;

      /* catenate all input dsets */
      if (N_iset == 1) {
         in_set = THD_open_dataset(filename[0]);
         CHECK_OPEN_ERROR(in_set,filename[0]) ;
         if (oc.voxdebug[0] >= 0) {
            /* setup for debugging */
            oc.voxdebug[3] = oc.voxdebug[0] + oc.voxdebug[1]*DSET_NX(in_set) +
                              oc.voxdebug[2]*DSET_NX(in_set)*DSET_NY(in_set);
         } else oc.voxdebug[3] = -1;
      } else {
         /* you'll need to read and catenate on the fly ... */
         ERROR_exit( "Not ready to deal with more than one input.\n"
                     "Consdier catenating the input externally.\n"
                     "Let me know if it becomes annoying ...\n");
      }

      /* load the set of distance files */
      im = mri_read_1D (signame);
      far = MRI_FLOAT_PTR(im);
      /* Now call distance function */
      if (!thd_Adist ( in_set,
                       mask,
                       far, im->ny,
                       &dist_set ,
                       oc)) {
         ERROR_exit("Failed in thd_Acluster");
      }
      if (im) mri_free(im); im = NULL; far = NULL;
      /* add history to output data and write them  out */
      if( oc.verb && dist_set)
            ININFO_message("\nWriting datasets: %s",prefix) ;
      if (dist_set) {
         EDIT_dset_items(  dist_set , ADN_prefix  , prefix, ADN_none);
         tross_Copy_History( in_set , dist_set ) ;
         tross_Make_History( "3dkmeans" , argc, argv , dist_set ) ;
         DSET_write(dist_set); DSET_unload(dist_set);
         DSET_delete(dist_set); dist_set = NULL;
      }
   } else {
      /* Doing clustering function */
      Ncoltot=0;
      /* Read in dset(s) and create D */
      for (iset = 0; iset < N_iset; ++iset) {
         if (oc.verb) fprintf(stderr,"Reading %s's header (%d/%d), ",
                                    filename[iset], iset+1, N_iset);
         in_set = THD_open_dataset(filename[iset]);
         CHECK_OPEN_ERROR(in_set,filename[iset]) ;
         if (oc.voxdebug[0] >= 0) {
            /* setup for debugging */
            oc.voxdebug[3] = oc.voxdebug[0] + oc.voxdebug[1]*DSET_NX(in_set) +
                              oc.voxdebug[2]*DSET_NX(in_set)*DSET_NY(in_set);
         } else oc.voxdebug[3] = -1;
         if (iset == 0) {
            ncol = DSET_NVALS(in_set);
            nrow = DSET_NVOX(in_set);
            nx = DSET_NX(in_set); ny = DSET_NY(in_set); nz = DSET_NZ(in_set);
            if (  mask                                              &&
                  ( ( (mnx >= 0 && mnx != DSET_NX(in_set)) ||
                      (mny >= 0 && mny != DSET_NY(in_set)) ||
                      (mnz >= 0 && mnz != DSET_NZ(in_set))  )    ||
                    ( mnxyz != nx*ny*nz ) ) ) {
               ERROR_exit("Dimension mismatch between mask (%d=%dx%dx%d)\n"
                          "                 and input dset (%d=%dx%dx%d)",
                          mnxyz, mnx, mny,mnz, nx*ny*nz, nx,ny,nz);
            }
            if (!mask) nmask = DSET_NVOX(in_set);
         } else { /* check for consistency with previous input */
            if (  (nx != DSET_NX(in_set) ||
                   ny != DSET_NY(in_set) ||
                   nz != DSET_NZ(in_set) ) ) {
               ERROR_exit( "Dimension mismatch between input dset"
                           " %s and preceding ones", filename[iset]);
            }
            ncol = DSET_NVALS(in_set);
         }
         if (oc.verb) fprintf(stderr," %d cols\n ",
                                    ncol);
         Ncoltot += ncol;
         /* get rid of dset */
         DSET_delete (in_set);
      }

      if (oc.distmetric == '*') {
         if (Ncoltot == 1) oc.distmetric = 'e';
         else oc.distmetric = 'u';
      }

      /* Now allocate for D */
      D = (float **)calloc(sizeof(float*), nmask);
      for (ii=0;ii<(nmask);++ii) {
         if (!(D[ii] = (float *)calloc(sizeof(float), Ncoltot))) {
            fprintf(stderr,"ERROR: Failed while allocating %dx%d float matrix\n",
                           nmask, Ncoltot);
            RETURN(1);
         }
      }

      dvec = (float * )malloc(sizeof(float)*Ncoltot) ;  /* array to hold series
                                                         longer than needed, but
                                                         less hassle*/
      nc0 = 0;
      for (iset = 0; iset < N_iset; ++iset) {
         if (oc.verb)
            fprintf(stderr,"Patience, rereading %s...\n", filename[iset]);
         in_set = THD_open_dataset(filename[iset]);
         DSET_load(in_set) ; ncol = DSET_NVALS(in_set);

         if (oc.verb) {
            ININFO_message("Filling cols [%d..%d] of D(%dx%d) (mask=%p).\n",
                              nc0,nc0+ncol-1, nmask, Ncoltot, mask);
         }
         ii = 0;
         for (nl=0; nl<DSET_NVOX(in_set); ++nl) {
            if (!mask || mask[nl]) {
               THD_extract_array( nl , in_set , 0 , dvec ) ;
               for (nc=0; nc<ncol; ++nc) D[ii][nc0+nc] = dvec[nc];
               ++ii;
            }
         }
         nc0 += ncol;
         if (iset != N_iset-1) DSET_delete(in_set);
         else DSET_unload(in_set);
      }
      free(dvec); dvec = NULL;

      /* Load initialization */
      if (clust_init_name) {
         if (oc.verb) fprintf(stderr,"Reading %s's header, ",
                                    clust_init_name);
         if (!(clust_init = THD_open_dataset(clust_init_name))) {
            ERROR_exit("Failed to read initialization dset %s\n",
                        clust_init_name);
         }
         DSET_load(clust_init) ;
         CHECK_OPEN_ERROR(clust_init,clust_init_name) ;
      }

      /* Now call clustering function */
      if (!thd_Acluster ( in_set,
                        mask, nmask,
                        &clust_set,
                        &dist_set ,
                        clust_init,
                        oc, D, Ncoltot)) {
         ERROR_exit("Failed in thd_Acluster");
      }

      /* freedom */
      if (D) {
         for (ii=0; ii<nmask; ++ii) if (D[ii]) free(D[ii]);
         free(D); D = NULL;
      }
      /* avovk; make prefix for other datasets, based on input prefix */

      n = 1 + strlen(prefix) + strlen("_vcd");
      prefixvcd = (char *)malloc(n*sizeof(char));
      sprintf (prefixvcd, "%s_vcd", prefix);


      /* add history to output data and write them  out */
      if( oc.verb &&
          (clust_set || dist_set))
        ININFO_message("\nWriting dataset: %s %s",prefix, prefixvcd) ;
      if (clust_set) {
         EDIT_dset_items(  clust_set , ADN_prefix  , prefix, ADN_none);
         tross_Copy_History( in_set , clust_set ) ;
         tross_Make_History( "3dkmeans" , argc, argv , clust_set ) ;
         DSET_write(clust_set); DSET_unload(clust_set);
         DSET_delete(clust_set); clust_set = NULL;
      }
      ININFO_message("\nWriting dataset: %s", prefixvcd) ;
      if (dist_set) {
         EDIT_dset_items(  dist_set , ADN_prefix  , prefixvcd, ADN_none);
         tross_Copy_History( in_set , dist_set ) ;
         tross_Make_History( "3dkmeans" , argc, argv , dist_set ) ;
         DSET_write(dist_set);
         DSET_unload(dist_set);
         DSET_delete(dist_set);
         dist_set = NULL;
      }

   }

   if (prefix) free(prefix); prefix=NULL;
   if (mask) free(mask); mask = NULL;
   if (oc.jobname) free(oc.jobname); oc.jobname = NULL;

   fprintf (stderr,"\n");
   RETURN(0);
}