Delaunay三角化

点集的三角剖分（Triangulation），对数值分析（比如有限元分析）以及图形学来说，都是极为重要的一项预处理技术。

尤其是Delaunay三角剖分，由于其独特性，关于点集的很多种几何图都和Delaunay三角剖分相关，如Voronoi图，EMST树，Gabriel图等。

Delaunay剖分所具备的优异特性：

1.最接近：以最近的三点形成三角形，且各线段(三角形的边)皆不相交。

2.唯一性：不论从区域何处开始构建，最终都将得到一致的结果。

3.最优性：任意两个相邻三角形形成的凸四边形的对角线如果可以互换的话，那么两个三角形六个内角中最小的角度不会变大。

4.最规则：如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大。

5.区域性：新增、删除、移动某一个顶点时只会影响临近的三角形。

6.具有凸多边形的外壳：三角网最外层的边界形成一个凸多边形的外壳。

概念及定义
二维实数域(二维平面)上的三角剖分

定义1：假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段, E为e的集合。
那么该点集V的一个三角剖分T=(V,E)是一个平面图G，该平面图满足条件：
1.除了端点，平面图中的边不包含点集中的任何点。
2.没有相交边。
3.平面图中所有的面都是三角面，且所有三角面的合集就是点集V的凸包。
那什么是Delaunay三角剖分呢?不过是一种特殊的三角剖分罢了。从Delaunay边说起。

Delaunay边
定义2：假设E中的一条边e（两个端点为a,b），e若满足下列条件，则称之为Delaunay边：
存在一个圆经过a,b两点，圆内不含点集V中任何的点，这一特性又称空圆特性。

Delaunay三角剖分

定义3：如果点集V的一个三角剖分T只包含Delaunay边，那么该三角剖分称为Delaunay三角剖分。

定义4：假设T为V的任一三角剖分，则T是V的一个Delaunay三角剖分，当前仅当T中的每个三角形的外接圆的内部不包含V中任何的点。

定义5：V的Voronoi图是由多边形区域的集合（有些区域可能不是闭合的），该区域仅含点集中的一个点v，区域中的任何位置到v的距离都比该位置到点集中其它所有点的距离短。

由Voronoi图和Delaunay三角剖分的关系，可以引出另一个Delaunay三角剖分的定义：
定义6：将Voronoi图相邻区域（共边的区域）中的点连接起来构成的图，称为Delaunay三角剖分。
如下图：


概念部分到此，下面看看怎么求Delaunay三角剖分。

计算Delaunay三角剖分

问题1：计算二维Delaunay三角剖分
问题输入：二维实数域上的点集V
问题输出：Delaunay三角剖分DT = (V, E).

算法

由不同的定义对应有不同的算法。用得较多的是基于定义3或4的算法。
目前常用的算法又分为好几种，被不同的家伙发现。什么扫描线法（Sweepline），随机增量法（Incremental），分治法（Divide and Conquer）

c++实现：（转载）

[cpp]  view plain  copy

1. //adapted by the example of leanring opencv by crazy_007  
2. //adapted by the OpenCV2.0\samples\c\delaunay.c  
3. //2010-4-22  
4. #include "StdAfx.h"  
5. #include <cv.h>  
6. #include <highgui.h>  
7. void draw_subdiv_edge( IplImage* img, CvSubdiv2DEdge edge, CvScalar color ); //为了查看代码方便，不然需调整调用函数顺序  
8.   
9. /* the script demostrates iterative construction of 
10.    delaunay triangulation and voronoi tesselation */  
11.   
12. CvSubdiv2D* init_delaunay( CvMemStorage* storage,  
13.                            CvRect rect )  
14. {  
15.     CvSubdiv2D* subdiv;  
16.   
17.     subdiv = cvCreateSubdiv2D( CV_SEQ_KIND_SUBDIV2D, sizeof(*subdiv),  
18.                                sizeof(CvSubdiv2DPoint),  
19.                                sizeof(CvQuadEdge2D),  
20.                                storage );  
21.     cvInitSubdivDelaunay2D( subdiv, rect );//矩形确定的边界  
22.   
23.     return subdiv;  
24. }  
25.   
26. //draw subdiv point  
27. void draw_subdiv_point( IplImage* img, CvPoint2D32f fp, CvScalar color )  
28. {  
29.     cvCircle( img, cvPoint(cvRound(fp.x), cvRound(fp.y)), 5, color, CV_FILLED, 8, 0 );//画fp为圆心，5为半径的实心圆表示delaunay顶点  
30. }  
31. //use an external point to locate an edge or vertex or step around the edges of a delaunay tirangle  
32. //画出delaunay 顶点  
33. void locate_point( CvSubdiv2D* subdiv, CvPoint2D32f fp, IplImage* img,  
34.                    CvScalar active_color )  
35. {  
36.     CvSubdiv2DEdge e;  
37.     CvSubdiv2DEdge e0 = 0;  
38.     CvSubdiv2DPoint* p = 0;  
39.   
40.     cvSubdiv2DLocate( subdiv, fp, &e0, &p );//使用一个外部的点定位边缘或顶点  
41.                                             //该函数填充三角形的边缘和顶点或者填充该点所处在的Voronoi面  
42.     if( e0 )  
43.     {  
44.         e = e0;  
45.         do  
46.         {  
47.             draw_subdiv_edge( img, e, active_color );//调用下面函数：画出红色直线  
48.             e = cvSubdiv2DGetEdge(e,CV_NEXT_AROUND_LEFT);//遍历Delaunay图：返回左区域的下一条的边缘  
49.         }  
50.         while( e != e0 );  
51.     }  
52.   
53.     draw_subdiv_point( img, fp, active_color );//调用上面的函数：实现在该点处画半径为5的圆  
54. }  
55.   
56.   
57.   
58.   
59. /*************                 分割线                ************************************/  
60. //draw subdiv edge  
61. void draw_subdiv_edge( IplImage* img, CvSubdiv2DEdge edge, CvScalar color )  
62. {  
63.     CvSubdiv2DPoint* org_pt;  
64.     CvSubdiv2DPoint* dst_pt;  
65.     CvPoint2D32f org;  
66.     CvPoint2D32f dst;  
67.     CvPoint iorg, idst;  
68.   
69.     org_pt = cvSubdiv2DEdgeOrg(edge);//Delaunay或者Voronoi边缘的原始点  
70.     dst_pt = cvSubdiv2DEdgeDst(edge);//Delaunay或者Voronoi边缘的终点  
71.   
72.     if( org_pt && dst_pt )  
73.     {  
74.         org = org_pt->pt;  
75.         dst = dst_pt->pt;  
76.   
77.         iorg = cvPoint( cvRound( org.x ), cvRound( org.y ));  
78.         idst = cvPoint( cvRound( dst.x ), cvRound( dst.y ));  
79.   
80.         cvLine( img, iorg, idst, color, 1, CV_AA, 0 );//画红色直线  
81.     }  
82. }  
83.   
84. //draw subdiv:遍历所有的Delaunay边  
85. void draw_subdiv( IplImage* img, CvSubdiv2D* subdiv,  
86.                   CvScalar delaunay_color, CvScalar voronoi_color )  
87. {  
88.     CvSeqReader reader;//使用cvSeqReader逐步遍历边：获得细分结构  
89.     int i, total = subdiv->edges->total;  
90.     int elem_size = subdiv->edges->elem_size;  
91.   
92.     cvStartReadSeq( (CvSeq*)(subdiv->edges), &reader, 0 );//initialize reader of the sequence  
93.   
94.     for( i = 0; i < total; i++ )//total是边数目  
95.     {  
96.         CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr);  
97.                                                             //CvQuadEdge2D平面划分中的Quad-edge(四方边缘结构):四个边缘 (e, eRot(红色) 以及它们的逆（绿色）  
98.        if( CV_IS_SET_ELEM( edge ))  
99.         {  
100.             draw_subdiv_edge( img, (CvSubdiv2DEdge)edge + 1, voronoi_color ); //不知如何理解(CvSubdiv2DEdge)edge + 1   
101.                 //书中P346：voronoi_edge=(CvSubdiv2DEdge)edge + 1  
102.             //直接采用数组位移法进行各种边的对应的(即edge+1),  
103.             //cvSubdiv2DRotateEdge((CvSubdiv2DEdge)edge,1)=(CvSubdiv2DEdge)edge+1  
104.             //参考网址为：http://tech.ddvip.com/2007-12/119897724239781.html  
105.             draw_subdiv_edge( img, (CvSubdiv2DEdge)edge, delaunay_color ); //调用上面的子函数  
106.         }  
107.   
108.         CV_NEXT_SEQ_ELEM( elem_size, reader );  
109.     }  
110. }  
111.   
112.   
113.   
114. /*************                 分割线                ************************************/  
115. //draw the voronoi facet:遍历Voronoi面  
116. void draw_subdiv_facet( IplImage* img, CvSubdiv2DEdge edge )  
117. {  
118.     CvSubdiv2DEdge t = edge;  
119.     int i, count = 0;  
120.     CvPoint* buf = 0;  
121.   
122.     // count number of edges in facet  
123.     do  
124.     {  
125.         count++;  
126.         t = cvSubdiv2DGetEdge( t, CV_NEXT_AROUND_LEFT );//返回左区域的下一条的边缘  
127.     } while (t != edge );  
128.   
129.     buf = (CvPoint*)malloc( count * sizeof(buf[0]));  
130.   
131.     // gather points  
132.     t = edge;  
133.     for( i = 0; i < count; i++ )  
134.     {  
135.         CvSubdiv2DPoint* pt = cvSubdiv2DEdgeOrg( t );//获得边缘的起点  
136.         if( !pt ) break;  
137.         buf[i] = cvPoint( cvRound(pt->pt.x), cvRound(pt->pt.y));//点记录在buf中  
138.         t = cvSubdiv2DGetEdge( t, CV_NEXT_AROUND_LEFT );  
139.     }  
140.   
141.     if( i == count )  
142.     {  
143.         CvSubdiv2DPoint* pt = cvSubdiv2DEdgeDst( cvSubdiv2DRotateEdge( edge, 1 ));//获得边缘的终点  
144.         cvFillConvexPoly( img, buf, count, CV_RGB(rand()&255,rand()&255,rand()&255), CV_AA, 0 );//一次只能画一个多边形，而且只能画凸多边形  
145.         cvPolyLine( img, &buf, &count, 1, 1, CV_RGB(0,0,0), 1, CV_AA, 0);//一次调用中绘制多个多边形  
146.         draw_subdiv_point( img, pt->pt, CV_RGB(0,0,0));//画圆  
147.     }  
148.     free( buf );  
149. }  
150. //draw & paint voronoi graph  
151. void paint_voronoi( CvSubdiv2D* subdiv, IplImage* img )  
152. {  
153.     CvSeqReader reader;  
154.     int i, total = subdiv->edges->total;  
155.     int elem_size = subdiv->edges->elem_size;  
156.   
157.     cvCalcSubdivVoronoi2D( subdiv );//计算Voronoi图表的细胞结构  
158.   
159.     cvStartReadSeq( (CvSeq*)(subdiv->edges), &reader, 0 );  
160.   
161.     for( i = 0; i < total; i++ )  
162.     {  
163.         CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr);  
164.   
165.         if( CV_IS_SET_ELEM( edge ))  
166.         {  
167.             CvSubdiv2DEdge e = (CvSubdiv2DEdge)edge;  
168.             // left  
169.             draw_subdiv_facet( img, cvSubdiv2DRotateEdge( e, 1 ));//调用上面的子函数  
170.   
171.             // right  
172.             draw_subdiv_facet( img, cvSubdiv2DRotateEdge( e, 3 ));  
173.         }  
174.   
175.         CV_NEXT_SEQ_ELEM( elem_size, reader );  
176.     }  
177. }  
178. /*************                 分割线                ************************************/  
179.   
180. void run(void)  
181. {  
182.     char win[] = "source";  
183.     int i;  
184.     CvRect rect = { 0, 0, 600, 600 };//外界边界矩形大小  
185.     CvMemStorage* storage;//为delaunay申请内存空间  
186.     CvSubdiv2D* subdiv;//细分  
187.     IplImage* img;  
188.     CvScalar active_facet_color, delaunay_color, voronoi_color, bkgnd_color;  
189.   
190.     active_facet_color = CV_RGB( 255, 0, 0 );  
191.     delaunay_color = CV_RGB( 0,0,0);  
192.     voronoi_color = CV_RGB(0, 180, 0);  
193.     bkgnd_color = CV_RGB(255,255,255);  
194.   
195.     img = cvCreateImage( cvSize(rect.width,rect.height), 8, 3 );//创建白色背景的图像  
196.     cvSet( img, bkgnd_color, 0 );  
197.   
198.     cvNamedWindow( win, 1 );//窗口名字为"source"  
199.   
200.     storage = cvCreateMemStorage(0);//初始化内存空间  
201.     subdiv = init_delaunay( storage, rect );//initialization convenience function for delaunay subdivision：调用子函数确定矩形边界  
202.   
203.     printf("Delaunay triangulation will be build now interactively.\n"  
204.            "To stop the process, press any key\n\n");  
205.   
206.     for( i = 0; i < 20; i++ )  
207.     {  
208.         CvPoint2D32f fp = cvPoint2D32f( (float)(rand()%(rect.width-10)+5),  
209.                                         (float)(rand()%(rect.height-10)+5));//This is our point holder  
210.   
211.         locate_point( subdiv, fp, img, active_facet_color );//调用函数  
212.         cvShowImage( win, img );  
213.   
214.         if( cvWaitKey( 100 ) >= 0 ) //等待600ms  
215.             break;  
216.   
217.         cvSubdivDelaunay2DInsert( subdiv, fp );//向Delaunay三角测量中插入一个点  
218.         cvCalcSubdivVoronoi2D( subdiv );//计算Voronoi图表的细胞结构  
219.         cvSet( img, bkgnd_color, 0 ); // 给一个对象全部元素赋值:void cvSet( CvArr* arr, CvScalar value, const CvArr* mask=NULL );  
220.         draw_subdiv( img, subdiv, delaunay_color, voronoi_color );//调用子函数：cvSeqReader逐步遍历边来获得细分结构  
221.         cvShowImage( win, img );  
222.   
223.         if( cvWaitKey(100) >= 0 )  
224.             break;  
225.     }  
226.   
227.     cvSet( img, bkgnd_color, 0 );  
228.     paint_voronoi( subdiv, img );//调用子函数：  
229.     cvShowImage( win, img );  
230.   
231.     cvWaitKey(0);  
232.   
233.     cvReleaseMemStorage( &storage );  
234.     cvReleaseImage(&img);  
235.     cvDestroyWindow( win );  
236. }  
237.   
238. int main( int argc, char** argv )  
239. {  
240.     run();  
241.     return 0;  
242. }  

原文地址： http://blog.csdn.net/kuikuijia/article/details/46282959