 /*
  * Khoros: $Id: lvkmeans.c,v 1.2 1991/12/18 09:25:02 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvkmeans.c,v 1.2 1991/12/18 09:25:02 dkhoros Exp $";
#endif

 /*
  * $Log: lvkmeans.c,v $
 * Revision 1.2  1991/12/18  09:25:02  dkhoros
 * HellPatch3
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * 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 to 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: lvkmeans.c
 >>>>
 >>>>      Program Name: vkmeans
 >>>>
 >>>> Date Last Updated: Wed Dec  4 14:44:08 1991 
 >>>>
 >>>>          Routines: lvkmeans - the library call for vkmeans
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvkmeans - library call for vkmeans
*
* Purpose:
*    
*    Performs the kmeans clustering algorithm on an image
*    
*    

* Input:
*    
*    img            input image structure.
*    
*    img1_flag      Output cluster centers in the "img1" image.
*    
*    sqd_img_flag   Output cluster variance in the "sqd_img" image.
*    
*    iter           maximum number of iterations.
*    
*    cluster        number of desired clusters.
*    
*    border         specifies number of border pixels.
*    
*    init_cluster   array of initial cluster centers.
*    
*    cc_flg         indicator to use cluster center image as input.
*    
*    cc_img         cluster center image structure.
*    
*    map_flag       If true output cluster center values in the map of
*                   the  cluster  number image. If false do not output
*                   cluster center values in the map  of  the  cluster
*                   number image.
*    
*    

* Output:
*    
*    img1           image structure for cluster center image.
*    
*    sqd_img        image structure for square distance values  (vari-
*                   ance).
*    
*    printdev       output ASCII file for cluster centers.
*    
*    This routine was written with the help of and ideas from Dr.  Don
*    Hush, University of New Mexico, Dept. of EECE.
*    
*    

*
* Written By: Tom Sauer, Charlie Gage
*    
*    

****************************************************************/


/* -library_def */
lvkmeans (img, img1, img1_flag, sqd_img, sqd_img_flag, iter, cluster, border, 
          init_cluster, printdev, cc_flg, map_flag, cc_img)
struct xvimage *img,            /* input image */
               *img1,           /* output image */
               *sqd_img,        /* square distance image */
               *cc_img;         /* cluster center image */
int     cluster,                /* number of desired clusters */
        iter;                   /* max number iterations */
int    *init_cluster[],         /* array of initial cluster centers x-y
                                   pairs */
        border,                 /* specifies number of border pixels */
        cc_flg,                 /* use cluster center as input */
        img1_flag,              /* flag for outputting cc image */
        sqd_img_flag,           /* flag for outputting variance image */
        map_flag;               /* flag for outputting cc in img map */
FILE * printdev;                /* 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 */
            nvects,             /* dummy number of vectors */
            dimension,          /* size of vector - num_bands */
            dim,                /* dummy dimension */
            vector,             /* index to input_vectors */
           *output_imagedata,   /* pointer to resulting image char */
            i,                  /* loop indices */
            j;

    float  *ptr,                /* pointer to the image data char */
          **input_vectors,      /* pointer to resulting image char */
          **c_center;           /* cluster centers */

    char **load_vector();
    char *map;
    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;



/* Allocate space for cluster centers */

    c_center = (float **) malloc (sizeof (float *) * cluster);
    if (c_center == NULL) {
        (void) fprintf (stderr,
                "lvkmeans: insufficient memory available\n");
        return (0);
    }

    for (i = 0; i < cluster; i++) {
        c_center[i] = (float *) malloc (sizeof (float) * dimension);
        if (c_center[i] == NULL) {
            (void) fprintf (stderr,
                    "lvkmeans: insufficient memory available\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, "lvkmeans: insufficient memory available\n");
        return (0);
    }

    /* Zero out output imagedata */

    bzero (output_imagedata, (unsigned int) vect * sizeof (float));


    input_vectors = (float **) load_vector(img, border, &nvects, &dim);

    /* release memory from input image */
    (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_center[i][j] = ptr[i + (cc_nc * cc_nr * j)];
            }
        }
        (void) free( cc_img -> imagedata);
    }
    else {              /* otherwise use input vectors */
        for (i = 0; i < cluster; i++) {
            vector = ((init_cluster[i][1] * (nc - (border * 2)))  + init_cluster[i][0]);
            for (j = 0; j < dimension; j++) {
                c_center[i][j] = input_vectors[vector][j];
            }
        }
    }

                        /* call kmeans routine */

    if (!kmeans (output_imagedata, img1, img1_flag, sqd_img, sqd_img_flag, 
                 input_vectors, cluster, iter, dimension, vect, c_center, 
                 printdev, nr, nc, border, &map, &map_flag)) 
    {
        (void) fprintf (stderr, "lvkmeans: kmeans Algorithm failed\n");
        return (0);
    }

       /* adjust header in the output image */
    img -> imagedata = (char *) output_imagedata;
    img -> data_storage_type = VFF_TYP_4_BYTE;
    img -> num_data_bands = 1;

     /*
      * if the map_flag is set then put the cluster center
      * values in the cluster number image's map. This
      * makes it spectrum compatable
      */

    if (map_flag == TRUE)
    {
       img -> maps = map;
       img -> map_storage_type = VFF_MAPTYP_FLOAT;
       img -> map_scheme = VFF_MS_ONEPERBAND;
       img -> map_row_size = dimension;
       img -> map_col_size = cluster;
       img -> map_subrow_size = 0;
       img -> ispare1 = 0;             /* to be compatable with spectrum */
       img -> ispare2 = dimension;      /* to be compatable with spectrum */
    }
    else
    {
       img -> maps = NULL;
       img -> map_storage_type = VFF_MAPTYP_NONE;
       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);
}


/******************************************************************

        K-Means Clusterng Algorithm 

        revision          date                  programmer
        --------          ----                  ----------
           0.0          08/11/88                Don Hush
           1.0          12/02/88                Tom Sauer
           1.1          07/18/89                Tom Sauer, Charlie Gage
           1.2          09/12/89                Tom Sauer

******************************************************************

        description :

        This routine implements the K-Means clustering algorithm.
    
        Revision 1.0
        Uses the same code for the algorithm, but the input and
        output are taken care of using khoros. The output
        will be an image. The maximum number of identifiable
        clusters is 256 (can only display 256 grey levels on
        an 16 plane machine). 

******************************************************************/

int
        kmeans (class, img1, img1_flag, sqd_img, sqd_img_flag, input_vectors, 
                n_clusters, max_n_iters, dimension, n_vects, c_center, 
                printdev, nr, nc, border, map, map_flag)

struct xvimage *img1, *sqd_img;

float **input_vectors;          /* pointer to resulting image char */
float **c_center;               /* c_center[i] ith cluster center */

int     n_clusters,             /* number of  initial cluster centers */
       *class,                  /* pointer to cluster center class */
        max_n_iters,            /* max number of iterations (user defined) */
        nr,
        nc,
        border,
        img1_flag,
        sqd_img_flag,
        *map_flag,
        n_vects;                /* size of image (nr * nc - border) */

char    **map;                  /* pointer to the map data for the cluster */
                                /* centers (map will be malloced */

FILE * printdev;                /* pointer to ASCII output stats file */
{

    /*---------- Program Variables ------------*/
    int     i,
            j,                  /* loop variables */
            k,
            x,
            z,
            count,
            iter,               /* iteration number */
            c,                  /* temporary class variable */
            zero_length,
            converge;           /* convergence flag */

    float  *n_pts,              /* n_pts[i] # points in cluster i */
           *sqdptr,
           *mapptr,
           *ptr,
            d, a,
            dmin,               /* distance & min distance */
            trace;
    double **avg,               /* avg[i] ith new cluster center value */
           **var;               /* var[i] ith cluster center value */

    float   dist ();


/*---------------------------START CODE------------------------------------*/


/*... allocate memory for data arrays ...*/


                        /* Allocate space for averages */

    avg = (double **) malloc (sizeof (double *) * n_clusters);
    if (avg == NULL) {
        (void) fprintf (stderr,
                "ERROR: cannot allocate memory for cluster centers\n");
        return (0);
    }

    for (i = 0; i < n_clusters; i++) {
        avg[i] = (double *) malloc (sizeof (double) * dimension);
        if (avg[i] == NULL) {
            (void) fprintf (stderr,
                    "ERROR: cannot allocate memory for cluster centers\n");
            return (0);
        }
    }

                        /* Allocate space for variances */

    var = (double **) malloc (sizeof (double *) * n_clusters);
    if (var == NULL) {
        (void) fprintf (stderr,
                "ERROR: cannot allocate memory for cluster centers\n");
        return (0);
    }

    for (i = 0; i < n_clusters; i++) {
        var[i] = (double *) malloc (sizeof (double) * dimension);
        if (var[i] == NULL) {
            (void) fprintf (stderr,
                    "ERROR: cannot allocate memory for cluster centers\n");
            return (0);
        }
    }

    n_pts = (float *) malloc (sizeof (float) * n_clusters);
    if (n_pts == NULL) {
        (void) fprintf (stderr,
                "ERROR: cannot allocate memory for data array\n");
        return (0);
    }


/*... initialize data arrays & variables ...*/


    iter = 0;

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

    do {

        /* reset parameters to 0 */
        zero_length = dimension * sizeof(double);
        for (i = 0; i < n_clusters; i++) {
            n_pts[i] = 0.0;
            bzero(avg[i], zero_length);
            bzero(var[i], zero_length);
        }

        /* assign each vector to a cluster */
        i = 0;
        for (x = border; x < nr - border; x++) {
            for (z = border; z < nc - border; z++) {

                dmin = XV_MAXFLOAT;

                /* find the cluster to which it is closest */
                for (j = 0; j < n_clusters; j++) {

                    /* compute the euclidean distance */

                    d = 0.0;
                    for (k = 0; k < dimension; k++) {
                       a = *(input_vectors[i] + k) - *(c_center[j] +k); 
                       d += a * a;
                    }
                    d = (float) sqrt ((double) d);

                   /* 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 */                         

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

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

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

        iter++;

    } while (!converge && (iter < max_n_iters));

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

    (void) fprintf (printdev, "vkmeans Statistics\n");
    (void) fprintf (printdev, "================\n\n");
    (void) fprintf (printdev, 
                  "\tTotal Number of K-MEANS Iterations : %d\n\n",iter);
    (void) fprintf(printdev, 
                    "\tTotal Number of Clusters = %d\n\n",n_clusters);

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

    (void) fprintf (printdev, "\n\nNumber of Vectors Per Cluster\n");
    (void) fprintf (printdev, "=============================\n\n");
    for  (i = 0; i< n_clusters; i++)
         (void) fprintf(printdev,"Cluster #%d = %d\n",i, (int)n_pts[i]);

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

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

    for (i = 0; i < n_clusters; i++) {
        count = 1;
        (void) fprintf (printdev, "\n\nCluster #%d\n\n", i);
        for (j = 0; j < dimension; j++) {
            (void) fprintf (printdev, "%f  ", 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", 
                   dimension);
    (void) fprintf (printdev, "==============================");
    for (i = 0; i < n_clusters; i++) {
        count = 1;
        (void) fprintf (printdev, "\n\nCluster #%d\n\n", i);

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


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

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

            }
        }
     }
 -------------------------------------------------------------------*/


    /*  
     *  update the cluster center image if the img1_flag is set
     *  update the cluster variance image if the sqdptr is set
     *  update the cluster number image`s map if the map_flag
     *  is set
     */

    if (img1_flag || sqd_img_flag || *map_flag )
    {
       if (img1_flag)
          ptr = (float *) img1->imagedata;
       if (sqd_img_flag)
          sqdptr = (float *) sqd_img->imagedata;
       if (map_flag)
       {
           *map = (char *)malloc (sizeof (float) * n_clusters * dimension);
           if (*map == NULL)
           {
             /* malloc space for the maps */ 
              fprintf(stderr, "lvkmeans: Warning, Cannot allocate enough space to store cluster centers in the \n");
              fprintf(stderr, "cluster number image map. The cluster number image will be written \n");
              fprintf(stderr, "out with without a map\n");
              *map_flag = FALSE;
           }
           mapptr = (float *) *map;
       }

       for (i = 0; i < n_clusters; i++) 
       {
           for (j = 0; j < dimension; j++) 
           {
               if (img1_flag)
               {
                 k = i + (j * img1->row_size * img1->col_size);
                 ptr[k] = c_center[i][j];
               }
               if (sqd_img_flag)
               {
                 k = i + (j * sqd_img->row_size * sqd_img->col_size);
                 sqdptr[k] = (float) var[i][j];
               }
               if (*map_flag)
                  mapptr[(j * n_clusters + i)] = c_center[i][j];
           }
       }
    }

    vfreelist(var, dimension);
    vfreelist(avg, dimension);
    


    return (1);
}
/* -library_code_end */
