_delaunay_tree.c

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

            /*      On utilise le fait que les flottants de l'accelerateur
                    flottant sont stockes avec une mantisse de 64 bits !!!
                    Pour faire du calcul entier juste, faite le en flottant !
            */

        a1=(double)z0*(double)n->c1 + (double)n->z1*(double)c2 +
                                    (double)n->z2*(double)c3;
        return (n->direct ? ( a1 <= 0.0)  : (a1 >= 0.0) );
    case Infini1    :
        if ( ( (p->x - n->s[U]->x)* (n->s[U]->y - n->s[V]->y)
                + (p->y - n->s[U]->y)* (n->s[V]->x - n->s[U]->x) ) == 0 )
        {               z0 = n->s[U]->x - p->x;
                        a1 = n->s[U]->y - p->y;
                        c2 = n->s[V]->x - p->x;
                        c3 = n->s[V]->y - p->y;
                        return ( z0*c2 + a1*c3 < 0);
        }else return (n->direct)
            ?   ( ( (p->x - n->s[U]->x)* (n->s[U]->y - n->s[V]->y)
                    + (p->y - n->s[U]->y)* (n->s[V]->x - n->s[U]->x) ) > 0 )
            :   ( ( (p->x - n->s[U]->x)* (n->s[V]->y - n->s[U]->y)
                    + (p->y - n->s[U]->y)* (n->s[U]->x - n->s[V]->x) ) > 0 ) ;
    case Infini2    :
                return ( (p->x - n->s[U]->x) * (n->s[V]->x + n->s[W]->x)
                    +(p->y - n->s[U]->y) * (n->s[V]->y + n->s[W]->y) > 0 );
    case Infini3    : return 1;
    case Impossible : return 0;
    }
    return 0;
}

//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 ************/
//	/* On utilise le fait que les flottants de l'accelerateur
//	   flottant sont stockes avec une mantisse de 64 bits !!!
//           Pour faire du calcul entier juste, faite le en flottant !
//	*/
// 
//	a1=(double)z0*(double)n->c1 + (double)n->z1*(double)c2 + (double)n->z2*(double)c3;
//		return (n->direct ? ( a1 <= 0.0)  : (a1 >= 0.0) );
//
//	case Infini1    : 
//		return (n->direct)
//			?	( ( (p->x - n->s[U]->x)* (n->s[U]->y - n->s[V]->y)
//					+ (p->y - n->s[U]->y)* (n->s[V]->x - n->s[U]->x) ) >= 0 )
//			:	( ( (p->x - n->s[U]->x)* (n->s[V]->y - n->s[U]->y)
//					+ (p->y - n->s[U]->y)* (n->s[U]->x - n->s[V]->x) ) >= 0 ) ;
//	case Infini2    :
//				return ( (p->x - n->s[U]->x) * (n->s[V]->x + n->s[W]->x)
//					+(p->y - n->s[U]->y) * (n->s[V]->y + n->s[W]->y) >= 0 );
//	case Infini3    : return 1;
//	case Impossible : return 0;
//	}
//	return 0;
//}



void delaunay_tree::clear()
{
  if (arbre != nouveau)
  { flag++;
    arbre->s[U]->visite = flag;
    arbre->s[V]->visite = flag;
    arbre->s[W]->visite = flag;
    clear_items(nouveau);
   }
  free_blocks();
  init();
}

delaunay_tree::~delaunay_tree() { clear(); }

delaunay_tree::delaunay_tree()  { init(); }

void delaunay_tree::init()
/* initialisation du Delaunay tree avec les trois points a l'infini */
{
        counter = 0;
        flag=1; 
        Star = NULL; 
        Reinsert = NULL; 
        Poubelle_noeud = Libre_noeud = NULL;

        used_blocks=NULL;
        Poubelle_point=NULL;
        Poubelle_listenoeud=NULL;
        Poubelle_noeud=NULL;
        Libre_noeud=NULL;
        Poubelle_star=NULL;
        Poubelle_reinserer=NULL;
        fl = 0;
        item_nombre = 0;
        item_next = 0;
        listenoeud_nombre = 0;
        listenoeud_next = 0;
        noeud_nombre = 0;
        noeud_next = 0;
        star_nombre = 0;
        star_next = 0;
        reinserer_nombre = 0;
        reinserer_next = 0;
        first_star=NULL;


	arbre = alloc_noeud();
	nouveau = arbre;

	arbre->type = Infini3;
	arbre->degenere = 0;
	arbre->visite = 0;
	arbre->stacked = 0;
	arbre->detruire = 0;
	arbre->meurtrier = NULL;
	arbre->fils[U] = arbre->fils[V] = arbre->fils[W] = NULL ;
	arbre->bfils = NULL;
	arbre->pere = NULL;
	arbre->bpere = NULL;
	arbre->voisin[U] =
	arbre->voisin[V] =
	arbre->voisin[W] =
	arbre->creation[U] =
	arbre->creation[V] =
	arbre->creation[W] = alloc_noeud();
	arbre->direct = 1;

	arbre->voisin[W]->visite = 0;
	arbre->voisin[W]->stacked = 0;
	arbre->voisin[W]->type = Impossible ;
	arbre->voisin[W]->degenere = 0;
	arbre->voisin[W]->meurtrier = NULL ;
	arbre->voisin[W]->fils[U] = arbre->voisin[W]->fils[V] =
							arbre->voisin[W]->fils[W] = NULL ;
	arbre->voisin[W]->bfils = NULL ;
	arbre->voisin[W]->pere = NULL ;
	arbre->voisin[W]->bpere = NULL ;
	arbre->voisin[W]->voisin[U] =
	arbre->voisin[W]->voisin[V] =
	arbre->voisin[W]->voisin[W] =
	arbre->voisin[W]->creation[U] =
	arbre->voisin[W]->creation[V] =
	arbre->voisin[W]->creation[W] = arbre ;
	arbre->voisin[W]->direct = 0;

	arbre->s[U] = arbre->voisin[W]->s[U] = alloc_point();
	arbre->s[V] = arbre->voisin[W]->s[V] = alloc_point();
	arbre->s[W] = arbre->voisin[W]->s[W] = alloc_point();
	arbre->s[U]->x =  2;
	arbre->s[U]->y =  0;
	arbre->s[V]->x = -1;
	arbre->s[V]->y =  2;
	arbre->s[W]->x = -1;
	arbre->s[W]->y = -2;
	arbre->s[U]->n = 0;
	arbre->s[V]->n = 0;
	arbre->s[W]->n = 0;

}

noeud delaunay_tree::Insert( noeud n, DT_item p)
/* recherche d'un conflit */
{
	listenoeud l;
	noeud nn;
	int i;

	if ( flag != n->visite ) {
	 n->visite = flag;
	 if ( conflit(n,p)){
		if ( n->meurtrier ) {
			for (i=U; i<=W; i++) if ( n->fils[i] )
				if ( nn = Insert( n->fils[i], p) ) return nn ;
		} else return n ;
		for ( l = n->bfils; l ; l = l->next )
				if ( nn = Insert( l->item, p) ) return nn ;
	 }
	}
	return NULL ;
}

noeud delaunay_tree::creenoeud( noeud pere, int i, DT_item p)
{
				nouveau = pere->fils[i] = alloc_noeud();
				beaufils(nouveau, pere->voisin[i] );
				nouveau->visite = flag;
				nouveau->stacked = flag;
				nouveau->detruire = 0;
				nouveau->meurtrier = NULL;
				nouveau->fils[U] = nouveau->fils[V] = nouveau->fils[W] = NULL ;
				nouveau->bfils = NULL;
				nouveau->s[U]=p;					/* createur */
				switch(i){
				case U: nouveau->direct = pere->direct;
						nouveau->s[V]= pere->s[V];
						nouveau->s[W]= pere->s[W]; break;
				case V: nouveau->direct = ! pere->direct;
						nouveau->s[V]= pere->s[U];
						nouveau->s[W]= pere->s[W]; break;
				case W: nouveau->direct = pere->direct;
						nouveau->s[V]= pere->s[U];
						nouveau->s[W]= pere->s[V]; break;
				}
				nouveau->pere = pere;
				nouveau->bpere =
				nouveau->voisin[U] =
				nouveau->creation[U] = pere->voisin[i];
										/* nouveau est-il fini ou infini */
				switch ( pere->type ) {
				case Fini       :
					cercle(nouveau);
					nouveau->type = Fini; break;
				case Infini1    :
					switch (i) {
					case U: 
					case V: demiplan(nouveau); nouveau->type = Infini1;
							nouveau->degenere = 0; break;
					case W: cercle(nouveau); nouveau->type = Fini; break;
					}break;
				case Infini2    :
					switch (i) {
					case U: nouveau->type = Infini2;
							nouveau->degenere = 0;break; 
					case V:
					case W: nouveau->type = Infini1; demiplan(nouveau);
							nouveau->degenere = 0;break; 
					}break;
				case Infini3    : nouveau->type = Infini2;
									nouveau->degenere = 0;break; 
				case Impossible : break;
				}
	return nouveau;
}

int delaunay_tree::creation( noeud detruit, DT_item p)
/* creation eventuelle des nouveaux noeuds et des liens et des voisinages*/
{
	DT_item first= NULL;
	noeud n1,father;
	int arete;		/* indice sur n1 de l'arete axe-createur */ 
	int a,b,c;		/* indice des trois sommets sur father */
	int a1,b1,c1;	/* indice des trois sommets sur le noeud precedent */
	int j;

	if ( ! detruit)
					return 0;
		for (j=U; j<= W - detruit->type; j++) if (confondu(p,detruit->s[j]))
									return 0;
		detruit->meurtrier = p;
		a = U; b = V; c = W;
		while ( detruit->voisin[c]->meurtrier
			|| conflit ( detruit->voisin[c],p) )
                {
		   detruit->voisin[c]->meurtrier = p;
		   a1 = a; b1 = b; c1 = c;
		   for (j=U; j<=W; j++)
		   if (detruit->voisin[c1]->s[j] == detruit->s[a1]) 
                       a=j;
		   else 
                      if (detruit->voisin[c1]->s[j] == detruit->s[b1]) 
                         c=j;
                      else
                         b=j;
		   detruit = detruit->voisin[c1] ;
		}
				p->cree = n1 = creenoeud(detruit,c,p);
				/* n1 remplace pere comme voisin de son beau-pere*/
				for (j=U; j<=W; j++)
					if (	(detruit->voisin[c]->s[j] != n1->s[V])
						&&	(detruit->voisin[c]->s[j] != n1->s[W])) break;
				detruit->voisin[c]->voisin[j] = n1 ;
				switch (c) {
					case U : first = detruit->s[a=(detruit->direct) ? W : V ];
									break;
					case V : first = detruit->s[a=(detruit->direct) ? U : W ];
									break;
					case W : first = detruit->s[a=(detruit->direct) ? V : U ];
									break;
				}
				b = c; c = U+V+W -b -a ;
				father = detruit;
				arete =  (n1->direct) ? V : W ;
	do {
				/* ac est l'arete par laquelle n1 est fils de father */
				/* on veut tourner autour de a en faisant bouger c vers b */
				/* i e on tourne dans le sens direct */
		while ( 	father->voisin[c]->meurtrier
				||	conflit ( father->voisin[c],p) 		){
			father->voisin[c]->meurtrier = p;
			a1 = a; b1 = b; c1 = c;
			for (j=U; j<=W; j++)
				if (father->voisin[c1]->s[j] == father->s[a1])			a=j ;
				else if (father->voisin[c1]->s[j] == father->s[b1])		c=j ;
				else													b=j;
			father = father->voisin[c1] ;
		}
		if (father->fils[c])
				n1->creation[arete] = n1->voisin[arete] = father->fils[c] ;
		else {  n1->creation[arete] = n1->voisin[arete] = creenoeud(father,c,p);
				/* le nouveau remplace pere comme voisin de son beau-pere*/
				for (j=U; j<=W; j++)
					if (	(father->voisin[c]->s[j] != father->fils[c]->s[V])
						&&	(father->voisin[c]->s[j] != father->fils[c]->s[W]))
																		break;
				father->voisin[c]->voisin[j] = father->fils[c] ;
		}
		arete =  (father->fils[c]->direct) ? V : W ;
		father->fils[c]->creation[V+W-arete] =
						father->fils[c]->voisin[V+W-arete] = n1 ;
		n1 = father->fils[c] ;
		j = a ;
		a = b ;
		b = c ;
		c = j ;
	} while ( father->s[a] != first );
	return 1;
}

DT_item delaunay_tree::localise( noeud detecte, DT_item p)
/* determination du plus proche voisin */
{
	DT_item first= NULL,trouve = NULL;
	coord dist, distance = 10000000.0;
	noeud father;
	int a,b,c;		/* indice des trois sommets sur father */
	int a1,b1,c1;	/* indice des trois sommets sur le noeud precedent */
	int j;

	if (!detecte)
					return NULL;

		for (j=U; j<= W - detecte->type; j++) if (confondu(p,detecte->s[j]))
									return detecte->s[j];
		detecte->visite = flag++ ;
		a = U; b = V; c = W;
		while ( 	( detecte->voisin[c]->visite == flag )
				||	conflit ( detecte->voisin[c],p) 		){
			detecte->voisin[c]->visite = flag;
			a1 = a; b1 = b; c1 = c;
			for (j=U; j<=W; j++)
				if (detecte->voisin[c1]->s[j] == detecte->s[a1])		a=j ;
				else if (detecte->voisin[c1]->s[j] == detecte->s[b1])	c=j ;
				else													b=j;
			detecte = detecte->voisin[c1] ;
		}
				switch (c) {
					case U : first = detecte->s[a=(detecte->direct) ? W : V ];
									break;
					case V : first = detecte->s[a=(detecte->direct) ? U : W ];
									break;
					case W : first = detecte->s[a=(detecte->direct) ? V : U ];
									break;
				}
				b = c; c = U+V+W -b -a ;
				father = detecte;
	do {
				/* on veut tourner autour de a en faisant bouger c vers b */
				/* i e on tourne dans le sens direct */
		while ( 	( father->voisin[c]->visite == flag )
				||	conflit ( father->voisin[c],p) 		){
			father->voisin[c]->visite = flag;
			a1 = a; b1 = b; c1 = c;
			for (j=U; j<=W; j++)
				if (father->voisin[c1]->s[j] == father->s[a1])			a=j ;
				else if (father->voisin[c1]->s[j] == father->s[b1])		c=j ;
				else													b=j;
			father = father->voisin[c1] ;
		}
		if (father->s[a]->n) {
			dist = (p->x - father->s[a]->x) * (p->x - father->s[a]->x) +
					(p->y - father->s[a]->y) * (p->y - father->s[a]->y)	;
			if (dist < distance) { distance = dist; trouve = father->s[a]; }
		}
		j = a ;
		a = b ;
		b = c ;
		c = j ;
	} while ( father->s[a] != first );
	return trouve;
}

void delaunay_tree::add_x1x2( reinserer r, noeud v, DT_item p)
/* ajoute v a comme triangle detruit de r */
{
	if ( v->direct )
					 if  (p==v->s[V]) r->detruit1 = v;
					 else			  r->detruit2 = v;
	else			
					 if  (p==v->s[W]) r->detruit1 = v;
					 else			  r->detruit2 = v;
}

void delaunay_tree::add_decroche( reinserer r, noeud v)
/* ajoute v a la liste des decroches de r */
{
	listenoeud nov;

	nov = alloc_listenoeud();
	nov->next = r->decroches;
	nov->item = v;
	r->decroches = nov;
}

void delaunay_tree::destroy( noeud u, DT_item p)
/* recherche recursive pour initialisation de Reinsert */
{
	int i;
	listenoeud l;
	reinserer rr;

	if ( ! u ) return;
	if (u->detruire == flag) return;
	for ( i=U ; i<=W; i++ ) if (u->s[i]==p) break;
	if ( i<=W ) {
		u->detruire = flag;
		if (u->meurtrier)
			for ( i=U ; i<=W; i++ ) if (u->fils[i]) destroy(u->fils[i], p);
		for ( l=u->bfils; l; l = l->next )      destroy(l->item   , p);
		if (u->s[U]==p) return;
		rr = locate(Reinsert,u->s[U]);
		add_x1x2(rr,u,p);
	} else {
		rr = locate(Reinsert,u->s[U]);
		add_decroche(rr,u);
		if (u->pere->detruire == flag)	u->pere		= NULL;
		else							u->bpere	= NULL;
	}
}

void delaunay_tree::recherche( DT_item p)
/* init des structures Star et Reinsert */
{
	DT_item first;
	noeud father;
	int arete;		/* indice sur n1 de l'arete axe-createur */ 
	int a,b,c;		/* indice des trois sommets sur father */
	int a1,b1,c1;	/* indice des trois sommets sur le noeud precedent */
	int j;

	clear_Reinsert(Reinsert); 
	clear_Star();
		arete = (  p->cree->direct  ) ? V : W ;
		first = p->cree->s[V + W - arete ];
		father = p->cree->pere;
		for (j=U; j<=W; j++)
			if (father->s[j] == p->cree->s[V + W - arete ])		a=j ;
			else if (father->s[j] == p->cree->s[arete])			c=j ;
			else								b=j;
	do {
				/* ac est l'arete par laquelle on est fils de father */
				/* on veut tourner autour de a en faisant bouger c vers b */
				/* et que on tourne autour de a dans le sens inverse */
		while ( ( father->voisin[c]->meurtrier == p )
				|| (father->voisin[c]->visite == flag) ) {
			a1 = a; b1 = b; c1 = c;
			for (j=U; j<=W; j++)
				if (father->voisin[c1]->s[j] == father->s[a1])			a=j ;
				else if (father->voisin[c1]->s[j] == father->s[b1])		c=j ;
				else													b=j;
			father = father->voisin[c1] ;
			father->visite = flag ;
			father->meurtrier = NULL;
		}
		father->visite = flag ;
		father->meurtrier = NULL;
		a1 = a; b1 = b; c1 = c;
		for (j=U; j<=W; j++)
				if (father->voisin[c1]->s[j] == father->s[a1])			a=j ;
				else if (father->voisin[c1]->s[j] == father->s[b1])		c=j ;
				else													b=j;
		father->voisin[c1]->special[b] = father;
		if ( father->voisin[c1]->voisin[b]->s[U] == p )
								father->voisin[c1]->voisin[b] = father;
		plusbeaufils(father->fils[c1], father->voisin[c1] );
		destroy(father->fils[c1],p);
		poubelle_noeud(father->fils[c1]);
		father->fils[c1] = NULL;
		init_Star(father,a1,b1,c1);
		a = b1 ;
		b = c1 ;
		c = a1 ;
	} while ( father->s[a] != first );
	init_Star( (noeud) NULL, 0, 0, 0);
}

void delaunay_tree::reinsertion( reinserer n, DT_item p)
{
	listenoeud l;
	noeud u,v,w,ww,w1,w2;
	DT_item x1,x2;
	int i,j,k,a,b,c,a1,b1,c1;
	star s1,s2;

	/* traitement des noeuds decroches */
	for ( l=n->decroches; l; l = l->next )
		if ( l->item->pere ) {
		/* il faut trouver un nouveau beau-pere et tous les voisinages 
                   par cette arete */
			u = l->item->pere;
			for( i=U; i<=W; i++) if (u->fils[i] == l->item) break;
			l->item->bpere = u->special[i];
			l->item->creation[U] = u->special[i];
			l->item->special[U] = u->special[i];
			if ( l->item->voisin[U]->detruire == flag )
								l->item->voisin[U] = u->special[i];
			beaufils(l->item, u->special[i]);
			for( j=U; j<=W; j++)
				if (	(l->item->bpere->s[j] != l->item->s[V] )