 /*
  * Khoros: $Id: lvlabel.c,v 1.2 1992/03/20 23:06:24 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvlabel.c,v 1.2 1992/03/20 23:06:24 dkhoros Exp $";
#endif

 /*
  * $Log: lvlabel.c,v $
 * Revision 1.2  1992/03/20  23:06:24  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * 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: lvlabel.c
 >>>>
 >>>>      Program Name: vlabel
 >>>>
 >>>> Date Last Updated: Tue Mar  5 22:20:58 1991 
 >>>>
 >>>>          Routines: lvlabel - the library call for vlabel
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define abs(x) ((x > 0) ? x : -x )

#define MAX_LOOP_NUMBER  50                /* maxi nbr of iteration for*/
                                           /* automatic research ......*/
#define UNDEFINED_REGION 0                 /* UNDEFINED  REGION NUMBER */
#define FIRST_REGION     1                 /* FIRST      REGION NUMBER */
#define MAX_REGION       1073741825        /* Maximum number of region */
#define CUR_REGION       -1     /* current region in a merging process */
#define INDEX(i,j) (i*nc + j)
#define IND(i,j,k) (i*nc + j + k*nc*nr)

/* global variables  that are available for any routine in this file   */

static int nc;    /* number of columns in theses images  ..............*/
static int nr;    /* number of rows in theses images ..................*/
 
static int i_start,i_end;   /* first and last row in image workspace ..*/
static int j_start,j_end;   /* first and last column in image workspace*/

static int band;      /* number of bands in the multi band image, ..   */
                      /*     can be 1 ..............................   */
static int num_clus;  /* number of cluster centers .................   */
static int num_reg;   /* current number of region for labelling.*/

static float thresh;  /* max amplitude in the distance between 2 neighb*/
static int   surf_min;/* minimum surface of a region ...............   */

static float *image;  /* pointer on multiband image data ...........   */
static int *label;    /* pointer on label image data ..........   */


static int   *clus_num;/* pointer on cluster numbers image data .....  */
static float *clus_cen;/* pointer on cluster centers images data ....  */



/* the following images gives the distance between two neighbours.     */
/* The number of different images should be 8 for 8 connectivity, and  */
/* 4 for the 4 connectivity, but in fact, 4 is enough for 8 connectiv. */
/* and 2 is enough for 4 connectivity.                                 */
/* Exemple: dist_neigh[index(i,j)][k] is the distance between the      */   
/* pixel(i,j) and its neighbour number k, knowing that the neighbours  */
/* number are :     0    for   i    ,   j+1                            */
/*                  1    for   i-1  ,   j                              */
/*                  2    for   i-1  ,   j+1                            */
/*                  3    for   i-1  ,   j-1                            */
/* this order is not conventional, but it help to compute in the same  */
/* program the 4 and 8 connectivity.
/* Also the 4 others neighbours 4,5,6,7 are not mentionned here because*/
/* of the symetricity of the distance metrics (dist(A,B) = dist(B,A)   */
/* Therefore if we consider that the neighbour number 4 is the symetric*/
/* of the neighbour number 0, the distance between (i,j) and its       */
/* neighbour 4 is equal to distance[0][neighbour(4),pixel(i,j)]........*/
/* the neighbour number 5 is symetric to number 1                      */
/* the neighbour number 6 is symetric to number 2                      */
/* the neighbour number 7 is symetric to number 3                      */
static float **distance;/* pointer on distance images (number of "band"*/
                        /* depend on connectivity: 2 for 4 connectivity*/
                        /*                         4 for 8 connectivity*/

static int metric;      /* type of distance use to link pixels ....... */
static float split;     /* split and merge factor .................... */
static int neigh_num2; /* Number of neighb. to check for distance comp.*/

static int neigh_i[4]= {0,-1,-1,-1};        
static int neigh_j[4]= {1,0,1,-1};           

/* Structure that link the pixel that can be label in the same   */
/* region ..........................................................   */
static struct LIST_CONNECT {
       int i;                /* i coordinate of the candidate pixel    */
       int j;                /* j coordinate of the candidate pixel    */
       struct LIST_CONNECT *next;  /* ptr on the next candidate ...    */
       }  *list_start;


/* the structure below describe all the region that have been labelled */
/* by the label routine .............................................. */
/* This can easily be completed for new purposes or for different      */
/* applications. ......  The structure is a structure of linked liste  */
/* Each element of the list described the region itself, and is linked */
/* via a pointer to the next region .................................. */
static struct REGION {
     int surface;               /* surface = number of pixel ......... */
     struct REGION *next;
     } *region_start, *reg;

/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvlabel - library call for vlabel
*
* Purpose:
*    
*    Performs  a  Labeling  on  Multiband  or   Cluster   xvimage
*    Structures
*    
*    
* Input:
* Output:
*
* Written By: Pascal ADAM
****************************************************************/


/* -library_def */
int lvlabel(ima,
            cc,
            cn,
            typ,
            dist,
            connec,
            bord,
            lab,
            surf,
            opt_auto,
            reg_num,
            split,
            merge,
            printdev)


struct xvimage *ima;  /* initial image if present                      */
struct xvimage *cc;   /* cluster center image if present               */
struct xvimage *cn;   /* cluster number image "   "                    */
struct xvimage *lab;  /* labelled image                                */

int typ;              /* case of labelling.............................*/
                      /*      if 0:  label from a multi/mono band only.*/
                      /*      if 1:  label from cluster images only....*/
                      /*      if 2:  label form cluster & mono/multi...*/
float surf;           /* nc * nr * surf_min = surface minimum under    */
                      /* which a labelled region is rejected during    */
                      /* the first pass of the algorithm...............*/
int dist;             /* 1: use euclidean distance to link pixel.......*/
                      /* 2: use city blocks distance...................*/
int connec;           /* type of connectivity used:                    */
                      /*        value = 1 for the ...... 4 connectivity*/
                      /*         ....   2 .............. 8 connectivity*/
int bord;             /* define the width of the border around the imag*/
                      /* that you do not process.......................*/
int opt_auto;         /* boolean, if true converge toward number of re-*/
                      /* gion expected by the user.....................*/
int reg_num;          /* number of labelled regions wanted.............*/
float split;          /* split and merge coefficient used to determine */
                      /* the threshold to label pixels in a connex reg.*/
int merge;            /* boolean, indicates if the small region should */
                      /* be merged together or with bigger ones, or if */
                      /* they should only be classified as undefined...*/
FILE *printdev;       /* Pointer on ASCII FIle for labelling info .....*/
                      /* storage.......................................*/

/* -library_def_end */

/* -library_code */
{
float f_distance();
void free_region();
void free_list();
struct LIST_CONNECT *new_list();
struct REGION *new_region();
int test_i();
int test_j();
void f_label();
void f_merge();
void converge();
void free();
char *malloc(),             /* function definition                     */
     *program ="lvabel";    /* Contains the library name               */

float sup_lim;        /* Superieur limit for the threshold label value */
int   sup_num;        /* number of region for the superior limit value */
float inf_lim;        /* Inferior  limit for the threshold label value */
int   inf_num;        /* number of region for the inferior limit value */
float bes_lim;        /* best value for the threshold label .......... */
int   bes_num;        /* best number of region for best threshold .... */
float cur_lim;        /* current threshold value ..................... */
int   cur_num;        /* number of region for the current thresh.value */
float las_lim;        /* last    threshold value ..................... */
int   las_num;        /* number of region for the last    thresh.value */

float step;           /* interval to find the proper threshold value   */
float distance_tmp;   /* temporary value                               */
int index_tmp;        /* temporary value                               */
int tmp_ave_reject;   /* temporary value                               */
int tmp_ave_kept;     /* temporary value                               */
int i,
    j,
    k,
    l,
    loop_num;/* scanning variables                            */


float dicho;          /* dichotomy variable............................*/
int fin;              /* test end of loop do while.....................*/

int *reg_resume;      /* array used to eliminate too small regions     */

int   clus_neig;      /* cluster number of  neighbor pixel of indice l */
int   clus_pix;       /* cluster number of current pixel               */

float *val_neigh;     /* Multiband value of the neighbours of the pixel*/
float *val_pixel;     /* Multiband value of the pixel.                 */ 


float *dist_cluster;  /* inter_cluster distances ..................... */
                      /* array [num_cluster][num_cluster] ............ */
                      /* dist_cluster[k][l] is the distance between    */
                      /* the cluster k and l.......................... */


float max,ave,min,var;/* statistic on distances .......................*/
float pix_num;        /* total number of pixel computed ...............*/


metric      = dist;
neigh_num2  = connec * 2;   /* due to the symetricity of the metrics   */
                            /* for the 4 connectivity we only need to  */
                            /* compute the distance for two neighbours */
                            /* to know the distance between any pixel  */



if ((typ == 0) || (typ == 2)) 
{
   image = (float *) ima->imagedata;
   band  = ima->num_data_bands;
   nc    = ima->row_size;
   nr    = ima->col_size;
}
if ((typ == 1))  
{
   clus_cen = (float *) cc->imagedata; 
   num_clus = cc->row_size;
   band     = cc->num_data_bands;
   dist_cluster=(float *)malloc((int)(num_clus*num_clus)*sizeof(float *));
   if (dist_cluster   == NULL) 
   {
      (void) fprintf(stderr,
      "%s:  No space for dist_cluster   _ malloc failed! \n", program);
      exit(1);
   }

}

if ((typ == 1) || (typ == 2)) 
{
   clus_num = (int *) cn->imagedata;
   nr       = cn->col_size;
   nc       = cn->row_size;
}

/* pointer on labelled image output   */
label = (int *) lab->imagedata; 


i_start = j_start = bord;
i_end   = nr-bord;
j_end   = nc-bord;

/* computes the minimum surface under which a region is rejected ........ */
surf_min    = (int) (surf * (float)(nc * nr) / 100.00);





/*....................................................................... */
/* computes the distance images dist_0, dist_2.. that gives the  distance */
/* between neighbours, respecting the chosen metric and connectivity .... */
/*....................................................................... */
/*....................................................................... */
/*....................................................................... */
/*....................................................................... */
/*....................................................................... */

/* dynamic allocation for dist_0 and dist_2 */
distance = (float **) malloc (neigh_num2 * sizeof(float *));
if (distance   == NULL)
{
   (void) fprintf(stderr,
   "lvabel:  No space for distance   _ malloc failed! \n");
   exit(1);
}

for (i = 0 ; i < neigh_num2 ; i++) 
{
   distance[i] = (float *) malloc ( nc * nr * sizeof(float));

   if (distance[i] == NULL) 
   {
      (void) fprintf(stderr,
      "lvabel:  No space for distance[%d] _ malloc failed! \n",i);
      exit(1);
   }
}


switch (typ) {
   /* case for which the labelling use only a mono or multi band image    */ 
   case 0:
      val_pixel = (float *) malloc (band * sizeof (float));
      if (val_pixel == NULL)
      {
         (void) fprintf(stderr,
         "lvlabel: No space for val_pixel _ malloc failed! \n");
         exit(1);
      }
      val_neigh = (float *) malloc (band * sizeof (float));
      if (val_neigh == NULL) 
      {
         (void) fprintf(stderr,
         "lvlabel: No space for val_neigh _ malloc failed! \n");
         exit(1);
      }


      for (i = i_start+1 ; i < i_end-1 ; i++) 
      {
         for (j = j_start+1 ; j < j_end-1 ; j++) 
         {
            for (l = 0; l<band ; l++) val_pixel[l] = image[IND(i,j,l)] ;  

            for (k = 0 ; k < neigh_num2 ; k ++) 
            { 
               for (l = 0; l<band ; l++) 
               {
                 val_neigh[l]=image[IND((i+neigh_i[k]),(j+neigh_j[k]),l)];
               }

               /* computes the distance between pixel and 2 or 4 neighbours */ 
               distance[k][INDEX(i,j)]=f_distance(val_pixel,val_neigh,band);
            }
         }
      }
      break;
   
   /* case for which the labelling use only a cluster number and cluster  */ 
   /* center image .....................................................  */
   case 1: 
      val_pixel = (float *) malloc (band * sizeof (float));
      if (val_pixel == NULL) 
      {
         (void) fprintf(stderr,
         "lvlabel: No space for val_pixel _ malloc failed! \n");
         exit(1);
      }
      val_neigh = (float *) malloc (band * sizeof (float));
      if (val_neigh == NULL) 
      {
         (void) fprintf(stderr,
         "lvlabel: No space for val_neigh _ malloc failed! \n");
         exit(1);
      }

      /* initialization of dist_cluster array ........................... */
      for (i = 0; i < num_clus; i++) 
      {
          for (j = i; j < num_clus; j++) 
          {
             for (l = 0; l < band ; l++) val_pixel[l] = clus_cen[(i+num_clus*l)];
             for (l = 0; l < band ; l++) val_neigh[l] = clus_cen[(j+num_clus*l)];

             dist_cluster[j*num_clus+i]=f_distance(val_pixel,val_neigh,band);
             dist_cluster[i*num_clus+j]= dist_cluster[j*num_clus+i];
          }
      }              

      for (i = i_start+1 ; i < i_end-1 ; i++) 
      {
         for (j = j_start+1 ; j < j_end-1 ; j++) 
         {
            clus_pix = clus_num[INDEX(i,j)];

            for (k = 0 ; k < neigh_num2 ; k ++) 
            {
               clus_neig = clus_num[INDEX((i+neigh_i[k]),(j+neigh_j[k]))];

               if (clus_pix == clus_neig)  distance[k][INDEX(i,j)] =0.0;
               else 
               {
                  /* computes distance between pixel and 2 or 4neighbours */
                  distance[k][INDEX(i,j)]=
                  dist_cluster[(clus_neig*num_clus+clus_pix)];
               }
            }
         }
      }
      break;

   /* case for which the labelling use a cluster number image and a    */ 
   /* multiband image ................................................ */ 
   case 2:
      val_pixel = (float *) malloc (band * sizeof (float));
      if (val_pixel == NULL) 
      {
         (void) fprintf(stderr,
         "lvlabel: No space for val_pixel _ malloc failed! \n");
         exit(1);
      }

      val_neigh = (float *) malloc (band * sizeof (float));
      if (val_neigh == NULL) 
      {
         (void) fprintf(stderr,
         "lvlabel: No space for val_neigh _ malloc failed! \n");
         exit(1);
      }

      for (i = i_start+1 ; i < i_end-1 ; i++) 
      {
         for (j = j_start+1 ; j < j_end-1 ; j++) 
         {
            clus_pix = clus_num[INDEX(i,j)];

            for (l = 0; l < band ; l++) val_pixel[l] = image[IND(i,j,l)] ;  

            for (k = 0 ; k < neigh_num2 ; k ++) 
            {
               clus_neig = clus_num[INDEX((i+neigh_i[k]),(j+neigh_j[k]))];

               if (clus_pix ==  clus_neig) distance[k][INDEX(i,j)] =0.0;
               else  
               {
                  for (l = 0; l<band ; l++) 
                  {
                     val_neigh[l]=image[IND((i+neigh_i[k]),(j+neigh_j[k]),l)];
                  }

                  /* computes the distance between pixel and 2 or 4 neighbours */
                  distance[k][INDEX(i,j)]=f_distance(val_pixel,val_neigh,band);
               }           
            }
         }
      }
      break;
}



/* At this point for the case  0, 1, 2  the information giving the distance */
/* between one pixel to any of its neighbour, whatever the connectivity is, */
/* is stored in the multiband image " d i s t a n c e".  the number of band */
/* varies from 2 to 4 depending on the connectivity.                        */
/* The next step will be to analyze the maximum variation or distance       */
/* between to neighbour to fix the threshold use by the labelling routine   */
/* to link pixels together in the same region number. This threshold will   */
/* be dependant also from the split and merge factor.                       */
/* 
/* The heuristic that will decide whether or not 2 pixels will be linked    */
/* is given below and is based on a threshold called   " m a x "            */
/* But this heuristic could be more "intelligent" and based on statistics   */
/* computed from the distance image ........................................*/
/* .........................................................................*/



/* search for maximum distance in distance image and computes average and   */
/* variance ................................................................*/

min = distance[0][INDEX((i_start+2),(j_start+2))]; max = 0.0; ave = 0.0; 
pix_num = (float)((i_end-i_start-4)*(j_end-j_start-4)); 

for (i = i_start+2 ; i < i_end-2 ; i++) 
{
   for (j = j_start+2 ; j < j_end-2 ; j++) 
   {
      index_tmp = INDEX(i,j);
      for (k = 0 ; k < neigh_num2 ; k ++) 
      {
         distance_tmp = distance[k][index_tmp];
         if (max < distance_tmp)   max = distance_tmp;
         if (min > distance_tmp)   min = distance_tmp;

         ave += distance_tmp;
      }
   }
}

/* computes average and variance...........................................*/
ave /= pix_num; 


/* print info on distance image in printdev ...............................*/
(void) fprintf (printdev, "vlabel Statistics\n");
(void) fprintf (printdev, "================\n\n");
(void) fprintf (printdev,
                "\tMinimum distance between pixels:    %f\n\n",min);
(void) fprintf (printdev,
                "\tMaximum distance between pixels:    %f\n\n",max);
(void) fprintf (printdev,
                "\tAverage distance between pixels:    %f\n\n",ave);

/* if automatic option selected, fixe parameters for iteration process ....*/
if (opt_auto == TRUE) 
{
                    
   /* fixes the initial step to increase or decrease the threshold during  */
   /* the iterative process ...............................................*/
   step = (min - max)/ (float) (MAX_LOOP_NUMBER);

   /* fixe superior and inferior limit for threshold used in labelling ....*/
   /* process. During the iteration research of a segmentation that gives..*/
   /* the expected number  or a number close from it, the threshold will...*/
   /* be scanning the interval [sup_lim, inf_lim]. It is not usefull to ...*/
   /* search outside of these limits because the segmentation would be to  */
   /* instable or unsucessfull.............................................*/
   /* sup_num and inf_num are respectively the number of region that would */
   /* be labelled when using a threshold of sup_lim, respectively inf_lim..*/

   sup_lim = min - step ; inf_lim = max + step ;  /* step is negatif ......*/
   sup_num = 0; inf_num = 0;


   /* cur_lim, cur_num  current threshold & current region nber associated */
   /* las_lim, las_num  last threshold used & region number associated.... */
   /* bes_lim, bes_num  best threshold & region number associated, gives   */
   /* the closest region number from the expected region number........... */
   /* Intialisation of these values .......................................*/

   cur_lim = las_lim = bes_lim = 0.0; cur_num = las_num = bes_num = 0;

   /* fixes maximum number of iterations for the best threshold research...*/   
   loop_num = MAX_LOOP_NUMBER;

   /* fin == end in french... boolean indicating if iteration should stop  */
   fin = FALSE;

   /* print info on automatic labelling in printdev .......................*/
   /* Attention for fprintf we have exchange sup & inf purposely because   */
   /* sup happened to be smaller than inf..................................*/
   (void) fprintf (printdev, "\n\n\nvlabel ~~ Automatic Labelling\n\n");
   (void) fprintf (printdev,
                   "\tSplit&Merge Factor Superior Limit ......:    %f\n",
                   inf_lim / (max * 0.5)); 
   (void) fprintf (printdev,
                   "\tSplit&Merge Factor Inferior limit ......:    %f\n\n",
                   sup_lim / (max * 0.5));
   
}

/* manual execution, the user has choosen a split and merge factor, only   */
/* one labelling process will be executed on the image ....................*/
else 
{
   /* init threshold by using max and the split_and_merge_factor ..........*/
   thresh = split * max * 0.5;

   loop_num = 1; fin = TRUE;


   /* print info on distance image in printdev ............................*/
   (void) fprintf (printdev, "\n\n\nvlabel ~~ Manual Labelling\n");
}

 

/* GLOBAL LOOP FOR RESEARCH OF A THRESHOLD THAT GIVES A REGION NUMBER AS   */
/* CLOSE AS POSSIBLE FROM THE EXPECTED ONE.................................*/
do
{
   /* if automatic option selected ........................................*/
   if (opt_auto == TRUE)
   {
      /* if number of loop reached 0 (30 have been excuted) ...............*/
      if (loop_num == 0) 
      {
         /* process the last iteration by labelling again with the best    */
         /* threshold value and indicates the END_OF ITERATION by setting  */
         /* the boolean fin to TRUE .......................................*/

         thresh = bes_lim; fin = TRUE;
      } 
      else 
      {
         /* beginning, first we have to determine what are the region nber */
         /* when using sup_lim and inf_lim.... These will be the first two */
         /* threshold used ............................................... */
         if (sup_num == 0)  thresh = cur_lim = sup_lim; 
         else 
         {
            if (inf_num == 0) thresh = cur_lim = inf_lim;  
            else thresh = cur_lim;
         }
      }
   }

   /* between two iterations the label image has to be reinitialise to 0 ..*/ 
   for (k=0;k < nc*nr; k++ ) label[k] = UNDEFINED_REGION;

   /* initialisation before labelling process .............................*/
   num_reg= FIRST_REGION; region_start=new_region(); reg=region_start;

   /* Labelling loop ......................................................*/
   for (i = i_start+2 ; i < i_end-2 ; i++) 
   { 
      for (j = j_start+2 ; j < j_end-2 ; j++) 
      { 
         if (label[INDEX(i,j)] ==  UNDEFINED_REGION) 
         {
            f_label(i,j); 

            reg->next = new_region(); reg = reg->next;

            num_reg +=  1;

            /* test if number of region is not greater than the maximum    */
            /* authorized ................................................ */
            if ( num_reg ==  MAX_REGION) 
            {
               (void) fprintf(stderr,
                      "lvlabel: Cannot label more than %d region\n\n",
                      MAX_REGION);
               return(1);
            }
         }
      }
   }

   /* print info on Threshold..... in printdev ............................*/
   if ((opt_auto == TRUE) && (fin == TRUE)) 
   {
      (void) fprintf (printdev,"\n\n\nvlabel ~~ Automatic Labelling\n");
      (void) fprintf (printdev,"    FINAL RESULT for BEST THRESHOLD\n");
   }

   (void) fprintf (printdev,
   "\n\n\tloop %d ... Split&Merge Factor = %f ...\n",
                   loop_num,
                   thresh / (max * 0.5));
   (void) fprintf (printdev,
   "\n\n\t~~~~ total number of region after labelling: %d\n",num_reg-1);


   /* Scanning of region information to mark the small regions and valid  */
   /* the big enough ones.................................................*/
   reg_resume = (int *)malloc(num_reg * sizeof(int));
   if (reg_resume == NULL) 
   {
      (void)fprintf(stderr,
      "lvlabel: No space for reg_resume _ malloc failed! \n");
      exit(1);
   }
   reg_resume[0] = 0;
   reg = region_start;
   j = 1; k = 0; l = 0;
   tmp_ave_reject = tmp_ave_kept = 0;
   for (i=1; i< num_reg;i++) 
   { 
      if (reg->surface > surf_min) 
      {
         reg_resume[i] = j++;

         /* keep track of average size of rejected small regions ........*/
         tmp_ave_kept += reg->surface; l++;
      }
      else 
      {
         reg_resume[i] =UNDEFINED_REGION;
    
         /* keep track of average size of rejected small regions ........*/
         tmp_ave_reject += reg->surface; k++;
      }

      reg = reg->next;
   }

   /* num_reg holds the final number of labelled regions ................*/
   num_reg = j;

   if (k != 0)
   {
      /* print info on region number  in printdev ............................*/
      (void) fprintf (printdev,
      "\t~~~~ Number of Rejected Region : %d (average size in pixel: %d)\n",
      k,tmp_ave_reject / k);
   }
   if (l != 0)
   {
      (void) fprintf (printdev,
      "\t~~~~ Number of Region kept     : %d (average size in pixel: %d)\n",
      l,tmp_ave_kept / l);
      (void) fprintf (printdev,
      "\t~~~~ Percentage of labelled pixels in the image : %d % \n\n",
      (tmp_ave_kept * 100) / (int)pix_num);
   }

   /* Gives the final label number to the regions kept during this first */
   /* part of the process. The small region are labelled UNDEFINED_REGION*/
   for (i = i_start ; i < i_end ; i++) { 
      for (j = j_start ; j < j_end ; j++) {
         label[INDEX(i,j)] = reg_resume[(label[INDEX(i,j)])]; 
      }
   }


   /* merging of small regions, if option merge selected.................*/
   if (merge == TRUE) 
   {
      for (i = i_start+2 ; i < i_end-2 ; i++) 
      {  
         for (j = j_start+2 ; j < j_end-2 ; j++) 
         {     
            if (label[INDEX(i,j)] == UNDEFINED_REGION) f_merge(i,j); 
         }
      } 

      tmp_ave_kept = 0;
      for (i = i_start+2 ; i < i_end-2 ; i++)
      {
         for (j = j_start+2 ; j < j_end-2 ; j++)
         {
            if (label[INDEX(i,j)] != UNDEFINED_REGION) tmp_ave_kept ++;
         }
      }


      /* print info on region number  in printdev ............................*/
      (void) fprintf (printdev, "\t~~~~ vlabel ~~ Merging of Small REGION\n");
      (void) fprintf (printdev,
      "\t~~~~ Final Number of Region : %d\n",num_reg-1 );

      (void) fprintf (printdev,
      "\t~~~~ Percentage of labelled pixels in the image : %d % \n",
      (tmp_ave_kept * 100) / (int)pix_num);



   } 

   /* free memory allocate for regions ..................................*/
   (void)free(reg_resume); (void)free_region();

   /* if thresh equals one of the limit initialise the limit numbers *****/
   if (opt_auto == TRUE) {
      if (thresh == sup_lim )  sup_num = cur_num = (int)num_reg;
      else 
      {
         if (thresh == inf_lim) inf_num = cur_num = (int)num_reg; 
         else cur_num = (int)num_reg;
      }

      /* if the number of region is the correct one it stops ****************/
      if ((int)num_reg == reg_num) fin = TRUE;
 
      /* other case .. call converge with both limit and decide the new val */
      else 
      {
         converge(&sup_lim,&sup_num,
                  &inf_lim,&inf_num,
                  &cur_lim,&cur_num,
                  &bes_lim,&bes_num,
                  &las_lim,&las_num,
                  reg_num,&step);
         if (abs(step) < 0.1) loop_num = 1;
      }

      /* decreases the number of loop by one ********************************/
      loop_num--;
   }
}  while (fin == FALSE);

return(1);
}
   
/***********************************************************************/
/* Function f_distance():                                              */
/*                     this function computes the distance between two */
/*                     vectors of dimension N.                         */
/*                     the metric use is either, euclidean, chessboard */
/*                     or cityblock depending on the value of the      */
/*                     argument dist of lvlabel                        */
/*                     N equal band if typ = 0 or 2                    */
/*                     or num_clus  if typ = 1                         */
/***********************************************************************/

/*
 *  Need to define ABS.   (John & Mark - Sat Jul 21 12:49:45 MDT 1990)
 */
#undef ABS
#define ABS(x) ((x) < 0 ? -(x) : (x))

float f_distance(vect1,vect2,n)
float *vect1, *vect2;
int n;
{
int k;
float result;

result = ABS((vect1[0]-vect2[0]));
switch (metric) {

   case 1:             /* euclidean metric ........................... */

      for (k=0;k<n;k++) result += (vect1[k]-vect2[k])*(vect1[k]-vect2[k]);
      result = (float) sqrt((double)(result));
      break;

   case 2:             /* City blocks metric.......................... */

      for (k=1;k<n;k++) result += ABS((vect1[k] - vect2[k]));
      break;
   }

return (result);
}


/***********************************************************************/
/* Function free_list():                                               */
/*                     free the first element of a list  and           */
/*                     replace the first element by the second one.    */
/***********************************************************************/

void free_region()
{
struct REGION *region_ptr;
void free();

do
{
   if (region_start->next == NULL) 
   {
      region_ptr = region_start;
      (void)free(region_ptr);
      region_start = NULL;
   }
   else
   {
      region_ptr = region_start;
      region_start = region_start->next;
      (void)free(region_ptr);
   }
} while (region_start != NULL);

}


/***********************************************************************/
/* Function free_list():                                               */
/*                     free the first element of a list  and           */
/*                     replace the first element by the second one.    */
/***********************************************************************/

void free_list()
{
struct LIST_CONNECT *list;
void free();

   if (list_start->next == NULL) 
   {
      list = list_start; (void)free(list); list_start = NULL;
   }
   else 
   { 
      list = list_start; list_start = list_start->next;
      (void)free(list);
   } 
}

/***********************************************************************/
/* Function new_list():                                                */
/*                    allocate memory for a new list                   */
/*                                                                     */
/***********************************************************************/
struct LIST_CONNECT *new_list()
{
struct LIST_CONNECT *list;
char *malloc();

list = (struct LIST_CONNECT *)malloc (sizeof(struct LIST_CONNECT));
if (list == NULL) {
   (void)fprintf(stderr,
          "lvabel:  No space for new list _ malloc failed! \n");
   exit(1);
   }
list->next = NULL;

return(list);
}

/***********************************************************************/
/* Function new_region():                                              */
/*                    allocate memory for a new region                 */
/*                                                                     */
/***********************************************************************/
struct REGION *new_region()
{
struct REGION *region;
char *malloc();

   region = (struct REGION *)malloc (sizeof(struct REGION));
   if (region == NULL) 
   {
      (void) fprintf(stderr,
             "lvabel:  No space for new region _ malloc failed! \n");
      exit(1);
   }
   region->next = NULL;

   return(region);
}

/***********************************************************************/
/* Function test_i(val) :                                              */
/*                   test if the integer val is between i_start & i_end*/
/*                                                                     */
/***********************************************************************/
int test_i(val) 
int val;
{
  if (( val >= i_start+2) && (val < i_end-2)) return(TRUE);
  else return(FALSE);
}


/***********************************************************************/
/* Function test_j(val) :                                              */
/*                   test if the integer val is between j_start & j_end*/
/*                                                                     */
/***********************************************************************/
int test_j(val) 
int val;
{
  if (( val >= j_start+2) && (val < j_end-2)) return(TRUE);
  else return(FALSE);
}



/***********************************************************************/
/*  Function label:                                                    */
/*               labelling of uniform regions relatively to the metric */
/*               computed before.                                      */
/*                                                                     */
/*               this function could be improved by adding heuristics  */
/*               for merging or split region.                          */
/***********************************************************************/
void f_label(i,j)
int i,j;
{
void free_list();
int test_i();
int test_j();
struct LIST_CONNECT *new_list();
int k;
int stop; 
struct LIST_CONNECT *list;
float dist_ij;
int ind_ij, ind_neigh;

       reg->surface = 1; 
      
       list_start = new_list();
        
       list = list_start; 
       list->i = i; list->j = j;

       stop = FALSE;

       ind_ij = INDEX(i,j);
       label[ind_ij] = num_reg;

       do { 


          /* If pixel located further than one pixel to the border   */
          if ((i < i_end-3) && (i >= i_start +3) && 
              (j < j_end-3) && (j >= j_start +3)) 
          { 

             /* Check neighbors and link them if dist < thresh       */
             /* The Number of neighbors depends of connectivity type */
             for (k = 0; k < neigh_num2 ; k++) 
             {
                dist_ij = distance[k][ind_ij];
                ind_neigh = INDEX((i+neigh_i[k]),(j+neigh_j[k]));

                if (label[ind_neigh]== UNDEFINED_REGION) 
                {
                   if (dist_ij < thresh) 
                   {
                      label[ind_neigh] = num_reg;
                      reg->surface +=  1;
                      list->next = new_list();
                      list = list->next;
                      list->i = i+neigh_i[k];
                      list->j = j+neigh_j[k];
                   }
                }

                ind_neigh = INDEX((i-neigh_i[k]),(j-neigh_j[k]));

                if (label[ind_neigh] == UNDEFINED_REGION) 
                {
                   if (distance[k][ind_neigh] < thresh) 
                   {
                      label[ind_neigh] = num_reg;
                      reg->surface +=  1;
                      list->next = new_list();
                      list = list->next;
                      list->i = i-neigh_i[k];
                      list->j = j-neigh_j[k];
                   }
                }
             }
          }

          /* If pixel has a location in the image such that one of   */
          /* its neighbor could be on the border, we check every ... */
          /* new pixel location to make sure that we will not write. */
          /* outside of the border defined by the user ............. */
          else 
          {                    
             for (k = 0; k < neigh_num2 ; k++) 
             {
                dist_ij = distance[k][ind_ij];

                if ((test_i(i+neigh_i[k]) == TRUE) &&
                   (test_j(j+neigh_j[k]) == TRUE) ) 
                {
                   ind_neigh = INDEX((i+neigh_i[k]),(j+neigh_j[k]));

                   if (label[ind_neigh] ==  UNDEFINED_REGION) 
                   { 
                      if (dist_ij < thresh) 
                      {
                         label[ind_neigh] = num_reg;
                         reg->surface +=  1;
                         list->next = new_list();
                         list = list->next;
                         list->i = i+neigh_i[k];
                         list->j = j+neigh_j[k];
                      }
                   } 
                }

                if ((test_i(i-neigh_i[k]) == TRUE) &&
                    (test_j(j-neigh_j[k]) == TRUE) ) 
                {
                   ind_neigh = INDEX((i-neigh_i[k]),(j-neigh_j[k]));
                   if (label[ind_neigh] ==  UNDEFINED_REGION ) 

                   {
                      if (distance[k][ind_neigh] < thresh) 
                      {
                         label[ind_neigh] = num_reg;
                         reg->surface +=  1;
                         list->next = new_list();
                         list = list->next;
                         list->i = i-neigh_i[k];
                         list->j = j-neigh_j[k];
                      }
                   }
                }
             }
          }


          /* the pixel (i,j) which was pointed by list_start has */
          /* been studied entirely, Now we can get rid of it and */
          /* study the next one in the list .................... */
          free_list();

          /* if the list is empty, every possible pixel has been */
          /* checked, We can quit f_label and go back to lvlabel */
          /* to study a new candidate region.................... */
          if (list_start == NULL) stop = TRUE;

          /* if the list is not empty, we start the process again*/
          /* taking for new pixel the following one in the list..*/
          else {
             i = list_start->i;
             j = list_start->j;
             ind_ij = INDEX(i,j);
             }
  
          } while (stop == FALSE);
       } 

   
/***********************************************************************/
/*  Function merge                                                     */
/***********************************************************************/
void f_merge(i,j)
int i,j;
{
void free_list();
int test_i();
int test_j();
struct LIST_CONNECT *new_list();

int k;
int new_label; 
int stop;
int ind_neigh;
int surf_region;     
struct LIST_CONNECT *list, *list_cursor;
  
   surf_region  = 1; 
      
   list_start = new_list();
        
   list_start->i = i; list_start->j = j;
   list = list_cursor = list_start; 


   label[INDEX(i,j)] = CUR_REGION;   

   do 
   { 
      /* Try to merge UNDEFINED_REGION pixel together form the */
      /* pixel which coordinate are defined below ............ */

      i = list_cursor->i; j = list_cursor->j;

      for (k = 0; k < neigh_num2 ; k++) 
      {
         if ((test_i(i+neigh_i[k])==TRUE)&&(test_j(j+neigh_j[k])==TRUE)) 
         {
            ind_neigh = INDEX((i+neigh_i[k]),(j+neigh_j[k]));

            if (label[ind_neigh] == UNDEFINED_REGION) 
            { 
               label[ind_neigh] = CUR_REGION;   
               surf_region +=  1;
               list->next = new_list();
               list = list->next;
               list->i = i+neigh_i[k];
               list->j = j+neigh_j[k];
            } 
         }

         if ((test_i(i-neigh_i[k])==TRUE)&&(test_j(j-neigh_j[k])==TRUE)) 
         {
            ind_neigh = INDEX((i-neigh_i[k]),(j-neigh_j[k]));
            if (label[ind_neigh] == UNDEFINED_REGION) 
            {
               label[ind_neigh] = CUR_REGION;   
               surf_region +=  1;
               list->next = new_list();
               list = list->next;
               list->i = i-neigh_i[k];
               list->j = j-neigh_j[k];
            }
         }
      }
      list_cursor =list_cursor->next;

   } while(list_cursor != NULL);
      

   /* Test if the new region (set of UNDEFINED PIXELS) is big...*/
   /* enough to be kept as a uniform region ....................*/
   /* If yes, a new region number is given to this set of pixel.*/
   /* If not, this new region will be merge to one of it connec-*/
   /* ted region ...............................................*/
   if (surf_region > surf_min) 
   {
      list = list_start;
      do 
      {
         label[INDEX((list->i),(list->j))] = num_reg;
         list = list->next;

      } while (list != NULL);

      num_reg += 1;
   }
  
   else 
   {
      list_cursor = list_start;
      new_label = CUR_REGION;   
      stop = FALSE;        
 
      do
      {
         i = list_cursor->i; j = list_cursor->j;

         for (k = 0; ((k < neigh_num2) && (stop == FALSE)) ; k++) 
         {
            if ((test_i(i+neigh_i[k])==TRUE)&&(test_j(j+neigh_j[k])==TRUE)) 
            {
               if (label[INDEX((i+neigh_i[k]),(j+neigh_j[k]))] != CUR_REGION) 
               {
                  new_label = label[INDEX((i+neigh_i[k]),(j+neigh_j[k]))];
                  stop = TRUE; 
               }
            }
            if ((test_i(i-neigh_i[k])==TRUE)&&
                (test_j(j-neigh_j[k])==TRUE)&&
                (stop == FALSE)) 
            {
               if (label[INDEX((i-neigh_i[k]),(j-neigh_j[k]))] != CUR_REGION) 
               {
                  new_label = label[INDEX((i-neigh_i[k]),(j-neigh_j[k]))];
                  stop = TRUE;
               }
            }
         }

         list_cursor = list_cursor->next;

      }  while ((list_cursor != NULL) || ( stop == FALSE));

      if (new_label == CUR_REGION) 
      {
         (void)fprintf(stderr,"lvlabel: internal bug in f_merge function\n");
      }
      else
      {
         list_cursor = list_start;

         do 
         {
            label[INDEX((list_cursor->i),(list_cursor->j))] = new_label;
            list_cursor = list_cursor->next;
         } while (list_cursor != NULL);
      }
   }
           

   /* free the entire list of linked pixels ...................... */
   while (list_start != NULL) free_list(); 
             
} 
      
/************************************************************************/
/* converge ............................................................*/
/************************************************************************/
void converge(sup_lim,sup_num,
              inf_lim,inf_num,
              cur_lim,cur_num,
              bes_lim,bes_num,
              las_lim,las_num,
              number,step)

int   *sup_num,*inf_num,*cur_num,*bes_num,*las_num,number;
float *sup_lim,*inf_lim,*cur_lim,*bes_lim,*las_lim,*step;
{


int  sup_dif, inf_dif, cur_dif, bes_dif;
int fin;


   sup_dif = *sup_num - number;
   inf_dif = *inf_num - number;
   cur_dif = *cur_num - number;
   bes_dif = *bes_num - number;
 
   if (*cur_lim != *sup_lim) 
   {
      /* if cur_lim = inf_lim, fixes the new thresh and computes best */
      if (*cur_lim == *inf_lim)
      {
         fin = FALSE;
         do 
         {
            if ((*cur_lim + *step) > *sup_lim) 
            {
               *cur_lim = *cur_lim + *step;
               fin = TRUE;
            }
            else *step = *step / 2.0;
         } while (fin == FALSE);  

         if (abs(inf_dif) < abs(sup_dif) )
         {
            *bes_num = *inf_num;
            *bes_lim = *inf_lim;
         }
         else
         {
            *bes_num = *sup_num;
            *bes_lim = *sup_lim;
         }
         *las_num = *cur_num;
         *las_lim = *cur_lim;
      }
      else
      /* if normal case,  determines the new thresh****************/
      {
         /* first modifies the new best values ********************/
         if ( (abs(cur_dif) < abs(bes_dif)) ||
              ( (abs(cur_dif) == abs(bes_dif)) && (*bes_lim < *cur_lim)) )
         {
            *bes_num = *cur_num;
            *bes_lim = *cur_lim;
         }

         /* increases cur_lim by the discrete value ****************/
         if (*las_num == *cur_num) 
         {
            *las_lim = *cur_lim;
            *las_num = *cur_num;
            *step = *step * 1.3;
         }
         else 
         {
            if ( ((*las_num > number) && (*cur_num < number)) ||
                 ((*las_num < number) && (*cur_num >  number)) )  
            {
               *step = (*cur_lim - *las_lim) / 2.0; 
               *cur_lim = *las_lim;
               *cur_num = *las_num;
            }
            else
            {
              *las_lim = *cur_lim;
              *las_num = *cur_num;
            }
         }    
 
         fin = FALSE;
         do 
         {
            if ((*cur_lim + *step) > *sup_lim) 
            {
               *cur_lim = *cur_lim + *step;
               fin = TRUE;
            }
            else *step = *step / 2.0;

         } while (fin == FALSE);
      }
   }
 
   /* case :   cur_lim = sup_lim 1st step of the algorithm **************/
   else 
   {
      *las_lim = *cur_lim;
      *las_num = *cur_num;
   }
}
/* -library_code_end */
