cdt.tar.gz_四叉树_四叉树算法_限制性四叉树

/***************************************************************************** ** ** quadedge.C: routines for QuadEdge and Mesh classes and for contstruction ** of constrained Delaunay triangulations (CDTs). ** ** Copyright (C) 1995 by Dani Lischinski ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. ** ** This program is distributed in the hope that it will be useful, ** but WITHOUT ANY WARRANTY; without even the implied warranty of ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ** GNU General Public License for more details. ** ** You should have received a copy of the GNU General Public License ** along with this program; if not, write to the Free Software ** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. ** ******************************************************************************/ #include "quadedge.h" /*********************** Basic Topological Operators ************************/ Edge* Mesh::MakeEdge(Boolean constrained = FALSE) { QuadEdge *qe = new QuadEdge(constrained); qe->p = edges.insert(edges.first(), qe); return qe->edges(); } Edge* Mesh::MakeEdge(Point2d *a, Point2d *b, Boolean constrained = FALSE) { Edge *e = MakeEdge(constrained); e->EndPoints(a, b); return e; } void Splice(Edge* a, Edge* b) // This operator affects the two edge rings around the origins of a and b, // and, independently, the two edge rings around the left faces of a and b. // In each case, (i) if the two rings are distinct, Splice will combine // them into one; (ii) if the two are the same ring, Splice will break it // into two separate pieces. // Thus, Splice can be used both to attach the two edges together, and // to break them apart. See Guibas and Stolfi (1985) p.96 for more details // and illustrations. { Edge* alpha = a->Onext()->Rot(); Edge* beta = b->Onext()->Rot(); Edge* t1 = b->Onext(); Edge* t2 = a->Onext(); Edge* t3 = beta->Onext(); Edge* t4 = alpha->Onext(); a->next = t1; b->next = t2; alpha->next = t3; beta->next = t4; } void Mesh::DeleteEdge(Edge* e) { // Make sure the starting edge does not get deleted: if (startingEdge->QEdge() == e->QEdge()) { Warning("Mesh::DeleteEdge: attempting to delete starting edge"); // try to recover: startingEdge = (e != e->Onext()) ? e->Onext() : e->Dnext(); Assert(startingEdge->QEdge() != e->QEdge()); } // remove edge from the edge list: QuadEdge *qe = e->QEdge(); edges.remove(edges.prev(qe->p)); delete qe; } QuadEdge::~QuadEdge() { // If there are no other edges from the origin or the destination // the corresponding data pointers should be deleted: if (e[0].data != NIL(Point2d) && e[0].next == e) delete e[0].data; if (e[2].data != NIL(Point2d) && e[2].next == (e+2)) delete e[2].data; // Detach edge from the rest of the subdivision: Splice(e, e->Oprev()); Splice(e->Sym(), e->Sym()->Oprev()); } Mesh::~Mesh() { QuadEdge *qp; for (LlistPos p = edges.first(); !edges.isEnd(p); p = edges.next(p)) { qp = (QuadEdge *)edges.retrieve(p); delete qp; } } /************* Topological Operations for Delaunay Diagrams *****************/ Mesh::Mesh(const Point2d& a, const Point2d& b, const Point2d& c) // Initialize the mesh to the triangle defined by the points a, b, c. { Point2d *da = new Point2d(a); Point2d *db = new Point2d(b); Point2d *dc = new Point2d(c); Edge* ea = MakeEdge(da, db, TRUE); Edge* eb = MakeEdge(db, dc, TRUE); Edge* ec = MakeEdge(dc, da, TRUE); Splice(ea->Sym(), eb); Splice(eb->Sym(), ec); Splice(ec->Sym(), ea); startingEdge = ec; } Mesh::Mesh(const Point2d& a, const Point2d& b, const Point2d& c, const Point2d& d) // Initialize the mesh to the Delaunay triangulation of // the quadrilateral defined by the points a, b, c, d. // NOTE: quadrilateral is assumed convex. { Boolean InCircle(const Point2d&, const Point2d&, const Point2d&, const Point2d&); Point2d *da = new Point2d(a); Point2d *db = new Point2d(b); Point2d *dc = new Point2d(c); Point2d *dd = new Point2d(d); Edge* ea = MakeEdge(da, db, TRUE); Edge* eb = MakeEdge(db, dc, TRUE); Edge* ec = MakeEdge(dc, dd, TRUE); Edge* ed = MakeEdge(dd, da, TRUE); Splice(ea->Sym(), eb); Splice(eb->Sym(), ec); Splice(ec->Sym(), ed); Splice(ed->Sym(), ea); // Split into two triangles: if (InCircle(c, d, a, b)) { // Connect d to b: startingEdge = Connect(ec, eb); } else { // Connect c to a: startingEdge = Connect(eb, ea); } } // The following typedefs are necessary to avoid stupid warnings // on some compilers: typedef Point2d* Point2dPtr; typedef Edge* EdgePtr; Mesh::Mesh( int numVertices, double *bdryVertices ) // Initialize the mesh to the Delaunay triangulation of the // convex polygon defined by the coordinates in the bdryVertices // array. The number of vertices (coordinate pairs) is specified // via numVertices. // NOTE: polygon is assumed convex. { Point2d orig, dest; Edge *e0, *e1; register int i; Assert(numVertices >= 3); Point2d **verts = new Point2dPtr[numVertices + 1]; Edge **edges = new EdgePtr[numVertices + 1]; // Create all vertices: for (i = 0; i < numVertices; i++) { verts[i] = new Point2d(bdryVertices[2*i], bdryVertices[2*i+1]); } verts[numVertices] = verts[0]; // Create all edges: for (i = 0; i < numVertices; i++) { edges[i] = MakeEdge(verts[i], verts[i+1], TRUE); } edges[numVertices] = edges[0]; // Connect edges together: for (i = 0; i < numVertices; i++) { Splice(edges[i]->Sym(), edges[i+1]); } // Triangulate: Triangulate(edges[0]); // Initialize starting edge: startingEdge = edges[0]; delete[] verts; delete[] edges; } Edge* Mesh::Connect(Edge* a, Edge* b) // Add a new edge e connecting the destination of a to the // origin of b, in such a way that all three have the same // left face after the connection is complete. // Additionally, the data pointers of the new edge are set. { Edge* e = MakeEdge(a->Dest(), b->Org()); Splice(e, a->Lnext()); Splice(e->Sym(), b); return e; } void Swap(Edge* e) // Essentially turns edge e counterclockwise inside its enclosing // quadrilateral. The data pointers are modified accordingly. { Edge* a = e->Oprev(); Edge* b = e->Sym()->Oprev(); Splice(e, a); Splice(e->Sym(), b); Splice(e, a->Lnext()); Splice(e->Sym(), b->Lnext()); e->EndPoints(a->Dest(), b->Dest()); } Point2d snap(const Point2d& x, const Point2d& a, const Point2d& b) { if (x == a) return a; if (x == b) return b; Real t1 = (x-a) | (b-a); Real t2 = (x-b) | (a-b); Real t = MAX(t1,t2) / (t1 + t2); return ((t1 > t2) ? ((1-t)*a + t*b) : ((1-t)*b + t*a)); } void Mesh::SplitEdge(Edge *e, const Point2d& x) // Shorten edge e s.t. its destination becomes x. Connect // x to the previous destination of e by a new edge. If e // is constrained, x is snapped onto the edge, and the new // edge is also marked as constrained. { Point2d *dt; if (e->isConstrained()) { // snap the point to the edge before splitting: dt = new Point2d(snap(x, e->Org2d(), e->Dest2d())); } else dt = new Point2d(x); Edge *t = e->Lnext(); Splice(e->Sym(), t); e->EndPoints(e->Org(), dt); Edge *ne = Connect(e, t); if (e->isConstrained()) ne->Constrain(); } /*************** Geometric Predicates for Delaunay Diagrams *****************/ inline Real TriArea(const Point2d& a, const Point2d& b, const Point2d& c) // Returns twice the area of the oriented triangle (a, b, c), i.e., the // area is positive if the triangle is oriented counterclockwise. { return (b[X] - a[X])*(c[Y] - a[Y]) - (b[Y] - a[Y])*(c[X] - a[X]); } Boolean InCircle(const Point2d& a, const Point2d& b, const Point2d& c, const Point2d& d) // Returns TRUE if the point d is inside the circle defined by the // points a, b, c. See Guibas and Stolfi (1985) p.107. { Real