 /*
  * Khoros: $Id: lviso2.c,v 1.2 1991/10/02 00:23:03 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lviso2.c,v 1.2 1991/10/02 00:23:03 khoros Exp $";
#endif

 /*
  * $Log: lviso2.c,v $
 * Revision 1.2  1991/10/02  00:23:03  khoros
 * HellPatch2
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lviso2.c
 >>>>
 >>>>      Program Name: viso2
 >>>>
 >>>> Date Last Updated: Tue Mar  5 22:20:34 1991 
 >>>>
 >>>>          Routines: lviso2 - the library call for viso2
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#include "vipl/lviso2.h"
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lviso2 - library call for viso2
*
* Purpose:
*    
*    Performs the ISO2 clustering algorithm on an image.
*    
*    
* Input:
*    
*    *img           input image
*    
*    *cc_img        input cluster center image
*    
*    min_pts_allowedmin number vectors/cluster
*    
*    max_cl_allowed max number cluster
*    
*    n_clusters     inital number cluster
*    
*    max_n_iters_isomax iso data iterations
*    
*    max_n_iters_kmeans
*                   max kmeans data iterations
*    
*    border         image border
*    
*    cc_flg         initial cluster locations
*    
*    *init_cluster[]array of initial cluster centers x-y pairs
*    
*    split_factor   splitting factor
*    
*    merge_factor   merging factor
*    
*    placement      splitting placement factor
*    
*    split_converge splitting convergence
*    
*    merge_converge merging convergence
*    
*    *printdev2     output ascii file descriptor
*    
*    
* Output:
*    
*    *img           output image
*    
*    **outcc_img    cluster center image
*    
*    **outvar_img   variance image
*    
*    
*
* Written By: Tom Sauer, Donna Koechner
*    
*    
****************************************************************/


/* -library_def */
lviso2 (img, outcc_img, outvar_img, min_pts_allowed, max_cl_allowed, n_clusters, split_factor, merge_factor, placement, split_converge, merge_converge, max_n_iters_iso, max_n_iters_kmeans, border, init_cluster, cc_flg, cc_img, printdev2  )

struct xvimage *img,            /* input image                 */
               **outcc_img,     /* cluster center image        */
               **outvar_img,    /* variance image              */
               *cc_img;         /* input cluster center image  */

int     min_pts_allowed,        /* min number vectors/cluster  */
        max_cl_allowed,         /* max number cluster          */
        n_clusters,             /* inital number cluster       */
        max_n_iters_iso,        /* max iso data iterations     */
        max_n_iters_kmeans,     /* max kmeans data iterations  */
        border,                 /* image border                */
        cc_flg,                 /* initial cluster locations   */
       *init_cluster[];       /* array of initial cluster centers x-y pairs */

float   split_factor,           /* splitting factor            */
        merge_factor,           /* merging factor              */
        placement,              /* splitting placement factor  */
        split_converge,         /* splitting convergence       */
        merge_converge;         /* merging convergence         */

FILE * printdev2;               /* output ascii file descriptor */
/* -library_def_end */

/* -library_code */
{
    int     nc,                 /* column size of image */
            nr,                 /* row size of image */
            cc_nc,              /* column size of image */
            cc_nr,              /* row size of image */
            vect,               /* number of vectors in image */
            dimension,          /* size of vector - num_bands */
            vector,             /* index to input_vectors */
           *classptr,           /* temp pointer to the class data structure */
            num_vects,          /* dummy variable */
            dim, 
            i,                  /* loop indices */
            j;

    float  *ptr;                /* pointer to the image data char */
    char **load_vector();      /* converts images to an array of vectors */

    int    *output_imagedata;   /* pointer to resulting image char */

    char   *program = "lviso2";

    struct c_struct c_data;     /* clustering information */
    struct iso_params i_parm;   /* iso2 paramenters */



    nr = img -> col_size;       /* number of rows in image */
    nc = img -> row_size;       /* number of columns in image */

    if (cc_flg) {               /* get row and col size of cc_img */
        cc_nr = cc_img -> col_size;
        cc_nc = cc_img -> row_size;
    }

     /* calculate size of image minus border ............................. */
    vect = (nr - (border * 2)) * (nc - (border * 2));
    dimension = img -> num_data_bands;

    /* initialize the iso data structures ................................. */

    c_data.dimension = dimension;
    c_data.n_clusters = n_clusters;
    c_data.max_n_iters_kmeans = max_n_iters_kmeans;
    c_data.n_vects = vect;

    i_parm.max_cl_allowed = max_cl_allowed;
    i_parm.max_n_iters_iso = max_n_iters_iso;
    i_parm.min_pts_allowed = min_pts_allowed;
    i_parm.split_factor = split_factor;
    i_parm.merge_factor = merge_factor;
    i_parm.placement = placement;
    i_parm.split_converge = split_converge;
    i_parm.merge_converge = merge_converge;

  
         /* allocate space for all the data structures used .................*/
    if ( ! get_space(&c_data, &i_parm))
    {
       (void)fprintf(stderr,"%s: Error - get_space() failed", program);  
       return(0);
    }

         /* zero out the class assignment array */
    bzero (c_data.class, (unsigned int) vect * sizeof (int));

    /* Assign image data address to ptr ................................. */

    c_data.input_vectors = (float **) load_vector(img,border,&num_vects,&dim);
    if (c_data.input_vectors == NULL)
    {
       fprintf(stderr,"%s: Error while loading image into vector format\n"
                      ,program);
       return(0);
    }

         /* free the unneeded memory ...................................*/

    (void) free(img -> imagedata);
    (void) free (img -> maps);

    /* Assign cluster center image data to c_center if used as input .......*/

    if (cc_flg) {
        ptr = (float *) cc_img -> imagedata;

        for (i = 0; i < cc_nr * cc_nc; i++) {
            for (j = 0; j < dimension; j++) {
                c_data.c_center[i][j] = ptr[i + (cc_nc * cc_nr * j)];
            }
        }
        (void) free( cc_img -> imagedata);
    }
    else {              /* otherwise use input vectors */
        for (i = 0; i < n_clusters; i++) {
            vector = (init_cluster[i][1] * nc + init_cluster[i][0]);
            for (j = 0; j < dimension; j++) {
                c_data.c_center[i][j] = c_data.input_vectors[vector][j];
            }
        }
    }

        /* ...............  call iso2 routine .........................*/

    if (! iso2(&c_data, &i_parm, printdev2, outcc_img, outvar_img)) {
                
        (void) fprintf (stderr, "lviso2: iso2 Algorithm failed\n");
        return (0);
    }

    /* Allocate space for output image data ........................... */

    output_imagedata = (int *) malloc ((unsigned int) nr * nc * sizeof (int));
    if (output_imagedata == NULL) {
        (void) fprintf (stderr, "%s: insufficient memory available\n", program);
        return (0);
    }

    bzero (output_imagedata, (unsigned int) nr * nc * sizeof (int));

      /* Need to transfer the class image to the output image. We want
       * to preserve the original size of the image, so must worry
       * about the border. Also since the 0 class is defined to be
       * the Null class, must add 1 to the current class assignment.
       * Note: a discarded vector in the iso2 routine is denoted as a -1
       * in the class assignment. Therefore, adding 1 make the class 0
       * which is the Null class, which is what we want
       */

    classptr = c_data.class;
    for (i = border; i < nr - border; i++) 
    {
        for (j = border; j < nc - border; j++) 
        {
           output_imagedata[i * nc + j] = (*classptr) + 1;
           classptr++;
        }
    }
    

       /* adjust header in the output image .........................*/
    img -> imagedata = (char *) output_imagedata;
    img -> data_storage_type = VFF_TYP_4_BYTE;
    img -> num_data_bands = 1;
    img -> map_scheme = VFF_MS_NONE;
    img -> map_row_size = 0;
    img -> map_col_size = 0;
    img -> map_subrow_size = 0;
    img -> color_space_model = VFF_CM_NONE;

    return (1);
}


/**************************************************************
*
* MODULE NAME: iso2
*
*     PURPOSE: 
*               Iso2 is a clustering algorithm which is derived from the 
*               ISODATA clustering algorithm.  Iso2 initially clusters the 
*               data using a standard K-means clustering algorithm. 
*               Then, any cluster, along with its associated vectors, is 
*               discarded if the membership of that cluster is below a user 
*               specified threshold. Next, clusters exhibiting large
*               variances are split in two, and clusters that are too 
*               close together are merged.
*               After these operations, K-means is reiterated, and the iso2
*               rules tested again.  This sequence is repeated until no 
*               clusters are discarded, split, or merged, or until the 
*               maximum number of iterations has been met.
*               
*               The merging and splitting rules used in iso2 are similar 
*               to those used in ISODATA, except that several modifications 
*               have been made which simplify 
*               parameter selection for the user.
*               The most significant change is that 
*               the splitting and merging rules are now relative to the
*               statistics of the data rather than being dependent 
*               on absolute limits, and these relative limits can be made to
*               converge so that splitting and merging are reduced as
*               the algorithm progresses.
*               Secondly, the splitting rules are no longer dependent on a user 
*               specified "desired number of clusters" nor on the iteration index.
*               Defining the splitting and merging requirements so that they
*               are relative to the
*               data set removes the need for the user to know the range
*               of the data, and makes characterization of the algorithm 
*               for the general case possible.  The methods of splitting and
*               merging are described below.
*               
*               
*               SPLITTING:  RELATIVE STANDARD DEVIATION APPROACH
*               
*               The standard deviation in each dimension of 
*               each cluster is computed, and the average of these 
*               standard deviations calculated.  The maximum standard
*               deviation for each cluster is then compared against 
*               the average standard deviation times a user specified
*               scaling factor (splitting factor).  
*               If the maximum standard deviation is greater than 
*               that value, the cluster is split in two in the dimension 
*               of maximum standard deviation.  Two other splitting 
*               conditions also exist: (1) the number of points in the
*               cluster must satisfy a minimum condition, and (2) the 
*               average distance between members of a cluster and their 
*               cluster center must be less than the average of this distance 
*               over all clusters.
*               After all eligible clusters have been split, the splitting
*               factor is updated.
*               A cluster can only be split once per iteration.
*               
*               
*               MERGING:  RELATIVE DISTANCE APPROACH
*               
*               The merging process is similar to the splitting process
*               described above.  The pairwise distances between cluster
*               centers are computed as is done in the original ISODATA 
*               algorithm and, in addition, the average pairwise distance
*               is computed.  If any distance between two cluster centers is 
*               less than the average pairwise distance multiplied by a scaling 
*               factor (merging factor), the two clusters are merged.
*               The user has the option to specify the maximum number of 
*               cluster pairs that can be merged.  Cluster pairs with the
*               smallest pairwise distances are merged first, and a cluster
*               can only be merged once per iteration.  After all eligible
*               clusters are merged, the merging factor is updated.
*               
*
*       INPUT:  c_data - clustering information structure
*               i_parm - iso information structure
*               printdev - file descriptor for output.
*               outcc_img - cluster center output image
*               outvar_img - cluster variance output image
*
*      OUTPUT:  updated c_data, i_parm, outcc_img, outvar_img
*
* CALLED FROM: lviso2()
*
* ROUTINES CALLED: discard(), kmean(), merge2(), split2()
*
**************************************************************/

int iso2 (c_data, i_parm, printdev, outcc_img,
          outvar_img)

   struct c_struct    *c_data;       /* clustering informatiom structure */
   struct iso_params  *i_parm;       /* iso information structure */
   struct xvimage     **outcc_img;   /* cluster center output image */
   struct xvimage     **outvar_img;  /* cluster variance output image */
   FILE               *printdev;     /* file descriptor for output */
{


    int i,j,k;           /* loop counters */
    int count;           
    float trace;         /* tarce of the covariance matrix */
    float *ccptr;
    float *varptr;
    char *comment = NULL;
    struct xvimage     *createimage();


    i_parm->iter_iso = 0;
    i_parm->discarded_vects = 0;

      /**********************************************************
      *  Start the Iso2 algorithm 
      ***********************************************************/

    do
    {
        if (! kmean(c_data, i_parm))
        {
          fprintf(stderr, "iso: kmean failed\n");
          return(0);
        }

        i_parm->change = 0;

/*  ISODATA RULE 1: discard clusters with too few elements ...........*/
/*   may want to keep track of discarded vectors */

        discard(c_data,i_parm);

/*  ISODATA RULE 2: split clusters with large standard deviation .....*/

        split2(c_data,i_parm);

/*  If no clusters were split up, then apply ISODATA RULE 3: merge    */
/*  cluster centers that are close.  Otherwise, the process is        */
/*  repeated starting with K-means. ..................................*/

        if (i_parm->change == FALSE)
        {
            merge2(c_data,i_parm);
        }

/*  Test convergence .................................................*/

        i_parm->iter_iso++;


    } while (i_parm->change && (i_parm->iter_iso < i_parm->max_n_iters_iso) );


    /* Print statistics for clustering routine....................... */

    (void) fprintf (printdev, "viso2 Statistics\n");
    (void) fprintf (printdev, "================\n\n");
    (void) fprintf (printdev, 
                  "\tTotal Number of ISODATA2 Iterations : %d\n\n",i_parm->iter_iso);
    (void) fprintf(printdev, 
                    "\tTotal Number of Clusters = %d\n\n",c_data->n_clusters);

    (void) fprintf (printdev, "\tTotal Number of Discarded Pixels: %d\n\n"
           ,i_parm->discarded_vects);


      /* print out the number of kmeans iters  per iso iteration .......... */

    (void) fprintf (printdev, "\n\nNumber of kmeans Iterations Per Iso Iteration\n");
    (void) fprintf (printdev, "=============================================\n\n");
    (void) fprintf (printdev, "Format: Iso iteration / number kmeans iterations\n");
    for ( i = 0; i < i_parm->iter_iso; i++)
    {
       if ( (i % 7) == 0) fprintf(printdev,"\n\t| ");
       (void) fprintf(printdev, "%d/%d  | ", i+1, c_data->kmean_iter[i]);
    }

      /* print out the total number of points per cluster................. */

    (void) fprintf (printdev, "\n\nNumber of Vectors Per Cluster\n");
    (void) fprintf (printdev, "=============================\n\n");
    (void) fprintf(printdev,
                "Cluster #0 = %d  (Null Cluster -- not including the border)\n",
                 i_parm->discarded_vects);
    for  (i = 0; i< c_data->n_clusters; i++)
         (void) fprintf(printdev,"Cluster #%d = %d\n",i+1,(int) c_data->n_pts[i]);

    
     /* print out the cluster center values............................... */

    (void) fprintf (printdev, 
          "\n\n\nCluster Center Values:        (Vector Size: %d)\n", c_data->dimension);
    (void) fprintf (printdev, "=====================");

    for (i = 0; i < c_data->n_clusters; i++) {
        count = 1;
        (void) fprintf (printdev, "\n\nCluster #%d\n\n", i);
        for (j = 0; j < c_data->dimension; j++) {
            (void) fprintf (printdev, "%f  ", c_data->c_center[i][j]);
            if (count++ % 5 == 0)
                (void) fprintf (printdev, "\n");
        }
    }


     /* print out the varaince for each  cluster center values.............*/

    (void) fprintf (printdev, 
          "\n\n\nCluster Center Variance Values:      (Vector Size: %d)\n", 
                   c_data->dimension);
    (void) fprintf (printdev, "==============================");
    for (i = 0; i < c_data->n_clusters; i++) {
        count = 1;
        (void) fprintf (printdev, "\n\nCluster #%d\n\n", i);

        for (j = 0; j < c_data->dimension; j++) {
            (void) fprintf (printdev, "%f  ", c_data->var[i][j]);
            if (count++ % 5 == 0)
                (void) fprintf (printdev, "\n");
        }
    }

     /* print out the trace of the Covariance Matrix........................*/
        
    (void) fprintf (printdev, "\n\n\nTrace of Covariance Matrix \n");
    (void) fprintf (printdev, "==========================");

    for (i = 0; i < c_data->n_clusters; i++) {
        (void) fprintf (printdev, "\n\nCluster #%d\t", i);
        trace = 0.0;
        for (j = 0; j < c_data->dimension; j++) {
            trace += (float) c_data->var[i][j];
        }
            (void) fprintf (printdev, "%f\n", trace);
    }



/* ------------  Include this code to print which image vectors belong 
                 to which cluster center class.  --------------------

    (void) fprintf (printdev, 
                        "\n\n\nVectors Cooresponding to Each Cluster\n");
    (void) fprintf (printdev, "=====================================");

    for (i = 0; i < c_data->n_clusters; i++)
    {
        count = 1;
        (void) fprintf (printdev, "\n\ncluster #%d\n", i+1);
        for (j = 0; j < c_data->n_vects; j++)
        {
            if (c_data->class[j] == i)
            {
                (void) fprintf (printdev, "%d  ", j);
                if (count++ % 10 == 0)
                    (void) fprintf (printdev, "\n");
            }
        }
     }

-------------------------------------------------------------------*/


    /* call createimage to create the cluster center value and variance
     * image. size = c_data->n_clusters
     */

   *outcc_img = createimage((unsigned long) c_data->n_clusters, 
                           (unsigned long) 1,
                           (unsigned long) VFF_TYP_FLOAT, 
                           (unsigned long) 1,
                           (unsigned long) c_data->dimension, 
                           comment, 
                           (unsigned long) 0, 
                           (unsigned long) 0, 
                           (unsigned long) VFF_MS_NONE,
                           (unsigned long) VFF_MAPTYP_NONE, 
                           (unsigned long) VFF_LOC_IMPLICIT, 
                           (unsigned long) 0);

   if ( outcc_img == NULL)
   {
       (void) fprintf(stderr,"iso2: cannot allocate space for cc image\n");
       return(0);
   }

   *outvar_img = createimage((unsigned long) c_data->n_clusters, 
                           (unsigned long) 1,
                           (unsigned long) VFF_TYP_FLOAT, 
                           (unsigned long) 1,
                           (unsigned long) c_data->dimension, 
                           comment, 
                           (unsigned long) 0, 
                           (unsigned long) 0, 
                           (unsigned long) VFF_MS_NONE,
                           (unsigned long) VFF_MAPTYP_NONE, 
                           (unsigned long) VFF_LOC_IMPLICIT, 
                           (unsigned long) 0);

   if ( outvar_img == NULL)
   {
       (void) fprintf(stderr,"iso2: cannot allocate space for var image\n");
       return(0);
   }
                           
    ccptr = (float *) outcc_img[0]->imagedata;
    varptr = (float *) outvar_img[0]->imagedata;
    for (i = 0; i < c_data->n_clusters; i++) {
        for (j = 0; j < c_data->dimension; j++) {
            k = i + (j * outcc_img[0]->row_size * outcc_img[0]->col_size);
            ccptr[k] = c_data->c_center[i][j];
            varptr[k] = (float) c_data->var[i][j];
        }
    }

    return(1);
}

/**************************************************************
*
* MODULE NAME: get_space()
*
*     PURPOSE: malloc all the space needed for data structures for the
*              lviso program.
*
*       INPUT: 
*                struct c_struct c_data -- clustering information
*                struct iso_params i_parm  -- iso2 paramenters
*
*      OUTPUT:
*             struct c_struct c_data - with malloced arrays
*             struct iso_params i_parm - with malloced arrays
*
* CALLED FROM:  lviso2()
*
* ROUTINES CALLED: None
*
**************************************************************/

int
get_space(c_data, i_parm)

    struct c_struct *c_data;     /* clustering information */
    struct iso_params *i_parm;   /* iso2 paramenters */

{

    int i;
    char *program = "get_space";

/*  allocate memory for number of points in each cluster array .......*/

    c_data->n_pts = (float *) malloc((unsigned int)sizeof(float) * i_parm->max_cl_allowed);
    if (c_data->n_pts == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing number points per cluster.\n");
        return(0);
    }

/*  allocate memory for distance of cluster pts from their cluster center */

    i_parm->d_avg = (float *) malloc((unsigned int)sizeof(float) * i_parm->max_cl_allowed);
    if (i_parm->d_avg == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing average distance of cluster pointsts.\n");
        return(0);
    }

/*  allocate memory for std dev array of (dim = c_data->dimension) */

    i_parm->sd = (float *) malloc((unsigned int)sizeof(float) * c_data->dimension);
    if (i_parm->sd == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing standard deviation information\n");
        return(0);
    }

/*  allocate memory for merge stat array (dim = n_clusters) */

    i_parm->merge_stat = (float *) malloc((unsigned int)sizeof(float) * i_parm->max_cl_allowed);
    if (i_parm->merge_stat == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing merge statistics information\n");
        return(0);
    }

/*  allocate memory for max std dev array (dim = n_clusters) */

    i_parm->max_sd = (float *) malloc((unsigned int)sizeof(float) * i_parm->max_cl_allowed);
    if (i_parm->max_sd == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing max standard deviation information\n");
        return(0);
    }

/*  allocate memory for # kmeans iterations/iso iteration  */

    c_data->kmean_iter = (int *) malloc((unsigned int)sizeof(int) * 10);
    if (c_data->kmean_iter == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,
                 "containing number of kmeans iterations per iso iteration\n");
        return(0);
    }

/*  allocate memory for max std dev array (dim = n_clusters) */

    i_parm->max_sd_dim = (int *) malloc((unsigned int)sizeof(int) * i_parm->max_cl_allowed);
    if (i_parm->max_sd_dim == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing dimension in which std dev is max\n");
        return(0);
    }

/*  allocate memory for merge map  */

    i_parm->merge_map = (int *) malloc((unsigned int)sizeof(int) * i_parm->max_cl_allowed);
    if (i_parm->merge_map == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing merge map\n");
        return(0);

    }

/*  allocate memory for merge map  */

    i_parm->shift_map = (int *) malloc((unsigned int)sizeof(int) * i_parm->max_cl_allowed);
    if (i_parm->shift_map == NULL )
    {
        (void)fprintf(stderr,"%s: cannot allocate memory for array", program);
        (void)fprintf(stderr,"containing shift_map\n");
        return(0);
    }



/* Allocate space for pairwise distance between c_centers */

    i_parm->dist_bet_c = (float **) malloc (sizeof (float *) * i_parm->max_cl_allowed);
    if (i_parm->dist_bet_c == NULL) {
        (void) fprintf (stderr,
                "%s: cannot allocate memory for pairwise distance array\n", program);
        return (0);
    }
    for (i = 0; i < i_parm->max_cl_allowed; i++) 
    {
        i_parm->dist_bet_c[i] = (float *) malloc (sizeof (float) * i_parm->max_cl_allowed);
        if (i_parm->dist_bet_c[i] == NULL) {
            (void) fprintf (stderr,
                    "%s: cannot allocate memory for pairwise distance array\n", program);
            return (0);
        }

     }

/* Allocate space for averages */

    c_data->avg = (double **) malloc (sizeof (double *) * i_parm->max_cl_allowed);
    if (c_data->avg == NULL) {
        (void) fprintf (stderr,
                "%s: cannot allocate memory for average array\n", program);
        return (0);
    }
    for (i = 0; i < i_parm->max_cl_allowed; i++) 
    {
        c_data->avg[i] = (double *) malloc (sizeof (double) * c_data->dimension);
        if (c_data->avg[i] == NULL) {
            (void) fprintf (stderr,
                    "%s: cannot allocate memory for average array\n", program);
            return (0);
        }

     }

/* Allocate space for variance */


    c_data->var = (double **) malloc (sizeof (double *) * i_parm->max_cl_allowed);
    if (c_data->var == NULL) {
        (void) fprintf (stderr,
                "%s: cannot allocate memory for variance array\n", program);
        return (0);
    }
    for (i = 0; i < i_parm->max_cl_allowed; i++) {
        c_data->var[i] = (double *) malloc (sizeof (double) * c_data->dimension);
        if (c_data->var[i] == NULL) {
            (void) fprintf (stderr,
                    "%s: cannot allocate memory for variance array\n", program);
            return (0);
        }
    }

/* Allocate space for cluster centers */

    c_data->c_center = (float **) malloc (sizeof (float *) * i_parm->max_cl_allowed);
    if (c_data->c_center == NULL) {
        (void) fprintf (stderr,
                "%s: insufficient memory available\n", program);
        return (0);
    }

    for (i = 0; i < i_parm->max_cl_allowed; i++) {
        c_data->c_center[i] = (float *) malloc (sizeof (float) * c_data->dimension);
        if (c_data->c_center[i] == NULL) {
            (void) fprintf (stderr,
                    "%s: insufficient memory available\n",program);
            return (0);
        }
    }

/* Allocate space for class data */

    c_data->class = (int *) malloc ((unsigned int) c_data->n_vects * sizeof (int));
    if (c_data->class == NULL) {
        (void) fprintf (stderr, "%s: insufficient memory available\n", program);
        return (0);
    }

    return(1);
}


/**************************************************************
*
* MODULE NAME:  kmean()
*
*     PURPOSE: 
*       The K-means algorithm is based on the minimization of the sum of
*       the squared distances from all points in a cluster to a cluster
*       center.  The K initial cluster centers and the image
*       vectors are iteratively distributed among the K cluster domains.
*       New cluster centers are computed from these results, such that
*       the sum of the squared distances from all points in a cluster to
*       the new cluster center is minimized.
*
*       INPUT:  c_data - clustering information structure.
*
*      OUTPUT: adjusted c_data structure.
*
* CALLED FROM:  iso2()
*
* ROUTINES CALLED: NONE
*
**************************************************************/

int
kmean (c_data, i_parm)

struct c_struct *c_data;
struct iso_params *i_parm;

{
    int c,i,j,iter,converge;
    int k;
    int zero_length;
    float a, d, dmin;

    iter = 0;

/*... START KMEANS LOOP ...*/

    do {

        /* reset parameters to 0 */
        zero_length = c_data->dimension * sizeof(double);

        bzero((char *) c_data->n_pts, (c_data->n_clusters * sizeof(float)));
        for (i = 0; i < c_data->n_clusters; i++) {
            bzero(c_data->avg[i], zero_length);
            bzero(c_data->var[i], zero_length);
        }

        /* assign each vector to a cluster */
        for (i = 0 ; i < c_data->n_vects; i ++ )
        {
           dmin = XV_MAXFLOAT;
           if (c_data->class[i] != OUT)
           {
                /* find the cluster to which it is closest */
                for (j = 0; j < c_data->n_clusters; j++) {

                    /* compute the euclidean distance */

                    d = 0.0;
                    for (k = 0; k < c_data->dimension; k++) {
                       a = *(c_data->input_vectors[i] + k) - *(c_data->c_center[j] +k); 
                       d += a * a;
                    }
                   /* test to see if we found a minimum */

                    if (d < dmin) {
                        dmin = d;
                        c = j;
                    }
                }

                /* assign vector to cluster, inc # points in that cluster 
                   and update new cluster center */                         

                c_data->class[i] = c;
                c_data->n_pts[c] += 1.0;
                for (j = 0; j < c_data->dimension; j++) {
                    *(c_data->avg[c] + j) += (double) *(c_data->input_vectors[i] + j);
                    *(c_data->var[c] + j) += (double) (*(c_data->input_vectors[i] + j) * 
                                               *(c_data->input_vectors[i] + j));
                }
           }
        }

        /* compute new cluster centers */
        for (i = 0; i < c_data->n_clusters; i++) {
            if (c_data->n_pts[i] > 0) {
                for (j = 0; j < c_data->dimension; j++) {
                    *(c_data->avg[i] + j) /= (double) c_data->n_pts[i];
                    *(c_data->var[i] + j) /= (double) c_data->n_pts[i];
                    *(c_data->var[i] + j) = *(c_data->var[i] + j) - (*(c_data->avg[i] + j) * *(c_data->avg[i] + j));
                }
            }
        }

        /* check for convergence */
        converge = TRUE;
        for (i = 0; i < c_data->n_clusters; i++) {
            for (j = 0; j < c_data->dimension; j++) {
                if (*(c_data->c_center[i] + j) != (float) *(c_data->avg[i] + j)) {
                    converge = FALSE;
                }
                *(c_data->c_center[i] + j) = (float) *(c_data->avg[i] + j);
            }
        }

        iter++;

    } while (!converge && (iter < c_data->max_n_iters_kmeans));

    if ( ((i_parm->iter_iso+1) % 10) == 0)
    {
       if(!(c_data->kmean_iter = (int *)realloc(c_data->kmean_iter,10 * sizeof(int))))
       {
          fprintf(stderr, "kmeans: Insufficient Memeory avaliable\n"); 
       }
    }
    c_data->kmean_iter[i_parm->iter_iso] = iter;
    return(1);
}


/**************************************************************
*
* MODULE NAME:  merge2()
*
*     PURPOSE: 
*       Compute pairwise distances between all cluster centers.
*       If the distance is less than a user specified merging
*       factor times the average of all the pairwise distances,
*       merge the two clusters and compute the new cluster center.
*
*       INPUT: c_data - clustering information structure
*              i_parm - iso2 paramenters structure
*
*      OUTPUT: adjusted c_data and i_parm.
*
* CALLED FROM: iso2()
*
* ROUTINES CALLED: shift()
*
**************************************************************/

merge2(c_data,i_parm)

    struct c_struct   *c_data;
    struct iso_params *i_parm;

{
    int   c,i,j,k,n,       /* loop variables */
          temp,
          count,           /* count clusters that qualify for merging */
          first,second;
    float min, a, d,
          avg_dist;        /* avg pairwise distance between c_centers */


    if (c_data->n_clusters/2 > 0 )
    {
        count = 0;
        temp = 0;
        avg_dist = 0.0;

        /* Set up array which gives merging status of clusters        */
        /* 0 = cannot be merged;    1 = can be merged  */

        bzero(i_parm->merge_stat, (c_data->n_clusters * sizeof(float)));

        /* calculate pair wise distances between all cluster centers  */
        /* and compute the average pairwise distance.  */

        n = 1;
        for (i=0 ; i < (c_data->n_clusters-1) ; i++)
        {
            for (j=n ; j < c_data->n_clusters ; j++)
            {
                d = 0.0;
                for (k=0 ; k < c_data->dimension ; k++)
                {
                    a = *(c_data->c_center[i] + k) - *(c_data->c_center[j] + k);
                    d += a*a;
                }
                *(i_parm->dist_bet_c[i]+j) = (float) sqrt((double) d);
                avg_dist += *(i_parm->dist_bet_c[i]+j);
                temp++;
            }
            n++;
        }
        avg_dist /= temp;


        /* count the number of clusters that meet merging requirement */

        n = 1;
        for (i=0 ; i < (c_data->n_clusters-1) ; i++)
        {
            for (j=n ; j < c_data->n_clusters ; j++)
            {
                if (*(i_parm->dist_bet_c[i]+j)<(i_parm->merge_factor*avg_dist))
                {
                    count += 1;
                    i_parm->merge_stat[i] = 1;
                    i_parm->merge_stat[j] = 1;
                }
            }
            n++;
        }
    
        /* if two or more clusters qualify, number clusters merged is */
        /* the smallest of the count and max clusters merged */

        if (count > 1)
        {
          /* set up the merge map */
            for ( c = 0; c < c_data->n_clusters; c++)
            {
              i_parm->merge_map[c] = c;
            }
      
            temp = c_data->n_clusters/2;
            if ( temp > count )
            {
                temp = count;
            }
    
            /* use the smallest pairwise distance sets first .........*/

            for (k=0 ; k<temp ; k++)
            {
                min = XV_MAXFLOAT;
                n = 1;
                for (i=0 ; i < (c_data->n_clusters-1) ; i++)
                {
                    for (j=n ; j < c_data->n_clusters ; j++)
                    {
                        if ( i_parm->merge_stat[i] == TRUE
                        &&   i_parm->merge_stat[j] == TRUE
                        &&   *(i_parm->dist_bet_c[i]+j) < min )
                        {
                             min = *(i_parm->dist_bet_c[i]+j);
                             first = i;
                             second = j;
                        }
                    }
                    n++;
                }

                /*  Once least distance is found, merge clusters and */
                /*  compute new cluster center  */

                for (j=0 ; j < c_data->dimension ; j++)
                {
                    *(c_data->c_center[first]+j) = 
                        ( c_data->n_pts[first]*(*(c_data->c_center[first]+j)) 
                        + c_data->n_pts[second]*(*(c_data->c_center[second]+j)))
                        / (c_data->n_pts[first] + c_data->n_pts[second]);
                }


                i_parm->merge_stat[first] = 0;
                i_parm->merge_stat[second] = 0;


                /*  Bookkeeping ......................................*/

                c_data->n_pts[first] += c_data->n_pts[second];
                c_data->n_pts[second] = 0;
                i_parm->merge_map[second] = first;
            }
            i_parm->change = TRUE;
            shift(c_data, i_parm);
            i_parm->merge_factor *= (1 - i_parm->merge_converge);
        }
    }
}



/**************************************************************
*
* MODULE NAME: split2()
*
*     PURPOSE: 
*       Splitting.  (relative standard deviation approach)
*       The standard deviation in each dimension of 
*       each cluster is computed and the average of these 
*       standard deviations calculated.  The maximum standard
*       deviation for each cluster is then compared against 
*       the average standard deviation times a user specified
*       scaling factor.  If the max std. dev. is greater than 
*       that value, the cluster is split in two in the dimension 
*       of maximum standard deviation.  Two other splitting 
*       conditions also exist: (1) the number of points in the
*       cluster have to satisfy a minimum condition and (2) the 
*       average distance between members of a cluster and their 
*       cluster center is less than this same average taken over 
*       all of the clusters.
*       A cluster can only be split once per isodata iteration.
*
*       INPUT: c_data - clustering information structure
*              i_parm - iso2 paramenters structure
*
*      OUTPUT: adjusted c_data and i_parm.
*
* CALLED FROM: iso2()
*
* ROUTINES CALLED: None
*
**************************************************************/

split2(c_data,i_parm)

    struct c_struct     *c_data;
    struct iso_params *i_parm;

{
    int   c,i,j,temp, k;        /* loop variables */
    float t, a, d,
          sd_avg;       /* avg std. dev. of clusters (all dimensions) */

                                

    if (c_data->n_clusters < i_parm->max_cl_allowed)  /*allocation requirement*/
    {
      /*
       * Compute the average distance of the samples in a cluster
       * from their corresponding cluster centers.  Then compute
       * the overall average dist. of samples from their centers.
       * (This is the measure that tou & gonzales use.  Might
       * consider looking at distance in each dimension??)
       */

        i_parm->d_overall_avg = 0.0;
        bzero(i_parm->d_avg, c_data->n_clusters * sizeof(float));
    
        for (c=0 ; c < c_data->n_clusters ; c++)
        {
            for (i=0 ; i < c_data->n_vects ; i++)
            {
                if (c_data->class[i] == c)
                {
                    d = 0.0;
                    for (k=0 ; k < c_data->dimension ; k++)
                    {
                      a = *(c_data->input_vectors[i] + k) - 
                                              *(c_data->c_center[c] + k);
                      d += a*a;
                    }
                    i_parm->d_avg[c] += (float) sqrt((double) d);
                }
            }
            i_parm->d_overall_avg += i_parm->d_avg[c];
            i_parm->d_avg[c] /= c_data->n_pts[c];
        }
        i_parm->d_overall_avg /= c_data->n_vects - i_parm->discarded_vects;

        sd_avg = 0;
        temp = c_data->n_clusters;

        bzero(i_parm->max_sd, c_data->n_clusters * sizeof(float));
        bzero(i_parm->max_sd_dim, c_data->n_clusters * sizeof(int));
        for (c=0 ; c < c_data->n_clusters ; c++)
        {
            bzero(i_parm->sd, c_data->dimension * sizeof(float));

            /* find standard deviation of cluster c in each dimension */

            for (i=0 ; i < c_data->n_vects ; i++)
            {
                if (c_data->class[i] == c)  
                {
                    for (j=0 ; j < c_data->dimension ; j++)
                    {
                        t = *(c_data->input_vectors[i]+j) - 
                                                    *(c_data->c_center[c]+j);
                        t *= t;
                        i_parm->sd[j] += t;
                    }
                }
            }
            for (j=0 ; j < c_data->dimension ; j++)
            {
                i_parm->sd[j] = sqrt(i_parm->sd[j]/c_data->n_pts[c]);

                sd_avg += i_parm->sd[j];       /* compute running average */

                if (i_parm->sd[j] > i_parm->max_sd[c])     /*find largest std.*/
                {                                  /*dev. for cluster */
                    i_parm->max_sd[c] = i_parm->sd[j];
                    i_parm->max_sd_dim[c] = j;
                }
            }
        }
        sd_avg /= (c_data->n_clusters * c_data->dimension);

/*  if splitting rules are satisfied, split the cluster in two ......... */

        c = 0;
        while ( temp < i_parm->max_cl_allowed && c < c_data->n_clusters)
        {
            if ( i_parm->max_sd[c] > (i_parm->split_factor * sd_avg)
            && i_parm->d_avg[c] >= i_parm->d_overall_avg
            && c_data->n_pts[c] > (2*(i_parm->min_pts_allowed+1)) )
            {
                i_parm->change = TRUE;
                for (j=0 ; j < c_data->dimension ; j++)
                {
                    *(c_data->c_center[temp]+j) = *(c_data->c_center[c]+j);
                }
                *(c_data->c_center[temp] + i_parm->max_sd_dim[c]) += 
                        i_parm->placement * i_parm->max_sd[c];
                *(c_data->c_center[c] + i_parm->max_sd_dim[c]) -= 
                        i_parm->placement * i_parm->max_sd[c];
                temp++;
            }
            c++;
        }
        c_data->n_clusters = temp;
        if ( i_parm->change == TRUE )
        {
            i_parm->split_factor *= (1 + i_parm->split_converge);
        }
    }
}

/**************************************************************
*
* MODULE NAME: shift()
*
*     PURPOSE: 
*       Starts with cluster zero.  Checks to see if cluster has 
*       zero points and, if so, shifts class[vector], n_pts[cluster],
*       c_center[cluster] accordingly.  Also decrements n_clusters.
*       NOTE: bases removal of a cluster on: n_pts[cluster] == 0
*
*       INPUT: c_data - clustering information structure
*              i_parm - iso2 paramenters structure
*
*      OUTPUT: adjusted c_data and i_parm.
*
* CALLED FROM: iso2()
*
* ROUTINES CALLED:
*
**************************************************************/

shift(c_data, i_parm)
struct c_struct     *c_data;
struct iso_params   *i_parm;
{
    int c,i,j;
    int move, new_pos;

    bzero(i_parm->shift_map, i_parm->max_cl_allowed * sizeof(int));

    move = 0;

    for (c=0 ; c < c_data->n_clusters ; c++)
    {
       if (c_data->n_pts[c] == 0)
       {
          move++;
          i_parm->shift_map[c] = OUT;
       }
       else 
       {
          new_pos = c - move;
          i_parm->shift_map[c] = new_pos;
          c_data->n_pts[new_pos] = c_data->n_pts[c];
          for (j=0 ; j < c_data->dimension ; j++)
          {
               *(c_data->c_center[new_pos]+j) = *(c_data->c_center[c]+j);
          }
       }
     }

     c_data->n_clusters -= move;

     for (i=0 ; i < c_data->n_vects ; i++)
     {
        if (c_data->class[i] != OUT)
           c_data->class[i] = 
                      i_parm->shift_map[i_parm->merge_map[c_data->class[i]]];
     }
}

/**************************************************************
*
* MODULE NAME:  discard
*
*     PURPOSE: 
*       description: isodata Rule 1:
*
*       Discard any cluster that does not have enough patterns
*       If cluster doesn't have enough points, it assigns the 
*       class of those points to an out of range class value (-1), 
*       and shifts the means, etc. of the other classes or vectors.
*
*       INPUT: c_data - clustering information structure
*              i_parm - iso2 paramenters structure
*
*      OUTPUT: adjusted c_data and i_parm.
*
* CALLED FROM: iso2()
*
* ROUTINES CALLED: shift
*
**************************************************************/

discard(c_data,i_parm)
struct c_struct     *c_data;
struct iso_params *i_parm;
{
    int c;
    int need_to_shift;
 
    need_to_shift = FALSE;

    for ( c=0 ; c < c_data->n_clusters ; c++ )
    {
        if (c_data->n_pts[c] < i_parm->min_pts_allowed )
        {
              /* need to shift is true */
            need_to_shift = TRUE;
            if (c_data->n_pts[c] > 0)
            {
               i_parm->change = TRUE;
               i_parm->discarded_vects += c_data->n_pts[c];
               c_data->n_pts[c] = 0;
            }
         }
      }
      if ( need_to_shift)
      {
             /* set up the merge map */
         for ( c = 0; c < c_data->n_clusters; c++)
         {
           i_parm->merge_map[c] = c;
         }

         shift(c_data,i_parm);
      }
}
/* -library_code_end */
