 /*
  * Khoros: $Id$
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id$";
#endif

 /*
  * $Log$
  */

/*
 *----------------------------------------------------------------------
 *
 * 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: lvrdclust.c
 >>>>
 >>>>      Program Name: vrdclust
 >>>>
 >>>> Date Last Updated: Wed Aug  7 09:39:02 1991 
 >>>>
 >>>>          Routines: lvrdclust - the library call for vrdclust
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */

#define MAX_FINAL_NUM 256

/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvrdclust - library call for vrdclust
*
* Purpose:
*
* Input:
* Output:
* Written By:
*
****************************************************************/


/* -library_def */

int lvrdclust (in_image,out_image,final_num,alg,scale,global_scale,calc_hist,mapping)

struct xvimage *in_image;          /* The input image */
struct xvimage **out_image;        /* The output image */
int final_num;                     /* Final number of clusters */
char *alg;                         /* Algorithm to use */
int scale;                         /* Scale mean vectors? */
int global_scale;                  /* Global Scaling? */
int calc_hist;                     /* Calculate Histogram? */
int **mapping;                     /* (OUTPUT) Mapping for the clusters */

/* -library_def_end */

/* -library_code */

{
     int num_entries;      /* Current number of mean vectors */
     float **means;        /* Mean vectors */
     int *count;           /* Number of elements */
     int *class;           /* Class assignment number */
     float ***cov_upper;   /* Upper triangular covariance matrix */
     float **cov_diag;     /* Diagonal covariance matrix */
     int r, c;             /* Row and column indices */
     int cluster;          /* Cluster index */
     float min, max;       /* Minimum, Maximum */
     int i, k;             /* Loop control */
     float *tmp_means;     /* For free'ing memory */
     float *tmp_cov_diag;  /* For free'ing memory */

     int *mapping_empty;     

     unsigned char *ptr_out;             /* TEMPORARY */
     unsigned char *ptr_1_BYTE;          /* TEMPORARY */
     unsigned short *ptr_2_BYTE;         /* TEMPORARY */
     int *ptr_4_BYTE;                    /* TEMPORARY */
     double *ptr_DOUBLE;                 /* TEMPORARY */

     extern float **contig_matrix();

     *mapping = NULL;
     mapping_empty = NULL;

     /*  CHECK NUMBER OF FINAL CLUSTERS  */

     if ( (final_num <= 0) || (final_num > MAX_FINAL_NUM) ) {
          fprintf(stderr,"Final number must be between 1 and %d.\n",
                    MAX_FINAL_NUM);
          return (FALSE);
     }

     /*  GET CLUSTER STATISTICS OF INPUT IMAGE  */

     means = NULL;
     count = NULL;
     class = NULL;
     cov_upper = NULL;
     cov_diag = NULL;

     num_entries = in_image->map_col_size;

     if (!in_image->maps){
          fprintf (stderr,"No map data present.\n");
          return (FALSE);
     }

     if (in_image->map_storage_type != VFF_MAPTYP_FLOAT) {

          fprintf(stderr,"Map is not floating point data:  assuming only\n");
          fprintf(stderr,"mean vectors for clusters exist\n");

          means = contig_matrix ( (int) in_image->map_col_size, 
                    (int) in_image->map_row_size );

          switch (in_image->map_storage_type) {
               case VFF_MAPTYP_1_BYTE:
                    ptr_1_BYTE = (unsigned char *) in_image->maps;
                    for (r=0;r<in_image->map_col_size;r++){
                    for (c=0;c<in_image->map_row_size;c++){
                         means [r][c] = 
                              (float) ptr_1_BYTE [r+c*(in_image->map_col_size)];
                    }
                    }
                    break;
               case VFF_MAPTYP_2_BYTE:
                    ptr_2_BYTE = (unsigned short *) in_image->maps;
                    for (r=0;r<in_image->map_col_size;r++){
                    for (c=0;c<in_image->map_row_size;c++){
                         means [r][c] = 
                              (float) ptr_2_BYTE [r+c*(in_image->map_col_size)];
                    }
                    }
                    break;
               case VFF_MAPTYP_4_BYTE:
                    ptr_4_BYTE = (int *) in_image->maps;
                    for (r=0;r<in_image->map_col_size;r++){
                    for (c=0;c<in_image->map_row_size;c++){
                         means [r][c] = 
                              (float) ptr_4_BYTE [r+c*(in_image->map_col_size)];
                    }
                    }
                    break;
               case VFF_MAPTYP_DOUBLE:
                    ptr_DOUBLE = (double *) in_image->maps;
                    for (r=0;r<in_image->map_col_size;r++){
                    for (c=0;c<in_image->map_row_size;c++){
                         means [r][c] = 
                              (float) ptr_DOUBLE [r+c*(in_image->map_col_size)];
                    }
                    }
                    break;
               default:
                    fprintf (stderr,"Cannot handle map data type <%d>\n",
                         in_image->map_storage_type);
                    return (FALSE);
          }

     } else if ( ! getstats (in_image->maps,in_image->ispare1,
          (int)in_image->ispare2, (int)in_image->map_col_size,&means,&count,
          &class,&cov_upper, &cov_diag) ) {

          fprintf (stderr,"Could not get cluster statistics from getstats()\n");
          return (FALSE);

     }

     /*  CALCULATE THE HISTOGRAM WHEN ASKED  */

     if ( calc_hist ) {

          /*  ALLOCATE MEMORY WHEN NEEDED  */

          if ( count == NULL ) {
               count = (int *) malloc (sizeof(int)* in_image->map_col_size);
               if ( count == NULL ) {
                    fprintf (stderr,"MEMORY ALLOCATION FAILURE\n");
                    return(FALSE);
               }
          }

          /*  INITIALIZE  */

          for ( i = 0 ; i < in_image->map_col_size ; i ++ )  count [i] = 0;

          /*  CALCULATE  */

          switch ( in_image->data_storage_type ) {
               case VFF_TYP_1_BYTE:
                    ptr_1_BYTE = (unsigned char *) in_image->imagedata;
                    for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                         count [(int)ptr_1_BYTE[k]] ++ ;
                    }
                    break;
               case VFF_TYP_2_BYTE:
                    ptr_2_BYTE = (unsigned short *) in_image->imagedata;
                    for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                         count [(int)ptr_2_BYTE[k]] ++ ;
                    }
                    break;
               case VFF_TYP_4_BYTE:
                    ptr_4_BYTE = (int *) in_image->imagedata;
                    for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                         count [(int)ptr_4_BYTE[k]] ++ ;
                    }
                    break;
               default:
                    fprintf (stderr,"Image data not char, short, or int...\n");
                    fprintf (stderr,"cannot continue\n");
                    return(FALSE);
          }

          /*  DELETE EMPTY CLUSTERS  */

          delete_empties ( &num_entries, means, count, cov_upper, cov_diag, 
               &mapping_empty );

     }

     /*  CLASS ASSIGNMENTS WILL NO LONGER BE VALID  */

     if ( class ) {
          free (class);
          class = NULL;
     }


     /*  IF SCALING WILL OCCUR, DON'T DEAL WITH COVARIANCE MATRICES  */

     if ( scale && ( (alg [0] == '1') || (alg [0] == '2') ) ) {

          if ( cov_upper ) {
               for ( r = 0 ; r < num_entries ; r ++ ) {
                    free ( cov_upper [r][0] );
                    free ( cov_upper [r] );
               }
               free ( cov_upper );
               cov_upper = NULL;
          }

          if ( cov_diag ) {
               for ( r = 0 ; r < num_entries ; r ++ ) {
                    free ( cov_diag [r] );
               }
               free ( cov_diag );
               cov_diag = NULL;
          }

     }

     /*  
      *  Since deleting clusters will re-arrange memory, and since memory
      *  is always allocated in contiguous chunks, keep track of some
      *  pointers.
      */
     tmp_means = NULL;
     tmp_cov_diag = NULL;
     if ( means != NULL ) tmp_means = means [0];
     if ( cov_diag != NULL ) tmp_cov_diag = cov_diag [0];

     /*  ALLOCATE AN OUTPUT BYTE IMAGE WITH NO MAP DATA */

     *out_image = createimage (
          in_image->col_size,
          in_image->row_size,
          (unsigned long) VFF_TYP_1_BYTE,
          (unsigned long) 1,
          (unsigned long) 1,
          "from cluster reduction program",
          (unsigned long) 0,
          (unsigned long) 0,
          (unsigned long) VFF_MS_NONE,
          (unsigned long) VFF_MAPTYP_NONE,
          (unsigned long) VFF_LOC_IMPLICIT,
          (unsigned long) 0
     );

     /*  CHECK FOR VALID MERGING ROUTINE  */

     if ( alg == NULL ) {
          fprintf (stderr,"No algorithm name given\n");
          return (FALSE);
     }

     if ( alg [0] == '1' ) {

          merge_alg_1 ( (int)in_image->ispare2, &num_entries, final_num,
               means, count, cov_upper, cov_diag, mapping );

     } else if ( alg [0] == '2' ) {

          if ( (means == NULL) || (count == NULL) ) {
               fprintf (stderr,"mean, count are not available\n");
               return(FALSE);
          }

          merge_alg_2 ( (int)in_image->ispare2, &num_entries, final_num,
               means, count, cov_upper, cov_diag, mapping );

     } else if ( alg [0] == '3' ) {

          if ( (means == NULL) || (count == NULL) || (cov_upper == NULL) ) {
               fprintf (stderr,"mean, count, and cov are not available\n");
               return(FALSE);
          }

          merge_alg_3 ( (int)in_image->ispare2, &num_entries, final_num,
               means, count, cov_upper, mapping );

     } else if ( alg [0] == '4' ) {
          fprintf (stderr,"Algorithm 4 not currently implemented.\n");
          return (FALSE);
     } else {
          fprintf (stderr,"Invalid algorithm name given\n");
          return (FALSE);
     }

     /*  USE THE MAPPING TO RE-MAP THE IMAGE DATA  */

     ptr_out = (unsigned char *) (*out_image)->imagedata;

     if ( mapping_empty != NULL ) {
     switch ( in_image->data_storage_type ) {
          case VFF_TYP_1_BYTE:
               ptr_1_BYTE = (unsigned char *) in_image->imagedata;
               for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                    ptr_out [k] = (unsigned char) 
                         (*mapping) [mapping_empty[ptr_1_BYTE [k]]];
               }
               break;
          case VFF_TYP_2_BYTE:
               ptr_2_BYTE = (unsigned short *) in_image->imagedata;
               for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                    ptr_out [k] = (unsigned char) 
                         (*mapping) [mapping_empty[ptr_2_BYTE [k]]];
               }
               break;
          case VFF_TYP_4_BYTE:
               ptr_4_BYTE = (int *) in_image->imagedata;
               for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                    ptr_out [k] = (unsigned char) 
                         (*mapping) [mapping_empty[ptr_4_BYTE [k]]];
               }
               break;
          default:
               fprintf (stderr,"Image data not char, short, or int...\n");
               fprintf (stderr,"output image data will not be valid\n");
     }
     free ( mapping_empty );
     mapping_empty = NULL;
     } else {
     switch ( in_image->data_storage_type ) {
          case VFF_TYP_1_BYTE:
               ptr_1_BYTE = (unsigned char *) in_image->imagedata;
               for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                    ptr_out [k] = (unsigned char) (*mapping) [ptr_1_BYTE [k]];
               }
               break;
          case VFF_TYP_2_BYTE:
               ptr_2_BYTE = (unsigned short *) in_image->imagedata;
               for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                    ptr_out [k] = (unsigned char) (*mapping) [ptr_2_BYTE [k]];
               }
               break;
          case VFF_TYP_4_BYTE:
               ptr_4_BYTE = (int *) in_image->imagedata;
               for (k=0;k<(in_image->col_size)*(in_image->row_size);k++){
                    ptr_out [k] = (unsigned char) (*mapping) [ptr_4_BYTE [k]];
               }
               break;
          default:
               fprintf (stderr,"Image data not char, short, or int...\n");
               fprintf (stderr,"output image data will not be valid\n");
     }
     }

     /*  SCALE THE MEAN VECTOR DATA TO BETWEEN 0 AND 255?  */

     if ( scale ) {

          if ( global_scale ) {

               min = means [0][0];
               max = means [0][0];

               for ( cluster = 0 ; cluster < num_entries ; cluster ++ ) {
               for ( c = 0 ; c < (int) in_image->ispare2 ; c ++ ) {
                    if ( means [cluster][c] < min ) min = means [cluster][c];
                    if ( means [cluster][c] > max ) max = means [cluster][c];
               }
               }

               if ( max == min ) {
                    fprintf (stderr,"Warning: min & max are equal....\n");
                    fprintf (stderr,"         no scaling done\n");
               } else {
                    scale = 255.0 / ( max - min ) ;
                    for ( cluster = 0 ; cluster < num_entries ; cluster ++ ) {
                    for ( c = 0 ; c < (int) in_image->ispare2 ; c ++ ) {
                         means [cluster][c] = (means [cluster][c]-min)*scale;
                    }
                    }
               }

          } else {

               for ( c = 0 ; c < (int) in_image->ispare2 ; c ++ ) {

                    min = means [0][c];
                    max = means [0][c];
                    for ( cluster = 0 ; cluster < num_entries ; cluster ++ ) {
                         if (means[cluster][c] < min) min = means [cluster][c];
                         if (means[cluster][c] > max) max = means [cluster][c];
                    }

                    if ( max == min ) {
                         fprintf (stderr,"Warning: min & max are equal....\n");
                         fprintf (stderr,"         no scaling done (%d)\n",c);
                    } else {
                         scale = 255.0 / ( max - min ) ;
                         for (cluster=0;cluster<num_entries;cluster++) {
                              means [cluster][c] = 
                                   (means [cluster][c]-min)*scale;
                         }
                     }

               }

          }

     }

     /*  STORE THE DATA BACK INTO MAP AREA FOR IMAGE  */

     (*out_image)->map_row_size = in_image->ispare2;
     if ( class ) ((*out_image)->map_row_size) ++ ;  
     if ( count ) ((*out_image)->map_row_size) ++ ;   
     if ( cov_diag ) ((*out_image)->map_row_size) += in_image->ispare2 ; 
     if ( cov_upper ) ((*out_image)->map_row_size) += 
          ((in_image->ispare2 * (in_image->ispare2 + 1))/2) ; 

     (*out_image)->map_col_size = num_entries;
     (*out_image)->map_storage_type = VFF_MAPTYP_FLOAT;
     (*out_image)->map_scheme = VFF_MS_SHARED;
     (*out_image)->ispare2 = in_image->ispare2;

     if ( ! putstats (&((*out_image)->maps),&((*out_image)->ispare1),
          (int)in_image->ispare2, num_entries, means, count, class, 
          cov_upper, cov_diag) ) {

          fprintf(stderr,"putstats Failed\n");
          return (FALSE);

     }

     /*  FREE MEMORY  */

     if ( means ) {
          free ( tmp_means );
          free ( means );
     }

     if ( count ) {
          free ( count );
     }

     if ( class ) {
          free ( class );
     }

     if ( cov_diag ) {
          free ( tmp_cov_diag );
          free ( cov_diag );
     }

     if ( cov_upper ) {
          for ( i = 0 ; i < num_entries ; i ++ ) {
               free ( cov_upper [i][0] );
               free ( cov_upper [i] );
          }
          free ( cov_upper );
     }

   return(TRUE);
}

/* -library_code_end */
