 /*
  * Khoros: $Id: lvclose.c,v 1.3 1992/03/20 23:04:44 dkhoros Exp $
  */

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

 /*
  * $Log: lvclose.c,v $
 * Revision 1.3  1992/03/20  23:04:44  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, 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 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvclose.c
 >>>>
 >>>>      Program Name: vclose
 >>>>
 >>>> Date Last Updated: Fri Feb 28 19:46:26 1992 
 >>>>
 >>>>          Routines: lvclose - the library call for vclose
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#include "vipl/lvclose.h"

static int  adjacents();
static int  rechercher_candidats();
static int  depiler_pile_cdats();
static int  extremite();
static int  max_grdt();
static int  sommet();
static void initialisation();
static void elimination();
static void prolongation();
static void empiler_pile_cdats();

/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvclose - library call for vclose
*
* Purpose:
*    
*    Close boundaries in images
*    
*    

* Input:
*    
*         edge - pointer on  the  data  of  a  xvimage  VFF_TYP_1_BYTE
*         "Boundaries" image
*         gradient - pointer on the data of a  xvimage  VFF_TYP_1_BYTE
*         "Gradient" image
*         nc - (int) number of columns of both images
*         nr - (int) number of rows of both images
*         length - (int) minimum gap in pixel to close two edges.
*         thresh - (int) Threshold used to close the  edges  from  the
*         gradient data
*    
*         columns as img1 and img2.
*    
*    

* Output:
*    
*         edge - holds the resulting image with closed boundaries
*    
*    

*
* Written By: Pascal ADAM
*    
*    Fixed indexing bug in function max_grdt() so that  no  nega-
*    tive indexes, same problem rechercher_candidats().
*    
*    John Rasure Fri Feb 28 19:46:13 MST 1992
*    
*    

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


/* -library_def */
int lvclose(edge,gradient,nc,nr,length,thresh)
unsigned char  *edge; 
unsigned char  *gradient;
int            length, 
               nc ,
               nr,
               thresh;
/* -library_def_end */

/* -library_code */
{
void creation_regions();
int            dim,
               i;
unsigned char  *out;

        dim = nc * nr;

        i_dim = nr;
        j_dim = nc;

        coef_seuil_grdt = thresh;

        lg_min_ch_front = length;

        out    = (unsigned char *) malloc(dim * sizeof(unsigned char));
        im_ctr = (short *) malloc(dim * sizeof(short));
        
        if (out == NULL || im_ctr == NULL)
        {
           (void)fprintf(stderr,"lvclose: No Space For Intermediate Images _ malloc failed! \n");
           return(0);
        }

        im_grdt = gradient;

        for (i=0; i<dim; i++)
           if (edge[i] > 0) im_ctr[i] = -1;
           else             im_ctr[i] =  0;

       
        /* Create a blank one pixel width frame on the border of the edge image im_ctr */
        for (i = 0; i < nc; i++) im_ctr[i]      = im_ctr[(nr - 1) * nc + i] = -1;
        for (i = 0; i < nr; i++) im_ctr[i * nc] = im_ctr[i * nc + nr - 1]   = -1;

        /* Create "region" */
        creation_regions();


        /* create result image and overwrite edge */
        for (i = 0; i < dim ; i++)
           if (im_ctr[i] == -1) edge[i] = 255;
           else edge[i] = 0;
 
       return(1);
}



/*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*\
|* Creation des regions d'une image par fermeture des chaines ouvertes  *|
|* entrees: fichier image des points contours                           *|
|*          fichier image des gradients                                 *|
|*          taille de l'image                                           *|
|*          coefficient du seuil du gradient adaptatif                  *|
|*          longueur minimale d'une frontiere                           *|
|* sortie: fichier image des regions                                    *|
\*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<*<>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*>*/
void creation_regions()
    {
     register int i,point;
          /* elimination(); */
/* Initialisation du tableau donnant la position des 8 voisins d'un point de l'image, relativement a celui-ci */
     initialisation();

/* Balayage de l'image et fermeture des contours                        */
     for (i=2*j_dim;i<(i_dim-2)*j_dim;i+=j_dim)
         for (point=i+2;point<i+j_dim-2;point++)
             if (POINT_CONTOUR(point))
                fermeture(point);

/* Balayage de l'image et elimination des chaines restees ouvertes      */
     for (i=j_dim;i<(i_dim-1)*j_dim;i+=j_dim)
         for (point=i+1;point<i+j_dim-1;point++)
             if (POINT_CONTOUR(point))
                elimination(point);

    }

static void initialisation()
      {
       voisinage[0]= -1;
       voisinage[1]= j_dim-1;
       voisinage[2]= j_dim;
       voisinage[3]= j_dim+1;
       voisinage[4]= 1;
       voisinage[5]= -j_dim+1;
       voisinage[6]= -j_dim;
       voisinage[7]= -j_dim-1;
      }

fermeture(point)
int point;
   {
    int dir,
        candidat;

    sommet_pile_cdats= NULL;
    dir=extremite(point,"prolongation");

/* Elimination des points isoles                                        */
    if (dir== -1)
       COLORIER(point,FOND);

/* Fermeture d'un contour par prolongation d'une extremite de chaine de points contours */
    else if (dir>0)
           {
            seuil_grdt=coef_seuil_grdt*GRADIENT(point);
            if ((candidat=rechercher_candidats(point,dir))!= NULL)
              {
               prolongation(point,candidat);
               return(couleur== CONTOUR);
              }
           }
    return(SUCCES);
   }

static int extremite(pt,traitement)
int pt;
char *traitement;
    {
     int k,
         entourage[8],
         nb_voisins=0;

     for (k=0;k<8;k++)
        {
         entourage[nb_voisins]=pt+voisinage[k];
         if (POINT_CONTOUR(entourage[nb_voisins]))
            nb_voisins++;
        }
     if (nb_voisins==0)
        return(-1);
     if (nb_voisins==1)
        return(cde_dir[(entourage[0]-pt+j_dim+1)/j_dim][(entourage[0]-pt+j_dim+1)%j_dim]);
     if ((nb_voisins==2)&&adjacents(entourage[0],entourage[1],4))
       {
        k=(((abs(entourage[0]-pt)==j_dim)||(abs(entourage[0]-pt)==1))?0:1);
        return(cde_dir[(entourage[k]-pt+j_dim+1)/j_dim][(entourage[k]-pt+j_dim+1)%j_dim]);
       }
     if (!strcmp(traitement,"elimination"))
        switch (nb_voisins)
              {
               case 2: if (adjacents(entourage[0],entourage[1],8))
                          return(-1);
                       else return(0);
               case 3:
               case 4:
                      {
                       int disjoint=0,
                           k=1;

                       while ((disjoint<2)&&(k<nb_voisins))
                            {
                             if (!adjacents(entourage[k-1],entourage[k],8))
                                disjoint++;
                             k++;
                            }
                       if (((disjoint==1)&&(!adjacents(entourage[nb_voisins-1],entourage[0],8)))||(disjoint>1))
                          return(0);
                      }
               case 5: if (!((POINT_CONTOUR(pt+voisinage[0])||((POINT_CONTOUR(pt+voisinage[1])||POINT_CONTOUR(pt+voisinage[2]))&&(POINT_CONTOUR(pt+voisinage[7])||POINT_CONTOUR(pt+voisinage[6]))))&&(POINT_CONTOUR(pt+voisinage[4])||((POINT_CONTOUR(pt+voisinage[5])||POINT_CONTOUR(pt+voisinage[6]))&&(POINT_CONTOUR(pt+voisinage[3])||POINT_CONTOUR(pt+voisinage[2]))))))
                          return(-1);
               default: return(0);
              }
     return(0);
    }

static int rechercher_candidats(pt,dir)
int pt,dir;
   {
    int dir1=dir%8,
        dir2=(dir-1+8)%8,   /* fixed by Rasure */
        dir3=(dir+1)%8,
        candidat1=pt+voisinage[dir1],
        candidat2=pt+voisinage[dir2],
        candidat3=pt+voisinage[dir3];

    if (GRADIENT(candidat1)>seuil_grdt)
       if ((GRADIENT(candidat1)+max_grdt(candidat1,dir1))<(GRADIENT(candidat2)+max_grdt(candidat2,dir2)))
          if ((GRADIENT(candidat2)+max_grdt(candidat2,dir2))<(GRADIENT(candidat3)+max_grdt(candidat3,dir3)))
            {
             empiler_pile_cdats(pt,candidat1);
             empiler_pile_cdats(pt,candidat2);
             return(candidat3);
            }
          else
              {
               if ((GRADIENT(candidat1)+max_grdt(candidat1,dir1))<(GRADIENT(candidat3)+max_grdt(candidat3,dir3)))
                 {
                  empiler_pile_cdats(pt,candidat1);
                  empiler_pile_cdats(pt,candidat3);
                 }
               else
                   {
                    if (GRADIENT(candidat3)>seuil_grdt)
                       empiler_pile_cdats(pt,candidat3);
                    empiler_pile_cdats(pt,candidat1);
                   }
               return(candidat2);
              }
       else if (GRADIENT(candidat2)>seuil_grdt)
               if ((GRADIENT(candidat2)+max_grdt(candidat2,dir2))<(GRADIENT(candidat3)+max_grdt(candidat3,dir3)))
                 {
                  empiler_pile_cdats(pt,candidat2);
                  if ((GRADIENT(candidat1)+max_grdt(candidat1,dir1))<(GRADIENT(candidat3)+max_grdt(candidat3,dir3)))
                    {
                     empiler_pile_cdats(pt,candidat1);
                     return(candidat3);
                    }
                  else
                      {
                       empiler_pile_cdats(pt,candidat3);
                       return(candidat1);
                      }
                 }
               else
                   {
                    if (GRADIENT(candidat3)>seuil_grdt)
                       empiler_pile_cdats(pt,candidat3);
                    empiler_pile_cdats(pt,candidat2);
                    return(candidat1);
                   }
            else if ((GRADIENT(candidat1)+max_grdt(candidat1,dir1))<(GRADIENT(candidat3)+max_grdt(candidat3,dir3)))
                   {
                    empiler_pile_cdats(pt,candidat1);
                    return(candidat3);
                   }
                 else
                     {
                      if (GRADIENT(candidat3)>seuil_grdt)
                         empiler_pile_cdats(pt,candidat3);
                      return(candidat1);
                     }
    else if (GRADIENT(candidat2)>seuil_grdt)
            if ((GRADIENT(candidat2)+max_grdt(candidat2,dir2))<(GRADIENT(candidat3)+max_grdt(candidat3,dir3)))
              {
               empiler_pile_cdats(pt,candidat2);
               return(candidat3);
              }
            else
                {
                 if (GRADIENT(candidat3)>seuil_grdt)
                    empiler_pile_cdats(pt,candidat3);
                 return(candidat2);
                }
         else if (GRADIENT(candidat3)>seuil_grdt)
                 return(candidat3);
              else return(NULL);
   }

static void prolongation(point,candidat)
int point,candidat;
   {
    int dir,
        nouv_candidat;
    void initialiser_pile_cdats();
    int            boucle();

    COLORIER(candidat,CDAT_CTR);
    if ((dir=extremite(candidat,"prolongation"))>0)
      {
       if (((nouv_candidat=rechercher_candidats(candidat,dir))!= NULL)&&(!boucle(nouv_candidat,dir)))
          prolongation(candidat,nouv_candidat);
       else couleur= FOND;
       COLORIER(candidat,couleur);
       if ((couleur== FOND)&& sommet(point))
          prolongation(point,depiler_pile_cdats());
      }
    else
        {
         initialiser_pile_cdats();
         COLORIER(candidat,CONTOUR);
         couleur= CONTOUR;
        }
   }

static int max_grdt(pt,dir)
int pt,dir;
   {
    int grdt_pt1= GRADIENT(pt+voisinage[dir%8]),
        /* grdt_pt2= GRADIENT(pt+voisinage[dir-1]),  John Rasure */
        grdt_pt2= GRADIENT(pt+voisinage[(dir-1+8)%8]),
        grdt_pt3= GRADIENT(pt+voisinage[(dir+1)%8]);

    return((grdt_pt1<grdt_pt2)?((grdt_pt2<grdt_pt3)?grdt_pt3:grdt_pt2):grdt_pt1);
   }

static void empiler_pile_cdats(point,candidat)
int point,
    candidat;
   {
    t_pile_cdats *noeud;
    
    if ((noeud= MALLOC(t_pile_cdats))== NULL)
    {
       fprintf(stderr,"max_grdt: can't allocate memory\n\n");
       exit(1);
    }
 
    noeud->point=point;
    noeud->candidat=candidat;
    noeud->suivant=sommet_pile_cdats;
    sommet_pile_cdats=noeud;
   }

static sommet(point)
int point;
   {
    if (sommet_pile_cdats== NULL)
       return(VIDE);
    return(point==sommet_pile_cdats->point);
   }

static depiler_pile_cdats()
      {
       t_pile_cdats *noeud;
       int candidat;

       noeud=sommet_pile_cdats;
       candidat=noeud->candidat;
       sommet_pile_cdats=noeud->suivant;
       free((char *)noeud);
       return(candidat);
      }

static void initialiser_pile_cdats()
      {
       t_pile_cdats *noeud;

       while ((noeud=sommet_pile_cdats)!= NULL)
            {
             sommet_pile_cdats=noeud->suivant;
             free((char *)noeud);
            }
      }

static void elimination(point)
int point;
   {
    int dir,
        lg_ch_front=0;

    sommet_pile_cdats= NULL;
    while ((dir=extremite(point,"elimination"))>0)
         {
          if (lg_ch_front>=lg_min_ch_front)
            {
             for (;lg_ch_front>0;lg_ch_front--)
                 COLORIER(depiler_pile_cdats(),CONTOUR);
             return;
            }
          empiler_pile_cdats(0,point);
          lg_ch_front++;
          COLORIER(point,FOND);
          point+=voisinage[(dir+4)%8];
         }
    if ((dir== -1)&&(lg_ch_front<lg_min_ch_front))
       COLORIER(point,FOND);
    if (lg_ch_front<=lg_min_ch_front)
       initialiser_pile_cdats();
   }

static int adjacents(pt1,pt2,type)
int pt1,
    pt2,
    type;
{
int res;
        if (type == 4)
        {
           res = (abs(pt1-pt2) == 1) || (abs(pt1-pt2) == j_dim);
        }
        else 
        {
           res = ( (abs(pt1-pt2) == 1) ||
                 ( (abs(pt1-pt2) >= j_dim-1) && (abs(pt1-pt2)<=j_dim+1)));
        }
        return(res);
}

static int boucle(nouv_candidat,dir)
int  dir,
     nouv_candidat;
{

     return (    (im_ctr[nouv_candidat+voisinage[(dir+6)%8]] == CDAT_CTR) ||
                 (im_ctr[nouv_candidat+voisinage[(dir+7)%8]] == CDAT_CTR) ||
                 (im_ctr[nouv_candidat+voisinage[dir%8]]     == CDAT_CTR) || 
                 (im_ctr[nouv_candidat+voisinage[(dir+1)%8]] == CDAT_CTR) ||
                 (im_ctr[nouv_candidat+voisinage[(dir+2)%8]]== CDAT_CTR)   ) ;
}

/* -library_code_end */
