quadedge.c

关于限制性四叉树实现算法的

// Inserts a new site into a CDT. This is basically the Guibas-Stolfi
// incremental site insertion algorithm, except that is does not flip
// constraining edges (in fact, it does not even test them.)
// Returns NIL(Edge) if insertion has failed; otherwise, returns an edge
// whose origin is the new site.
{
	Edge *e = Locate(x);
	if (e == NIL(Edge))
	    return e;

	if (coincide(x, e->Org2d(), dist))
		return e;
	if (coincide(x, e->Dest2d(), dist))
	    return e->Sym();

	Boolean hasLeft = hasLeftFace(e);
	Boolean hasRight = hasRightFace(e);

	if (!hasLeft && !hasRight) {
		Warning("InsertSite: edge does not have any faces");
		return NIL(Edge);
	}

	// x should be inside a face adjacent to e:
	Boolean onEdge = OnEdge(x, e);
	Boolean insideLeft  = hasLeft && (onEdge || LeftOf(x, e)) &&
	                      RightOf(x, e->Onext()) && RightOf(x, e->Dprev());
	Boolean insideRight = hasRight && (onEdge || RightOf(x, e)) &&
	                      LeftOf(x, e->Oprev()) && LeftOf(x, e->Dnext());
	if (!insideLeft && !insideRight) {
		Warning("InsertSite: point not in a face adjacent to edge");
		return NIL(Edge);
	}	                    

	if (insideLeft  && coincide(x, e->Onext()->Dest2d(), dist))
		return e->Lprev();
	if (insideRight && coincide(x, e->Oprev()->Dest2d(), dist))
		return e->Dnext();

	// Now we know, x lies within the quadrilateral whose diagonal is e
	// (or a triangle, if e only has one adjacent face). We also know x
	// is not one of e's endpoints.

	if (onEdge) {

		// snap x to e, and check for coincidence:
		Point2d xx = snap(x, e->Org2d(), e->Dest2d());
		if (coincide(xx, e->Org2d(), dist))
			return e;
		if (coincide(xx, e->Dest2d(), dist))
			return e->Sym();

		if (hasRight && hasLeft) {
			// has two faces
			Edge *a, *b, *c, *d;
			a = e->Oprev();
			b = e->Dnext();
			c = e->Lnext();
			d = e->Lprev();
			SplitEdge(e, x);
			Connect(e, e->Lprev());
			Connect(e->Oprev(), e->Sym());
			FixEdge(a);
			FixEdge(b);
			FixEdge(c);
			FixEdge(d);
		} else {
			// has only one face
			if (hasRight)
		        e = e->Sym();
			Edge *c = e->Lnext();
			Edge *d = e->Lprev();
			SplitEdge(e, x);
			Connect(e, e->Lprev());
			FixEdge(c);
			FixEdge(d);
		}
		startingEdge = e->Sym();
		return startingEdge;
	}

	// x is not on e, should be in face to the left of e:
	if (!insideLeft) {
		Warning("InsertSite: point is not to the left of edge");
		return NIL(Edge);
	}

	// x should be strictly inside the left face:
	if (OnEdge(x, e->Onext()) || OnEdge(x, e->Dprev())) {
		Warning("InsertSite: point is not strictly inside face");
		return NIL(Edge);
	}

	// Now, hopefully, x is strictly inside the left face of e,
	// so everything should proceed smoothly from this point on.

   	Point2d* first = e->Org();
	Edge* base = MakeEdge(e->Org(), new Point2d(x));
	Splice(base, e);
	// use this edge as the starting point next time:
	startingEdge = base->Sym();
	do {
		base = Connect(e, base->Sym());
		e = base->Oprev();
	} while (e->Dprev() != startingEdge);

	// Examine suspect edges to ensure that the Delaunay condition
	// is satisfied.
   	do {
		Edge* t = e->Oprev();
		if (!e->isConstrained() &&
			InCircle(e->Org2d(), t->Dest2d(), e->Dest2d(), x)) {
				Swap(e);
				e = e->Oprev();
			}
		else if (e->Lprev() == startingEdge)  // no more suspect edges
			return startingEdge;
		else  // pop a suspect edge
		    e = e->Onext()->Lprev();
	} while (TRUE);
}

void Mesh::InsertEdge(const Point2d& a, const Point2d& b)
// Inserts a constraining edge into a CDT.
{
	Edge *ea, *eb;
	Point2d aa, bb;
	if (ea = InsertSite(a))
	    aa = ea->Org2d();
	if (eb = InsertSite(b))
	    bb = eb->Org2d();
	if (ea == NIL(Edge) || eb == NIL(Edge)) {
		Warning("InsertEdge: could not insert endpoints of edge");
		return;
	}

	startingEdge = ea;
	if ((ea = Locate(aa)) == NIL(Edge)) {
		Warning("InsertEdge: could not locate an endpoint");
	    return;
	}
	if (!(aa == ea->Org2d()))
	    ea = ea->Sym();
	if (!(aa == ea->Org2d())) {
		Warning("InsertEdge: point a is not an endpoint of ea");
		return;
	}

	if (aa == bb) {
		Warning("InsertEdge: both ends map to same vertex");
		return;
	}

	// aa and bb are vertices in the mesh; our goal is to connect
	// them by a sequence of one or more edges.

	Line ab(aa, bb);
	Vector2d dd = bb - aa;
	Edge *t, *ne;
	Point2d *last = eb->Org(); 

	while (ea->Org() != last) {

		// Set ea to the first edge to the right of (or aligned with)
		// the segment (a,b), by moving ccw around the origin of ea:
		t = ea;
		do {
			if ((ea->Dest2d() == ab && ((ea->Dest2d() - aa)|dd) > 0))
				break;
			if (ea->Onext()->Dest2d() == ab &&
				((ea->Onext()->Dest2d() - aa)|dd) > 0) {
					ea = ea->Onext();
					break;
				}
			if (!cw(ea->Dest2d(), bb, aa) && cw(ea->Onext()->Dest2d(), bb, aa))
			    break;
			ea = ea->Onext();
			if (ea == t) {
				Warning("InsertEdge: infinite loop");
				return;
			}
		} while (TRUE);

		// check to see if an edge is already there:
		if (ea->Dest2d() == ab) {
			ea->Constrain();
			aa = ea->Dest2d();
			if (aa == bb)
				return;
			ab.set(aa, bb);
			dd = bb - aa;
			ea = ea->Sym()->Onext();
			continue;
		}

		t = ea;

		while (TRUE) {
			if (t->Lnext()->Dest2d() == ab) {
				if (t->Lnext()->Lnext()->Lnext() == ea) {
					// edge is already there
					t->Lnext()->Lnext()->Constrain();
					ea = t->Lnext()->Sym();
					break;
				} else {
					// need a new edge
					ne = Connect(t->Lnext(), ea); ne->Constrain();
					ea = t->Lnext()->Sym();
					Triangulate(ne->Lnext());
					Triangulate(ne->Oprev());
					break;
				}
			}
			if (t->Lnext()->Dest2d() < ab) {
				// edges cross
				if (!t->Lnext()->isConstrained()) { 
					DeleteEdge(t->Lnext());
				} else {
					// the crossing edge is also constrained
					// compute and insert the intersection
					Point2d x = Intersect(t->Lnext(), ab);
					// split t->Lnext() into two at x
					SplitEdge(t->Lnext(), x);
					// connect to the new vertex:
					ne = Connect(t->Lnext(), ea);
					ne->Constrain();
					ea = t->Lnext()->Sym();
					Triangulate(ne->Lnext());
					Triangulate(ne->Oprev());
					break;
				}
			} else {
				t = t->Lnext();
			}
		}
	}
}

Boolean BruteForceLocate = FALSE;

Edge* Mesh::Locate(const Point2d& x)
// The goals of this routine are as follows:
// - if x coincides with one or more vertices in the mesh, return an edge
//   whose origin is one such vertex;
// - otherwise, return the edge nearest to x oriented ccw around a face
//   containing x (up to numerical round-off); 
// In the event of failure, we pass the search to BruteForceLocate.
{
	if (::BruteForceLocate)
		return BruteForceLocate(x);

	Edge* e = startingEdge;
	int iterations = 0;

	while (TRUE) {
		iterations++;

		if (iterations > numEdges()) {
			Warning("Locate: infinite loop");
			return BruteForceLocate(x);
		}

		// first, examine the endpoints for coincidence:
		if (x == e->Org2d())
			return startingEdge = e;
		if (x == e->Dest2d())
			return startingEdge = e->Sym();
		
		// x isn't one of the endpoints.
		// orient e s.t. x is not to its right:
		if (RightOf(x, e))
			e = e->Sym();
		
		// x is not to the right of e.
		// now we need to test if x is in the left face of e:
		if (hasLeftFace(e)) {
			
			// does x coincide with the third vertex?
			if (x == e->Lprev()->Org2d())
				return startingEdge = e->Lprev();
			
			if (LeftOf(x, e->Onext())) {
				e = e->Onext();
				continue;
			}
			if (LeftOf(x, e->Dprev())) {
				e = e->Dprev();
				continue;
			}
			
			// x is inside the left face of e!
			// return the edge nearest to x:
			Line l1(e->Org2d(), e->Dest2d());
			Line l2(e->Dest2d(), e->Lprev()->Org2d());
			Line l3(e->Lprev()->Org2d(), e->Org2d());
			switch (imin(fabs(l1.eval(x)), fabs(l2.eval(x)), fabs(l3.eval(x)))) {
			  case 0:
				startingEdge = e; break;
			  case 1:
				startingEdge = e->Lnext(); break;
			  case 2:
				startingEdge = e->Lprev(); break;
			}
			return startingEdge;
		}
		
		// there is no left face!
		if (OnEdge(x, e))
			return startingEdge = e->Sym();

		// ok, let's try to get closer using the edges of the right face
		if (RightOf(x, e->Oprev())) {
			e = e->Oprev()->Sym();
			continue;
		}
		if (RightOf(x, e->Dnext())) {
			e = e->Dnext()->Sym();
			continue;
		}

		Warning("Locate: I'm stuck and I can't move on");
		return BruteForceLocate(x);
	}	
}

Edge* Mesh::BruteForceLocate(const Point2d& x)
// Same as Locate, but uses a brute force exhaustive search.
{
	Edge *e;

	for (LlistPos p = edges.first(); !edges.isEnd(p); p = edges.next(p)) {

		e = ((QuadEdge *)edges.retrieve(p))->edges();
		// first, examine the endpoints for coincidence:
		if (x == e->Org2d())
			return startingEdge = e;
		if (x == e->Dest2d())
			return startingEdge = e->Sym();

		// x isn't one of the endpoints.
		// orient e s.t. x is not to its right:
		if (RightOf(x, e))
			e = e->Sym();

		// x is not to the right of e.
		// now we need to test if x is in the left face of e:
		if (hasLeftFace(e)) {

			// does x coincide with the third vertex?
			if (x == e->Lprev()->Org2d())
				return startingEdge = e->Lprev();

			if (!RightOf(x, e->Lprev()) && !RightOf(x, e->Lnext())) {
				// x is in this face!
				// return the edge nearest to x:
				Line l1(e->Org2d(), e->Dest2d());
				Line l2(e->Dest2d(), e->Lprev()->Org2d());
				Line l3(e->Lprev()->Org2d(), e->Org2d());
				switch(imin(fabs(l1.eval(x)), fabs(l2.eval(x)), fabs(l3.eval(x)))) {
				  case 0:
					startingEdge = e; break;
				  case 1:
					startingEdge = e->Lnext(); break;
				  case 2:
					startingEdge = e->Lprev(); break;
				}
				return startingEdge;
			} else
				continue;
		} 

		// there is no left face!
		if (OnEdge(x, e))
			return startingEdge = e->Sym();
		else
			continue;

	}

	// if we're here, it must be because we failed :-{
	Warning("Locate: Could not locate using brute force");
	return NIL(Edge);
}


/************************* Mesh Traversal Routines **************************/

#define UNVISITED 0
#define VISITED   1

inline void Mesh::MarkEdges(Edge *e)
{
	register Edge *t = e;
	do {
		t->mark = VISITED;
		t = t->Onext();
	} while (t != e);
}

void Mesh::ApplyEdges( void (*f)( void *arg, void *edge, Boolean ), void *arg )
// Applies the function f to every edge in the edge-list. The function
// takes application dependent arguments pointed to by arg, a pointer
// to the edge, and a flag that specifies whether or not the edge is
// constrained. This last argument allows the application to handle
// constrained edges differently.
{
	register QuadEdge *qp;
	register LlistPos p;

	for (p = edges.first(); !edges.isEnd(p); p = edges.next(p)) {
	    qp = (QuadEdge *)edges.retrieve(p);
	    f(arg, qp->edges(), qp->isConstrained()); 
	} 	
}

void Mesh::ApplyVertices( void (*f)( void *arg, void *vertex ), void *arg )
// Applies the function f to every vertex in the mesh. This is done
// similarly to ApplyEdges, but care must be taken to apply only once
// to each vertex. Every time f is applied to a vertex, all edges
// originating at that vertex are marked as VISITED. When we traverse
// the edge-list we only apply f to the origins of UNVISITED edges.
// After the traversal is complete, we go over the edges once more to
// reset the marks.
{
	register Edge *e;
	register LlistPos p;

	for (p = edges.first(); !edges.isEnd(p); p = edges.next(p)) {
	    e = ((QuadEdge *)edges.retrieve(p))->edges();
		if (e->mark == UNVISITED) {
			f(arg, e->Org());
			MarkEdges(e);
		}
		if (e->Sym()->mark == UNVISITED) {
			f(arg, e->Sym()->Org());
			MarkEdges(e->Sym());
		}
	} 	

	for (p = edges.first(); !edges.isEnd(p); p = edges.next(p)) {
	    e = ((QuadEdge *)edges.retrieve(p))->edges();
		e->mark = e->Sym()->mark = UNVISITED;
	}
}