逐点插入法建立Delaunay三角网

<P>逐点插入法进行构网，分析该算法中制约运行速度的因素，在三角网拓扑关系、三角形的查找、LOP算法（Local Optimization Procedure）等方面进行了优化处理，使之有了较高的执行效率。<BR>　　1　数据结构<BR>　　在采用逐点插入进行Delaunay三角剖分的过程中，存在大量的查询操作，利用数据结构提供三角形之间的拓扑信息，能够有效地提高算法效率。边结构应包含点与三角形的信息，三角形结构应包含点与边的信息，再考虑到构网过程中动态的数据更新，算法中采用了双向链表结构，以方便于剖分过程中新旧边（三角形）的添加、删除工作。<BR>　　2　算法原理<BR>　　逐点插入法是Lawson在1977年提出的，该算法思路简单，易于编程实现。基本原理为：对于已有的三角网络，向其中插入一点，该点与包含它的三角形三个顶点相连，形成三个新的三角形，然后逐个对它们进行空外接圆检测，通过交换对角线的方法来保证所形成的三角网为Delaunay三角网。<BR>　　3　算法基本步骤</P>
<P>      a: 点在三角形内，b:点在三角形边上。<BR>　　在按a进行剖分时，添加了三条新边及三个新三角形NT，删除了旧三角形T。同样，在b的情况中，记录点所在的边E，根据拓扑关系，找出该边的左右相邻三角形T1，T2，添加四条新边和四个新三角形NT，删除T1，T2以及边E。<BR>　　构网的关键是对新剖分出的三角形用LOP算法进行优化，其基本过程为：新三角形与周围的三角形构成共用同一条对角线的四边形，逐个对四边形中的两个三角形进行空外接圆检测，如果满足空外接圆准则，则略过。如果不满足，则用另一对角线替换现有对角线，在交换对角线后，又会产生相应的新四边形，继续进行空外接圆检测。这种方法可连续进行,直到全部满足空外接准则为止。这个过程是随着数据点P的逐次插入，与三角剖分同时进行的，可以通过递归调用优化程序来实现。<BR>　　</P>

<img src="images/post/smile/dvbbs/em01.gif" />

我以前做过Delaunay三角网方面的东西，谁要这方面的资料可以发邮件给我：<a href="mailtzlzhang1982@163.com" target="_blank" >zlzhang1982@163.com</A>

<P>这个算法20万点内可直接建网，可以的，以用在一个公路三维项目上了</P>
<P>当然，还可以分区建模，然后合并。</P>

<P>我有生成等高线程序：不过没有进行三次样条光滑，当散点多的时候速度不是很好。指1700点以上的时候。</P>

<P>这个程序我做过了，可是当点太多的时候效果也不是太少。</P>
<P>算法的这个东东，到处都有，关键实现却没有，其实以上判断都是用4点共圆函数搞定的。</P>
<P>有兴趣加QQ:575167467</P>

<P>提供的思路不错，谢了！</P>
<img src="images/post/smile/dvbbs/em01.gif" /><img src="images/post/smile/dvbbs/em02.gif" />

不好意思，忘了分页了。希望斑竹帮忙，分一下页，另外还有一点说明资料，欢迎交流。

<P>#include <geom2d.h><BR>class QuadEdge;<BR>class Edge {<BR>friend QuadEdge;<BR>friend void Splice(Edge*, Edge*);<BR>private:<BR>int num;<BR>Edge *next;<BR>Point2d *data;<BR>public:<BR>Edge() { data = 0; }<BR>Edge* Rot();</P>
<P>Edge* invRot();<BR>Edge* Sym();<BR>Edge* Onext();<BR>Edge* Oprev();<BR>Edge* Dnext();<BR>Edge* Dprev();<BR>Edge* Lnext();<BR>Edge* Lprev();<BR>Edge* Rnext();<BR>Edge* Rprev();<BR>Point2d* Org();<BR>Point2d* Dest();<BR>const Point2d; Org2d() const;<BR>const Point2d; Dest2d() const;<BR>void EndPoints(Point2d*, Point2d*);<BR>QuadEdge* Qedge() { return (QuadEdge *)(this - num); }<BR>};<BR>class QuadEdge {<BR>friend Edge *MakeEdge();<BR>private:<BR>Edge e[4];<BR>public:<BR>QuadEdge();<BR>};<BR>class Subdivision {<BR>private:<BR>Edge *startingEdge;<BR>Edge *Locate(const Point2d;);<BR>public:<BR>Subdivision(const Point2d;, const Point2d;, const Point2d;);<BR>void InsertSite(const Point2d;);<BR>void Draw();<BR>};<BR>inline QuadEdge::QuadEdge()<BR>{<BR>e[0].num = 0, e[1].num = 1, e[2].num = 2, e[3].num = 3;<BR>e[0].next = ;(e[0]); e[1].next = ;(e[3]);<BR>e[2].next = ;(e[2]); e[3].next = ;(e[1]);<BR>}<BR>/************************* Edge Algebra *************************************/<BR>inline Edge* Edge::Rot()<BR>// Return the dual of the current edge, directed from its right to its left.<BR>{<BR>return (num < 3) ? this + 1 : this - 3;<BR>}<BR>inline Edge* Edge::invRot()<BR>// Return the dual of the current edge, directed from its left to its right.</P>
<P>{<BR>return (num > 0) ? this - 1 : this + 3;<BR>}<BR>inline Edge* Edge::Sym()<BR>// Return the edge from the destination to the origin of the current edge.<BR>{<BR>return (num < 2) ? this + 2 : this - 2;<BR>}<BR>inline Edge* Edge::Onext()<BR>// Return the next ccw edge around (from) the origin of the current edge.<BR>{<BR>return next;<BR>}<BR>inline Edge* Edge::Oprev()<BR>// Return the next cw edge around (from) the origin of the current edge.<BR>{<BR>return Rot()->Onext()->Rot();<BR>}<BR>inline Edge* Edge::Dnext()<BR>// Return the next ccw edge around (into) the destination of the current edge.<BR>{<BR>return Sym()->Onext()->Sym();<BR>}<BR>inline Edge* Edge::Dprev()<BR>// Return the next cw edge around (into) the destination of the current edge.<BR>{<BR>return invRot()->Onext()->invRot();<BR>}<BR>inline Edge* Edge::Lnext()<BR>// Return the ccw edge around the left face following the current edge.<BR>{<BR>return invRot()->Onext()->Rot();<BR>}<BR>inline Edge* Edge::Lprev()<BR>// Return the ccw edge around the left face before the current edge.<BR>{<BR>return Onext()->Sym();<BR>}<BR>inline Edge* Edge::Rnext()<BR>// Return the edge around the right face ccw following the current edge.<BR>{<BR>return Rot()->Onext()->invRot();<BR>}<BR>inline Edge* Edge::Rprev()</P>
<P>// Return the edge around the right face ccw before the current edge.<BR>{<BR>return Sym()->Onext();<BR>}<BR>/************** Access to data pointers *************************************/<BR>inline Point2d* Edge::Org()<BR>{<BR>return data;<BR>}<BR>inline Point2d* Edge::Dest()<BR>{<BR>return Sym()->data;<BR>}<BR>inline const Point2d; Edge::Org2d() const<BR>{<BR>return *data;<BR>}<BR>inline const Point2d; Edge::Dest2d() const<BR>{<BR>return (num < 2) ? *((this + 2)->data) : *((this - 2)->data);<BR>}<BR>inline void Edge::EndPoints(Point2d* or, Point2d* de)<BR>{<BR>data = or;<BR>Sym()->data = de;<BR>}<BR>/*********************** Basic Topological Operators ************************/<BR>Edge* MakeEdge()<BR>{<BR>QuadEdge *ql = new QuadEdge;<BR>return ql->e;<BR>}<BR>void Splice(Edge* a, Edge* b)<BR>// This operator affects the two edge rings around the origins of a and b,<BR>// and, independently, the two edge rings around the left faces of a and b.<BR>// In each case, (i) if the two rings are distinct, Splice will combine<BR>// them into one; (ii) if the two are the same ring, Splice will break it<BR>// into two separate pieces.<BR>// Thus, Splice can be used both to attach the two edges together, and<BR>// to break them apart. See Guibas and Stolfi (1985) p.96 for more details<BR>// and illustrations.<BR>{<BR>Edge* alpha = a->Onext()->Rot();<BR>Edge* beta = b->Onext()->Rot();</P>
<P>Edge* t1 = b->Onext();<BR>Edge* t2 = a->Onext();<BR>Edge* t3 = beta->Onext();<BR>Edge* t4 = alpha->Onext();<BR>a->next = t1;<BR>b->next = t2;<BR>alpha->next = t3;<BR>beta->next = t4;<BR>}<BR>void DeleteEdge(Edge* e)<BR>{<BR>Splice(e, e->Oprev());<BR>Splice(e->Sym(), e->Sym()->Oprev());<BR>delete e->Qedge();<BR>}<BR>/************* Topological Operations for Delaunay Diagrams *****************/<BR>Subdivision::Subdivision(const Point2d; a, const Point2d; b, const Point2d; c)<BR>// Initialize a subdivision to the triangle defined by the points a, b, c.<BR>{<BR>Point2d *da, *db, *dc;<BR>da = new Point2d(a), db = new Point2d(b), dc = new Point2d(c);<BR>Edge* ea = MakeEdge();<BR>ea->EndPoints(da, db);<BR>Edge* eb = MakeEdge();<BR>Splice(ea->Sym(), eb);<BR>eb->EndPoints(db, dc);<BR>Edge* ec = MakeEdge();<BR>Splice(eb->Sym(), ec);<BR>ec->EndPoints(dc, da);<BR>Splice(ec->Sym(), ea);<BR>startingEdge = ea;<BR>}<BR>Edge* Connect(Edge* a, Edge* b)<BR>// Add a new edge e connecting the destination of a to the<BR>// origin of b, in such a way that all three have the same<BR>// left face after the connection is complete.<BR>// Additionally, the data pointers of the new edge are set.<BR>{<BR>Edge* e = MakeEdge();<BR>Splice(e, a->Lnext());<BR>Splice(e->Sym(), b);<BR>e->EndPoints(a->Dest(), b->Org());<BR>return e;<BR>}<BR>void Swap(Edge* e)<BR>// Essentially turns edge e counterclockwise inside its enclosing<BR>// quadrilateral. The data pointers are modified accordingly.</P>
<P>{<BR>Edge* a = e->Oprev();<BR>Edge* b = e->Sym()->Oprev();<BR>Splice(e, a);<BR>Splice(e->Sym(), b);<BR>Splice(e, a->Lnext());<BR>Splice(e->Sym(), b->Lnext());<BR>e->EndPoints(a->Dest(), b->Dest());<BR>}<BR>/*************** Geometric Predicates for Delaunay Diagrams *****************/<BR>inline Real TriArea(const Point2d; a, const Point2d; b, const Point2d; c)<BR>// Returns twice the area of the oriented triangle (a, b, c), i.e., the<BR>// area is positive if the triangle is oriented counterclockwise.<BR>{<BR>return (b.x - a.x)*(c.y - a.y) - (b.y - a.y)*(c.x - a.x);<BR>}<BR>int InCircle(const Point2d; a, const Point2d; b,<BR>const Point2d; c, const Point2d; d)<BR>// Returns TRUE if the point d is inside the circle defined by the<BR>// points a, b, c. See Guibas and Stolfi (1985) p.107.<BR>{<BR>return (a.x*a.x + a.y*a.y) * TriArea(b, c, d) -<BR>(b.x*b.x + b.y*b.y) * TriArea(a, c, d) +<BR>(c.x*c.x + c.y*c.y) * TriArea(a, b, d) -<BR>(d.x*d.x + d.y*d.y) * TriArea(a, b, c) > 0;<BR>}<BR>int ccw(const Point2d; a, const Point2d; b, const Point2d; c)<BR>// Returns TRUE if the points a, b, c are in a counterclockwise order<BR>{<BR>return (TriArea(a, b, c) > 0);<BR>}<BR>int RightOf(const Point2d; x, Edge* e)<BR>{<BR>return ccw(x, e->Dest2d(), e->Org2d());<BR>}<BR>int LeftOf(const Point2d; x, Edge* e)<BR>{<BR>return ccw(x, e->Org2d(), e->Dest2d());<BR>}<BR>int OnEdge(const Point2d; x, Edge* e)<BR>// A predicate that determines if the point x is on the edge e.<BR>// The point is considered on if it is in the EPS-neighborhood<BR>// of the edge.<BR>{<BR>Real t1, t2, t3;<BR>t1 = (x - e->Org2d()).norm();</P>
<P>t2 = (x - e->Dest2d()).norm();<BR>if (t1 < EPS || t2 < EPS)<BR>return TRUE;<BR>t3 = (e->Org2d() - e->Dest2d()).norm();<BR>if (t1 > t3 || t2 > t3)<BR>return FALSE;<BR>Line line(e->Org2d(), e->Dest2d());<BR>return (fabs(line.eval(x)) < EPS);<BR>}<BR>/************* An Incremental Algorithm for the Construction of *************/<BR>/************************ Delaunay Diagrams *********************************/<BR>Edge* Subdivision::Locate(const Point2d; x)<BR>// Returns an edge e, s.t. either x is on e, or e is an edge of<BR>// a triangle containing x. The search starts from startingEdge<BR>// and proceeds in the general direction of x. Based on the<BR>// pseudocode in Guibas and Stolfi (1985) p.121.<BR>{<BR>Edge* e = startingEdge;<BR>while (TRUE) {<BR>if (x == e->Org2d() || x == e->Dest2d())<BR>return e;<BR>else if (RightOf(x, e))<BR>e = e->Sym();<BR>else if (!RightOf(x, e->Onext()))<BR>e = e->Onext();<BR>else if (!RightOf(x, e->Dprev()))<BR>e = e->Dprev();<BR>else<BR>return e;<BR>}<BR>}<BR>void Subdivision::InsertSite(const Point2d; x)<BR>// Inserts a new point into a subdivision representing a Delaunay<BR>// triangulation, and fixes the affected edges so that the result<BR>// is still a Delaunay triangulation. This is based on the<BR>// pseudocode from Guibas and Stolfi (1985) p.120, with slight<BR>// modifications and a bug fix.<BR>{<BR>Edge* e = Locate(x);<BR>if ((x == e->Org2d()) || (x == e->Dest2d())) // point is already in<BR>return;<BR>else if (OnEdge(x, e)) {<BR>e = e->Oprev();<BR>DeleteEdge(e->Onext());<BR>}<BR>// Connect the new point to the vertices of the containing<BR>// triangle (or quadrilateral, if the new point fell on an<BR>// existing edge.)</P>
<P>Edge* base = MakeEdge();<BR>base->EndPoints(e->Org(), new Point2d(x));<BR>Splice(base, e);<BR>startingEdge = base;<BR>do {<BR>base = Connect(e, base->Sym());<BR>e = base->Oprev();<BR>} while (e->Lnext() != startingEdge);<BR>// Examine suspect edges to ensure that the Delaunay condition<BR>// is satisfied.<BR>do {<BR>Edge* t = e->Oprev();<BR>if (RightOf(t->Dest2d(), e) ;;<BR>InCircle(e->Org2d(), t->Dest2d(), e->Dest2d(), x)) {<BR>Swap(e);<BR>e = e->Oprev();<BR>}<BR>else if (e->Onext() == startingEdge) // no more suspect edges<BR>return;<BR>else // pop a suspect edge<BR>e = e->Onext()->Lprev();<BR>} while (TRUE);<BR>}<BR>/*****************************************************************************/</P>
<P>老外的代码，我看不懂，希望高手指点！QQ:47827889</P>

可否再就数据结构，以及点的查找提供些资料！