_delaunay_tree.c

数据类型和算法库LEDA 数据类型和算法库LEDA

/*******************************************************************************
+
+  LEDA  3.0
+
+
+  _delaunay_tree.c
+
+
+  Copyright (c) 1992  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 6600 Saarbruecken, FRG     
+  All rights reserved.
+ 
*******************************************************************************/


/*****************************************************************************
+
+	Delaunay tree    by Olivier Devillers
+
+		Fully dynamic construction of the Delaunay triangulation with
+		good randomized complexity
+
+	Copyright (c) 1990
+	INRIA, 2004 route des Lucioles, 06 565 Valbonne Cedex, France
+
+
******************************************************************************/


#include <LEDA/impl/delaunay_tree.h>


#define Fini		0
#define Infini1 	1		/* 1 point a l'infini */
#define Infini2		2		/* 2 points a l'infini */
#define Infini3		3		/* 3 points a l'infini, c'est la racine */
#define Impossible	99		/* Le voisin du precedent */
#define U 0				/* indices pour les trois sommets */
#define V 1
#define W 2


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

			Declarations de type

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


struct LISTENOEUD{
	NOEUD *item;
	LISTENOEUD *next;
};



struct REINSERER {
	DT_item	x;				/* a point devant etre reinserer */
	NOEUD   *detruit1,			/* triangles a detruire si ils existent */
		*detruit2;			/* NULL sinon xpx1 est direct et xpx2 indirect*/
	listenoeud
		decroches;			/* triangles decroches */
	REINSERER
		*finf,*fsup;		/* pour l'arbre binaire */
};



struct STAR{
	STAR
			*next,*previous;/*arete suivante et precedente sur le polygone*/
	NOEUD
			*triangle;		/* triangle adjacent */
	int		n,p,arete;		/*point suivant, precedent, arete dans le triangle*/
};




struct NOEUD{

	int		type;	/* le triangle est Fini ou Infini[123] 
				   dans le cas d'un triangle infini, les sommets
				   infinis sont la direction a` l'infini */
	
	int		degenere;  /* booleen, 3 pts alignes ? */
	int		direct;	   /* booleen, le triangle est-il direct ? */
        int             stacked;
	int		visite;	   /* triangle marque visite si ce flag
			              est le flag courant.
				      pour la suppression le triangle est
	                              marque si il appartient ou a
	                              appartenu a A 
	                            */
	
	int		detruire;   /* le triangle est a detruire si ce flag
				       est le flag courant */
	coord 	cx,cy;		    /* centre du cercle 
	                                      (uniqu^t pour tracer voronoi)*/
	coord	a2,b2,a3,b3,c1,z0,z1,z2; /* pour les tests de conflit */
	
	DT_item	meurtrier; /* NULL si le triangle n'est pas mort */
	DT_item	s[3];	   /* sommets du triangle */
				   /* U est le createur */
				   /* si Infini1 alors W est a l'infini */
				   /* si Infini2 alors V et W sont a l'infini */
				   /* si Infini3 alors U,V et W sont a l'infini*/
	NOEUD
			*pere,		/* pere */
			*bpere,		/* beau-pere */
			*fils[3];	/* fils */
	listenoeud
			bfils;		/* liste des beaux-fils */
	NOEUD
			*voisin[3],	/* 3 voisins (courant ou a la mort) */
			*creation[3],	/* 3 voisins a la creation */
			*special[3];	/* 3 voisins speciaux! */
	star	star[3];		/* pointeur reciproque vers Star */


};

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

			Allocation memoire

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

#define SIZEALLOC 100


char* delaunay_tree::alloc_block(int size)
{
   char* p =  (char*)malloc( size + sizeof(double) ); 

   if (!p) error_handler(1,"out of memory");

   *((char**)p) = used_blocks;
   used_blocks = p;

   return (used_blocks+sizeof(double)) ; 
 }


void delaunay_tree::free_blocks()
{
  char* p;
  while (used_blocks)
  { p = used_blocks;
    used_blocks = *((char**)used_blocks);
    free(p);
   }
}

void delaunay_tree::poubelle_point(DT_item p)
{ p->cree = (noeud) Poubelle_point;
  Poubelle_point = p;
}

DT_item delaunay_tree::vrai_alloc_point()
{
 if ( ! item_nombre ) {
		item_next = (DT_item) alloc_block( (item_nombre = SIZEALLOC) * sizeof(DT_POINT) ) ; 
	}
	item_nombre--;
	return item_next++;
}

DT_item delaunay_tree::alloc_point()
{
	if ( Poubelle_point ) { DT_item p = Poubelle_point;
				Poubelle_point = (DT_item) Poubelle_point->cree;
				return p;  }
	return vrai_alloc_point();
}


void delaunay_tree::poubelle_listenoeud(listenoeud p)
{
  p->next =  Poubelle_listenoeud;
  Poubelle_listenoeud = p;
}

listenoeud delaunay_tree::vrai_alloc_listenoeud()
{
 if ( ! listenoeud_nombre ) {
		listenoeud_next = (listenoeud) alloc_block((listenoeud_nombre = SIZEALLOC) * sizeof(LISTENOEUD) ) ; 
	}
	listenoeud_nombre--;
	return listenoeud_next++;
}

listenoeud delaunay_tree::alloc_listenoeud()
{
	if ( Poubelle_listenoeud ) 
                              { listenoeud p = Poubelle_listenoeud;
			        Poubelle_listenoeud = Poubelle_listenoeud->next;
				return p;  }
	return vrai_alloc_listenoeud();
}


void delaunay_tree::poubelle_noeud(noeud p)
/* on ne se resert que des noeud dont le champs detruire n'est pas flag */
{ 
  p->fils[0] =  Poubelle_noeud;
  Poubelle_noeud = p;
  if (fl != p->detruire ) { fl = p->detruire; Libre_noeud = Poubelle_noeud ;}
}

noeud delaunay_tree::vrai_alloc_noeud()
{
 if ( ! noeud_nombre ) {
		noeud_next = (noeud) alloc_block((noeud_nombre = SIZEALLOC) * sizeof(NOEUD) ) ; 
	}
	noeud_nombre--;
	return noeud_next++;
}

noeud delaunay_tree::alloc_noeud()
/* on ne se resert que des noeud dont le champs detruire n'est pas flag */
{	
	if (Libre_noeud && Libre_noeud->fils[0] )
			{ noeud p = Libre_noeud->fils[0];
			  Libre_noeud->fils[0] = Libre_noeud->fils[0]->fils[0];
			  return p;  }
	return vrai_alloc_noeud();
}


void delaunay_tree::poubelle_star(star p)
{
  p->next =  Poubelle_star;
  Poubelle_star = p;
}

star delaunay_tree::vrai_alloc_star()
{
 if ( ! star_nombre ) {
		star_next = (star) alloc_block((unsigned int) (star_nombre = SIZEALLOC) * sizeof(STAR) ) ; 
	}
	star_nombre--;
	return star_next++;
}


star delaunay_tree::alloc_star()
{
	if ( Poubelle_star ) { star p = Poubelle_star;
			       Poubelle_star = Poubelle_star->next;
			       return p;  }
	return vrai_alloc_star();
}


void delaunay_tree::poubelle_reinserer( reinserer p)
{
  p->finf =  Poubelle_reinserer;
  Poubelle_reinserer = p;
}
 
reinserer delaunay_tree::vrai_alloc_reinserer()
{
 if ( ! reinserer_nombre ) {
		reinserer_next = (reinserer) alloc_block((reinserer_nombre = SIZEALLOC) * sizeof(REINSERER) ) ; 
	}
	reinserer_nombre--;
	return reinserer_next++;
}


reinserer delaunay_tree::alloc_reinserer()
{
	if ( Poubelle_reinserer ) {reinserer p = Poubelle_reinserer;
			   	   Poubelle_reinserer = Poubelle_reinserer->finf;
				   return p;  }

	return vrai_alloc_reinserer();
}

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

			Structure Reinsert

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

reinserer delaunay_tree::locate( reinserer n, DT_item x)
{
	reinserer nn;

	if ( ! Reinsert ) {
		nn = alloc_reinserer();
		nn->x = x;
		nn->decroches = NULL;
		nn->detruit1 = NULL;
		nn->detruit2 = NULL;
		nn->finf = NULL;
		nn->fsup = NULL;
		Reinsert = nn;
	 return nn;
	}
	if ( n->x == x) return n;
	if ( x->n < n->x->n ) {
		if (n->finf) return locate(n->finf,x);
		nn = alloc_reinserer();
		nn->x = x;
		nn->decroches = NULL;
		nn->detruit1 = NULL;
		nn->detruit2 = NULL;
		nn->finf = NULL;
		nn->fsup = NULL;
		n->finf = nn;
	} else {
		if (n->fsup) return locate(n->fsup,x);
		nn = alloc_reinserer();
		nn->x = x;
		nn->decroches = NULL;
		nn->detruit1 = NULL;
		nn->detruit2 = NULL;
		nn->finf = NULL;
		nn->fsup = NULL;
		n->fsup = nn;
	}
	return nn;
}

void delaunay_tree::clear_Reinsert( reinserer n)
{	listenoeud l,ll;
	if ( ! Reinsert ) return;
	if (n->finf) clear_Reinsert(n->finf);
	if (n->fsup) clear_Reinsert(n->fsup);
	for ( l = n->decroches; l; l = ll ) {ll = l->next; poubelle_listenoeud(l);}
	poubelle_reinserer(n);
	if ( n == Reinsert ) Reinsert = NULL;
}

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

			Structure Star

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

void delaunay_tree::clear_Star()
{
	star s;

	if ( ! Star ) return;
	for ( s = Star->next ; s != Star ;)
				{s = s->next; poubelle_star(s->previous);}
	poubelle_star(Star);
	Star = NULL;
}

void delaunay_tree::init_Star( noeud n, int i, int j, int k)
/* l'arete i,j de n est sur le bord de star */
{
	star s;

	if (! n) {
		Star->next = first_star;
		first_star->previous = Star;
		first_star = NULL;
		return;
	}
	s = alloc_star();
	s->triangle = n;
	s->p = i;
	s->n = j;
	s->arete = k;
	n->star[k] = s ;
	s->previous = Star;
	if (! Star ) first_star = s;
	else		 Star->next = s;
	Star = s;
}

star delaunay_tree::loc_Star( DT_item xx)
/* trouve l'arete telle que xx soit le point precedent */
{
	star s;

	s = Star;
	while ( s->triangle->s[s->p] != xx ) s = s->next ;
	return s;
}

void delaunay_tree::split_Star( star n1,star n2)
/* raccroche les n1 a n2 en viarant la partie intermediare */
{
	star s,ss;

	if (n1==n2) {
		s = alloc_star();
		s->previous = n1;
		s->next = n1->next;
		n1->next = s;
		s->next->previous = s;
		return;
	}
	for ( s = n1->next; s != n2 ; ){ss=s->next; poubelle_star(s); s=ss; }
	n1->next = n2;
	Star = n2->previous = n1;
}

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

			Procedures

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

void delaunay_tree::beaufils( noeud bf, noeud bp)
/* ajoute bf a la liste des beaufils de bp */
{
	listenoeud nov;

	nov = alloc_listenoeud();
	nov->next = bp->bfils;
	nov->item = bf;
	bp->bfils = nov;
}

void delaunay_tree::plusbeaufils( noeud bf, noeud bp)
/* enleve bf a la liste des beaufils de bp */
{
	listenoeud n,nn;

	if ( bp->bfils->item == bf ) {
		nn = bp->bfils;
		bp->bfils = bp->bfils->next ;
		poubelle_listenoeud(nn);
		return;
	}
	for ( n = bp->bfils; n->next->item != bf; n = n->next );
	nn = n->next ;
	n->next = n->next->next ;
	poubelle_listenoeud(nn);
}

void delaunay_tree::demiplan( noeud t)
/* calcul de la normale */
{
	t->type = Infini1;
	t->degenere = 0;
	if ( t->direct ) {
			t->cx= t->s[U]->y - t->s[V]->y ;
			t->cy= t->s[V]->x - t->s[U]->x ;
	} else {
			t->cx= t->s[V]->y - t->s[U]->y ;
			t->cy= t->s[U]->x - t->s[V]->x ;
	}
}

void delaunay_tree::cercle( noeud t)
{
	register double a,b,c,d;
/* calcul pour les tests de conflits */
   	t->a2 = t->s[1]->x - t->s[0]->x ;
   	t->b2 = t->s[1]->y - t->s[0]->y ;
   	t->a3 = t->s[2]->x - t->s[0]->x ;
   	t->b3 = t->s[2]->y - t->s[0]->y ;
 
   	t->c1 = t->a2 * t->b3 - t->a3 * t->b2 ;
   	t->z1 = t->a2 * (t->s[1]->x + t->s[0]->x) + t->b2 * (t->s[1]->y+t->s[0]->y);
   	t->z2 = t->a3 * (t->s[2]->x + t->s[0]->x) + t->b3 * (t->s[2]->y+t->s[0]->y);
/* calcul du centre et du rayon du cercle  */
	a = t->s[U]->x * t->s[U]->x + t->s[U]->y * t->s[U]->y ;
	b = t->s[V]->x * t->s[V]->x + t->s[V]->y * t->s[V]->y - a ;
	c = t->s[W]->x * t->s[W]->x + t->s[W]->y * t->s[W]->y - a ;
	d = 2*( (t->s[V]->x - t->s[U]->x) * (t->s[W]->y - t->s[U]->y)
				- (t->s[W]->x - t->s[U]->x) * (t->s[V]->y - t->s[U]->y));
	if ( (t->degenere = (d==0) ) ) { 
						t->cx= t->pere->cx ;
						t->cy= t->pere->cy ;
					return;
				}
	t->cx=(b*(t->s[W]->y - t->s[U]->y)-c*(t->s[V]->y - t->s[U]->y))/d;
	t->cy=(c*(t->s[V]->x - t->s[U]->x)-b*(t->s[W]->x - t->s[U]->x))/d;
}

int delaunay_tree::confondu( DT_item a, DT_item b)
{
	return ((a->x == b->x)&&(a->y==b->y));
}

int delaunay_tree::direct( noeud n)
/* teste si un triangle est direct */
{
	switch(n->type){
	case Fini       : 
		return ((n->s[V]->x - n->s[U]->x)*(n->s[W]->y - n->s[U]->y)
				- (n->s[V]->y - n->s[U]->y)*(n->s[W]->x - n->s[U]->x) > 0);
	case Infini1	 :
		return ((n->s[V]->x - n->s[U]->x)*n->s[W]->y
				- (n->s[V]->y - n->s[U]->y)*n->s[W]->x > 0);
	case Infini2     :
		return (n->s[V]->x*n->s[W]->y - n->s[V]->y*n->s[W]->x > 0);
	case Infini3     : return 1;
	case Impossible  : return 1;
	}
	return 2;
}

 
int delaunay_tree::conflit(noeud n, DT_item p)
{
  /* teste si un point est en conflit avec une region */

    register coord a1,c2,c3,z0;

    switch(n->type){
    case Fini       :
        a1 = p->x - n->s[0]->x ;
        c3 = p->y - n->s[0]->y ;

        z0 = a1 * (p->x + n->s[0]->x) + c3 * (p->y + n->s[0]->y) ;

        c2 = n->a3 * c3 - a1 * n->b3 ;
        c3 = a1 * n->b2 - n->a2 * c3 ;

            /********** Arithmetique double precision ************/