 /*
  * Khoros: $Id: lvquant.c,v 1.3 1991/12/18 09:26:02 dkhoros Exp $
  */

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

 /*
  * $Log: lvquant.c,v $
 * Revision 1.3  1991/12/18  09:26: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: lvquant.c
 >>>>
 >>>>      Program Name: vquant
 >>>>
 >>>> Date Last Updated: Thu Nov 14 16:03:11 1991 
 >>>>
 >>>>          Routines: lvquant - the library call for vquant
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
struct vector
  {
    float *vec;                           /* Vector components */
    int count;                            /* Pixel count */
    struct vector *prev;                  /* Ptrs to previous & next vector */
    struct vector *next;
    struct pixel *pixlist;                /* List of pixels of this vector */
    struct pixel *pixtail;                /* End of pixel list */
  };

struct pixel
  {
    short x,y;                            /* Location of pixel */
    struct pixel *next;                   /* Pointer to next one */
  };

struct vector *vhead,*vtail;               /* Composite vector list */
struct vector *vh[65536],*vt[65536];       /* Hash index vector lists */

struct box
  {
    float *vmin,*vmax;                    /* Bounds on box */
    struct vector *vectors;               /* List of vectors in box */
    struct box *next;                     /* Next box */
    int count;                            /* Number of pixels in box */
  };

struct box *boxhead,*boxtail;

int ndim;                                 /* Dimensionality of space */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvquant - library call for vquant
*
* Purpose:
*    
*    Perform N-dimensional vector quantization or classification
*    
*    

* Input:
*    
*    i              pointer to vector image to be quantized.
*    
*    nvecs          number of vectors to generate
*    
*    mapflg         output image map enable
*    
*    ccimage        cluster center output image structure
*    
*    cvimage        cluster variance output image structure
*    
*    sp             split point method: 0 is mid-span, 1 is mean
*    
*    axis           split axis: 1 is max-span axis, 2 is max  variance
*                   axis, 3 is principal eigenvector
*    
*    

* Output:
*    
*    i              the image pointed to by i on input is modiifed  to
*                   contain the single-band classification pixel plane
*                   and the map is modified so that  it  contains  the
*                   representative vectors.
*    
*    

*
* Written By: Scott Wilson
*    
*    23-Jul-91 Scott Wilson Fixed uninitialized next pointer (Tom
*    Lane)
*    
*    

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


/* -library_def */
int
lvquant(i,nvecs,mapflg,ccimage,ccflg,cvimage,cvflg,sp,axis)
struct xvimage *i,**ccimage,**cvimage;
int nvecs,mapflg,ccflg,cvflg,sp,axis;
/* -library_def_end */

/* -library_code */
{
  ndim = i->num_data_bands;            /* Set the dimensionality */
  load_vector_list(i);                 /* Make vector list */
  compress_vectors(nvecs,sp,axis);     /* Squeeze the vector list */
  quant_adjust_image(i,nvecs,mapflg,ccimage,ccflg,cvimage,cvflg); /* quantize */
/* Note: The "free" routine is commented out because it takes
         durn near as long for the system to take the data structure
         back as it does to do the entire quantization!  If you
         need to use this routine in a program, you'll want to uncomment
         the "free" call.    SRW */
/*    free_vecs();  */                      /* Give back all memory */

  return(1);                          /* success, BUT other routines exit() */
}

load_vector_list(i1)                      /* Build the vector list */
struct xvimage *i1;
  {
    int i,j,k,l,len;
    float *vec;

    vhead = NULL;                         /* Initialize vector list */
    vtail = NULL;
    for (i=0; i<65536; i++)               /* Initialize hash index lists */
      {
        vh[i] = NULL; vt[i] = NULL;
      }

    vec = (float *)malloc(ndim*sizeof(float)); /* Work vector */
    if (vec == NULL)
      {
        (void)fprintf(stderr,"lvquant: Not enough memory for work vector!\n");
        exit(1);
      }

    len = i1->row_size*i1->col_size;     /* Number of pixels per plane */
    for (j=0; j<i1->col_size; j++)
      {
        for (i=0; i<i1->row_size; i++)
          {
            l = j*i1->row_size+i;              /* Index of pixel in plane 0 */
            for (k=0; k<ndim; k++)
              vec[k] = *((float *)(i1->imagedata)+k*len+l);
            update_vlist(vec,i,j);
          }
      }

    /* Now append all vector lists together */
    for (i=0; i<65536; i++)
      {
        if (vh[i] != NULL)     /* If have something to add */
          {
            if (vtail != NULL) /* If main list not empty */
              {
                vtail->next = vh[i];
                vh[i]->prev = vtail;
                vtail = vt[i];
              }
            else               /* Main list is empty */
              {
                vhead = vh[i];
                vtail = vt[i];
              }
          }
      }
  }

update_vlist(vec,x,y)
float *vec;
int x,y;
  {
    struct vector *pv;

    find_vector(vec,&pv);
    if (pv == NULL) add_vector(vec,x,y);
    else update_vector(pv,x,y);
  }

find_vector(vec,pv)
float *vec;
struct vector **pv;
  {
    int i,k,found;
    struct vector *p;

    k = f(vec);         /* Compute hash index */
    p = vh[k];
    while (p != NULL)
      {
        found = 1;
        for (i=0; i<ndim; i++) if (p->vec[i] != vec[i]) found = 0;
        if (found) break;
        else p = p->next;
      }
    *pv = p;
  }

add_vector(vec,x,y)
float *vec;
int x,y;
  {
    int k;
    struct vector *p;

    p = (struct vector *)malloc(sizeof(struct vector));
    p->vec = (float *)malloc(ndim*sizeof(float));
    if (p == NULL || p->vec == NULL)
      {
        (void)fprintf(stderr,"lvquant: Not enough memory for additional vector!\n");
        exit(1);
      }
    k = f(vec);
    if (vh[k] == NULL && vt[k] == NULL) /* Empty vector list */
      {
        vh[k] = p;
        vt[k] = p;
        p->prev = NULL;
        p->next = NULL;
      }
    else                                /* Not empty, add to head */
      {
        p->prev = NULL;
        p->next = vh[k];
        vh[k]->prev = p;
        vh[k] = p;
      }
    for (k=0;  k<ndim; k++) p->vec[k] = vec[k];
    p->count = 1;
    p->pixlist = (struct pixel *)malloc(sizeof(struct pixel));
    if (p->pixlist == NULL)
      {
        (void)fprintf(stderr,"lvquant: Not enough memory for additional pixel!\n");
        exit(1);
      }
    p->pixlist->next = NULL;
    p->pixlist->x = x;
    p->pixlist->y = y;
    p->pixtail = p->pixlist;
  }

update_vector(p,x,y)
struct vector *p;
int x,y;
  {
    int k;

    k = f(p->vec);
    p->count++;
    p->pixtail->next = (struct pixel *)malloc(sizeof(struct pixel));
    if (p->pixtail->next == NULL)
      {
        (void)fprintf(stderr,"lvquant: Not enough memory for additional pixel!\n");
        exit(1);
      }
    p->pixtail = p->pixtail->next;
    p->pixtail->x = x;
    p->pixtail->y = y;
    p->pixtail->next = NULL;
    if (p != vh[k]) quant_head_swap(p);
  }

quant_head_swap(p)
struct vector *p;
  {
    int k;

    k = f(p->vec);

    /* Delete current vector from list */
    if (p == vt[k])
      {
        p->prev->next = p->next;
        vt[k] = p->prev;
      }
    else
      {
        p->prev->next = p->next;
        p->next->prev = p->prev;
      }

    /* Re-insert vector at head of list*/
    p->next = vh[k];
    p->prev = NULL;
    vh[k]->prev = p;
    vh[k] = p;
  }
 
compress_vectors(nvecs,sp,axis)
int nvecs,sp,axis;
  {
    int i,j,n,k;
    float vmax,vmean,mean,var;
    double *m,*pe,*mv,inner(),mu,mtest;
    struct vector *p1,*p2;
    struct box *b1,*b2,*find_biggest_box();

    /* Initialize box structure */
    boxhead = NULL;
    boxtail = NULL;

    /* Make first box and attach all vectors to it */
    b1 = (struct box *)malloc(sizeof(struct box));
    b1->vmin = (float *)malloc(ndim*sizeof(float));
    b1->vmax = (float *)malloc(ndim*sizeof(float));
    if (b1 == NULL || b1->vmin == NULL || b1->vmax == NULL)
      {
        (void)fprintf(stderr,"lvquant: Not enough memory for root box!\n");
        exit(1);
      }
    boxhead = b1;
    boxtail = b1;
    b1->vectors = vhead;
    quant_set_box_bounds(b1);
    b1->next = NULL;

    /* Now go through and perform a median cut on the box list */
    /* until the desired number of boxes (vectors) are formed. */
    for (i=0; i<nvecs-1; i++)
      {
        /* Hunt down the biggest box */
        b1 = find_biggest_box();
        if (b1 == NULL)
          {
            (void)fprintf(stderr,"lvquant: Could not find splittable vector subspace!\n");
            exit(1);
          }

        /* First make new box and attach to box list */
        b2 = (struct box *)malloc(sizeof(struct box));
        b2->vmin = (float *)malloc(ndim*sizeof(float));
        b2->vmax = (float *)malloc(ndim*sizeof(float));
        if (b2 == NULL || b2->vmin == NULL || b2->vmax == NULL)
          {
            (void)fprintf(stderr,"lvquant: Not enough memory for additional box!\n");
            exit(1);
          }
        b2->next = NULL;
        b2->vectors = NULL;
        boxtail->next = b2;
        boxtail = boxtail->next;

        /* Copy the vector bounds to the new box. */
        for (j=0; j<ndim; j++) 
          {
            b2->vmax[j] = b1->vmax[j];
            b2->vmin[j] = b1->vmin[j];
          }

        /* Find the axis for splitting the space */
        if (axis == 0 || axis == 1)
          {
            /* Split on axis with largest span */
            vmax = b1->vmax[0]-b1->vmin[0];
            n = 0;
            for (j=1; j<ndim; j++)
              {
                if (b1->vmax[j]-b1->vmin[j] > vmax)
                  {
                    vmax = b1->vmax[j]-b1->vmin[j];
                    n = j;
                  }
              }
          }
        else if (axis == 2)
          {
            /* Split on axis with largest variance */
            vmax = -1;
            vmean = 0;
            n = 0;
            for (j=0; j<ndim; j++)
              {
                mean = 0.0; k = 0;
                p1 = b1->vectors;
                while (p1 != NULL)
                  {
                    mean += p1->vec[j];
                    p1 = p1->next;
                    k++;
                  }
                mean /= (float)k;
                var = 0.0;
                p1 = b1->vectors;
                while (p1 != NULL)
                  {
                    var += (p1->vec[j]-mean)*(p1->vec[j]-mean);
                    p1 = p1->next;
                  }
                var /= (float)k;
                if (var > vmax)
                  {
                    vmax = var;
                    vmean = mean;
                    n = j;
                  }
              }
            }
          else
            {
              /* Split on principal eigenvector */
              /* Obtain this object by constructing the covariance
                 matrix for the vectors and obtaining the eigenvector
                 associated with the largest eigenvalue of this matrix */
              /* First obtain mean vector */
              mv = (double *)malloc(ndim*sizeof(double));
              bzero(mv,ndim*sizeof(double));
              k = 0;
              p1 = b1->vectors;
              while (p1 != NULL)
                {
                  for (j=0; j<ndim; j++) mv[j] += p1->vec[j];
                  p1 = p1->next;
                  k++;
                }
              for (j=0; j<ndim; j++) mv[j] = mv[j] / k;
              /* Now compute covariance matrix */
              m = (double *)malloc(ndim*ndim*sizeof(double));
              bzero(m,ndim*ndim*sizeof(double));
              p1 = b1->vectors;
              while (p1 != NULL)
                {
                  for (j=0; j<ndim; j++)
                    for (k=0; k<ndim; k++)
                      m[j*ndim+k] = (p1->vec[k]-mv[k]) * (p1->vec[j]-mv[j]);
                  p1 = p1->next;
                }
              pe = (double *)malloc(ndim*sizeof(double));
              get_princ_axis(m,ndim,pe);
              free(m);
            }
        
        /* Partition the vector space of the first box */
        p1 = b1->vectors;           /* Vector list in original box */
        p2 = b2->vectors;           /* End of vector list in new box */
        if ((axis == 1 || axis == 2) && sp == 0)
          {
            b1->vmax[n] = (b1->vmax[n]+b1->vmin[n])/2;
            b2->vmin[n] = b1->vmax[n];
            while (p1 != NULL)
              {
                if (p1->vec[n] > b1->vmax[n]) box_move_vector(b1,b2,&p1,&p2);
                else p1 = p1->next;
              }
          }
        else
          {
            if (axis == 1 || axis == 2)
              {
                if (axis==1) b1->vmax[n] = vmean;
                else b1->vmax[n] = (b1->vmax[n]+b1->vmin[n])/2;
                b2->vmin[n] = b1->vmax[n];
                while (p1 != NULL)
                  {
                    if(p1->vec[n] > b1->vmax[n]) box_move_vector(b1,b2,&p1,&p2);
                    else p1 = p1->next;
                  }
              }
            else
              {
                /* Split along principal eigenvalue axis */
                /* Split the population by checking against the
                   inner product of the mean and p-evect .vs. the
                   inner product of the vector and the p-evect */
                mu = inner(mv,pe,ndim);
                free(mv);
                p1 = b1->vectors;           /* Vector list in original box */
                p2 = b2->vectors;           /* End of vector list in new box */
                mv = (double *)malloc(ndim*sizeof(double));
                while (p1 != NULL)
                  {
                    for (j=0; j<ndim; j++) mv[j]=p1->vec[j]; /* To DOUBLE */
                    mtest = inner(mv,pe,ndim);
                    if (mtest > mu) box_move_vector(b1,b2,&p1,&p2);
                    else p1 = p1->next;
                  }
                free(pe);
                free(mv);
              }
          }

        /* Re-compute the vector bounds for both boxes */
        quant_set_box_bounds(b1);
        quant_set_box_bounds(b2);
      }
  }

box_move_vector(b1,b2,p1,p2)
struct box *b1,*b2;
struct vector **p1,**p2;
  {
    /* Move the vector p1 from box b1's vector list to box b2's vector list,
       which has the tail of its vector list marked by p2. p1 must me left
       pointing to the next vector in b1's vector list, which p2 must be left
       pointing at the last vector in b2's vector list. */

    struct vector *v1;

    v1 = (*p1)->next;           /* Remember where "next" is. */
    unlink_vector(b1,p1);
    append_vector(b2,p2,p1);
    *p1 = v1;                   /* Regurgitate "next". */
  }

unlink_vector(b,v)
struct box *b;
struct vector **v;
  {
    if (b->vectors == *v && (*v)->next != NULL)     /* v is head of vector list */
      {
        b->vectors = (*v)->next;
        (*v)->next->prev = NULL;
      }
    else if (b->vectors == *v && (*v)->next == NULL) /* v is only vector */
      {
        b->vectors = NULL;
      }
    else if (b->vectors != *v && (*v)->next == NULL) /* v is tail of list */
      {
        (*v)->prev->next = NULL;
      }
    else                                            /* v in middle of list */
      {
        (*v)->next->prev = (*v)->prev;
        (*v)->prev->next = (*v)->next;
      }
    (*v)->next = NULL;
    (*v)->prev = NULL;
  }

append_vector(b,v,vnew)
struct box *b;
struct vector **v,**vnew;
  {
    if (b->vectors == NULL)                          /* If empty vector list */
      {
        b->vectors = *vnew;
        *v = *vnew;
      }
    else                                            /* Non-empty vectors list */
      {
        (*v)->next = *vnew;
        (*vnew)->prev = *v;
        *v = *vnew;
      }
  }

struct box *find_biggest_box()
 {
   /* Find the box that represents the largest subspace in the 2-norm
      sense. */
   struct box *b,*b1;
   float max,sum,x;
   int i;

   b = boxhead;
   b1 = NULL;
   max = 0.0;
   while (b != NULL)
     {
       sum = 0.0;
       for (i=0; i<ndim; i++) /* Accumulate 2-norm of subspace */
         {
           x = b->vmax[i]-b->vmin[i];
           sum += x*x;
         }
       if (sum > max)
         {
           max = sum;
           b1 = b;
         }
       b = b->next;
     }
   return(b1);
 }

quant_set_box_bounds(b)
struct box *b;
  {
    struct vector *v;
    int i,n;

    for (i=0; i<ndim; i++)
      {
        b->vmax[i] = -1.0e38;
        b->vmin[i] =  1.0e38;
      }

    n = 0;
    v = b->vectors;
    while (v != NULL)
      {
        for (i=0; i<ndim; i++)
          {
            if (v->vec[i] > b->vmax[i]) b->vmax[i] = v->vec[i];
            if (v->vec[i] < b->vmin[i]) b->vmin[i] = v->vec[i];
          }
        v = v->next;
        n++;
      }
    b->count = n;
  }

quant_adjust_image(image,nvecs,mapflg,ccimage,ccflg,cvimage,cvflg)
struct xvimage *image,**ccimage,**cvimage;
int nvecs,mapflg,ccflg,cvflg;
  {
    int i,j,k;
    static float *vec;
    float *f;
    unsigned int *d;
    struct box *bp;
    struct vector *v;
    struct pixel *p;
    struct xvimage *createimage();

    /* Grab space for the vector sum array */
    if (vec == NULL) vec = (float *)malloc(ndim*sizeof(float));

    /* Free up the existing image data area and allocate a new one of
       type int.  Then kill off the old map (if there was one) and set up for
       a new one to contain the quantized vectors. */
    free(image->imagedata);
    image->imagedata = malloc(image->row_size*image->col_size*sizeof(unsigned int));
    if (image->imagedata == NULL)
      {
        (void)fprintf(stderr,"lvquant: Unable to allocate memory for new image!\n");
        exit(1);
      }
    image->data_storage_type = VFF_TYP_4_BYTE;
    image->num_data_bands = 1;
    if (mapflg)
      {
        /* Output map desired */
        image->map_row_size = ndim;
        image->map_col_size = nvecs;
        image->map_storage_type = VFF_MAPTYP_FLOAT;
        free(image->maps);
        image->maps = malloc(nvecs*ndim*sizeof(float));
        if (image->maps == NULL)
          {
            (void)fprintf(stderr,"lvquant: Not enough memory for maps!\n");
            exit(1);
          }
        image->map_enable = VFF_MAP_OPTIONAL;
        image->map_scheme = VFF_MS_ONEPERBAND;
        image->color_space_model = VFF_CM_NONE;
        image->map_subrow_size = image->subrow_size; /* Transfer subrow info */
        image->subrow_size = 0;
        image->ispare1 = 0;     /* for compatability with spectrum */
        image->ispare2 = ndim;  /* for compatability with spectrum */
      }
    else
      {
        /* No output map is desired */
        image->map_row_size = 0;
        image->map_col_size = 0;
        image->map_storage_type = VFF_MAPTYP_NONE;
        free(image->maps);
        image->map_enable = VFF_MAP_OPTIONAL;
        image->map_scheme = VFF_MS_NONE;
        image->color_space_model = VFF_CM_NONE;
        image->map_subrow_size = 0;
        image->subrow_size = 0;
      }

    /* Read thru the box list stepping thru its vector list and placing
       the pixels in the output image */
    d = (unsigned int *)(image->imagedata);
    i = 0;
    bp = boxhead;
    while (bp != NULL)
      {
        v = bp->vectors;
        while (v != NULL)
          {
            p = v->pixlist;
            while (p != NULL)
              {
                d[p->y*image->row_size + p->x] = i;
                p = p->next;
              }
            v = v->next;
          }
        bp = bp->next;
        i++;
      }
 
    if (mapflg || ccflg || cvflg)
      {
        /* Compute the vector for each box as the weighted average of the
           vector components inside the box */
        bp = boxhead;
        while (bp != NULL)
          {
            i = 0;
            bzero(vec,ndim*sizeof(float));
            v = bp->vectors;
            while (v != NULL)
              {
                for (k=0; k<ndim; k++) vec[k] += v->count*v->vec[k];
                i += v->count;
                v = v->next;
              }
            for (k=0; k<ndim; k++) bp->vmin[k] = vec[k]/i;
            bp = bp->next;
          }
      }

    if (cvflg)
      {
        /* Compute the variance vector for each box as the variance of the
           vector components inside the box */
        bp = boxhead;
        while (bp != NULL)
          {
            i = 0;
            bzero(vec,ndim*sizeof(float));
            v = bp->vectors;
            while (v != NULL)
              {
                for (k=0; k<ndim; k++) vec[k] += v->count*v->vec[k]*v->vec[k];
                i += v->count;
                v = v->next;
              }
            for (k=0; k<ndim; k++) bp->vmax[k] = vec[k]/i;
            bp = bp->next;
          }
      }

    if (mapflg)
      {
        /* Put the vectors in the map */
        j = 0;
        bp = boxhead;
        while (bp != NULL)
          {
            for (i=0; i<ndim; i++)
              *((float *)(image->maps)+i*image->map_col_size+j) = bp->vmin[i];
            bp = bp->next;
            j++;
          }
      }

    if (ccflg)
      {
        /* Create an image structure to hold the cluster centers, fill
           it in, and set the pointer to point at this image */
        *ccimage = createimage(1,nvecs,
                        (unsigned long)
                        VFF_TYP_FLOAT,         /* data_storage_type */
                        (unsigned long) 1,      /* num_of_images */
                        (unsigned long) ndim,      /* num_data_bands */
                        "Cluster center image",        /* comment */
                        (unsigned long) 0,      /* map_row_size */
                        (unsigned long) 0,      /* map_col_size */
                        (unsigned long)
                        VFF_MS_NONE,             /* map_scheme */
                        (unsigned long)
                        VFF_MAPTYP_NONE,        /* map_storage_type */
                        (unsigned long)
                        VFF_LOC_IMPLICIT,       /* location_type */
                        (unsigned long) 0);     /* location_dim */

        if (ccimage == NULL)
          {
            (void)fprintf(stderr,"vquant: Unable to create CC image!\n");
            exit(1);
          }

        f = (float *)((*ccimage)->imagedata);
        /* Put the cluster center vectors in the image */
        bp = boxhead;
        while (bp != NULL)
          {
            for (i=0; i<ndim; i++)
              *(f+i*nvecs) = bp->vmin[i];
            bp = bp->next;
            f++;
          }
      }
        
    if (cvflg)
      {
        /* Create an image structure to hold the cluster centers, fill
           it in, and set the pointer to point at this image */
        *cvimage = createimage(1,nvecs,
                        (unsigned long)
                        VFF_TYP_FLOAT,         /* data_storage_type */
                        (unsigned long) 1,      /* num_of_images */
                        (unsigned long) ndim,      /* num_data_bands */
                        "Cluster variance image",        /* comment */
                        (unsigned long) 0,      /* map_row_size */
                        (unsigned long) 0,      /* map_col_size */
                        (unsigned long)
                        VFF_MS_NONE,             /* map_scheme */
                        (unsigned long)
                        VFF_MAPTYP_NONE,        /* map_storage_type */
                        (unsigned long)
                        VFF_LOC_IMPLICIT,       /* location_type */
                        (unsigned long) 0);     /* location_dim */

        if (cvimage == NULL)
          {
            (void)fprintf(stderr,"vquant: Unable to create CV image!\n");
            exit(1);
          }

        f = (float *)((*cvimage)->imagedata);
        /* Put the cluster variance vectors in the image */
        bp = boxhead;
        while (bp != NULL)
          {
            for (i=0; i<ndim; i++)
              *(f+i*nvecs) = bp->vmax[i]-(bp->vmin[i]*bp->vmin[i]);
            bp = bp->next;
            f++;
          }
      }
  }

free_vectors()
  {
    /* Release all malloc()'ed items except the new colormap. */
    struct box *b,*b1;
    struct vector *v,*v1;
    struct pixel *p,*p1;

    b = boxhead;
    while (b != NULL)
      {
        v = b->vectors;
        while (v != NULL)
          {
            p = v->pixlist;
            while (p != NULL)
              {
                p1 = p->next;
                free((char *)p);
                p = p1;
              }
            v1 = v->next;
            free((char *)v);
            v = v1;
          }
        b1 = b->next;
        free((char *)b);
        b = b1;
      }
  }

f(vec)
float *vec;
  {
    /* Compute a hash-index function that is 16 bits long from the
       components of the supplied vector. This is done by summing
       the set of 16 bit short ints formed by splitting each float
       into its high and low 16-bit words. */
    int i;
    short *s;
    int sum;

    s = (short *)vec;
    sum = 0;
    for (i=0; i<2*ndim; i++) sum += *s++;

    return(sum & 0xffff);
  }
/* -library_code_end */
