Delaunay三角剖分.docx

Delaunay三角剖分来源：相关文章：OpenCV三角剖分的遍历和纹理映射：Delaunay三角剖分是1934年发明的将空间点连接为三角形，使得所有三角形中最小角最大的一个技术。如果你熟悉计算机图形学，你便会知道Delaunay三角剖分是变现三维形状的基础。如果我们在三维空间渲染一个，我们可以通过这个物体的投影来建立二维视觉图，并用二维Delaunay三角剖分来分析识别该物体，或者将它与实物相比较。Delaunay剖分是连接计算机视觉与计算机图形学的桥梁。然而使用OpenCV实现三角剖分的不足之处就是OpenCV只实现了二维的Delaunay剖分。如果我们能够对三维点进行三角剖分，也就是说构成立体视觉，那么我们可以在三维的计算机图形和计算机视觉进行无缝的转换。然而二维三角剖分通常用于计算机视觉中标记空间目标的特征或运动场景跟踪，目标识别，或两个不同的摄像机的场景匹配（如图从立体图像中获得深度信息）。下面内容摘自：1 三角剖分与Delaunay剖分的定义如何把一个离散几何剖分成不均匀的三角形网格，这就是离散点的三角剖分问题，散点集的三角剖分，对数值分析以及图形学来说，都是极为重要的一项处理技术。该问题图示如下：1.1 三角剖分定义【定义】三角剖分：假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段，E为e的集合。那么该点集V的一个三角剖分T=(V,E)是一个平面图G,该平面图满足条件：1、除了端点，平面图中的边不包含点集中的任何点。2、没有相交边。（边和边没有交叉点）3、平面图中所有的面都是三角面，且所有三角面的合集是散点集V的凸包。1.2 Delaunay三角剖分的定义在实际中运用的最多的三角剖分是Delaunay三角剖分，它是一种特殊的三角剖分。先从Delaunay边说起：【定义】Delaunay边：假设E中的一条边e(两个端点为a,b)，e若满足下列条件，则称之为Delaunay边：存在一个圆经过a,b亮点，圆内（注意是圆内，圆上最多三点共圆）不含点集V中任何其他的点，这一特性又称空圆特性。【定义】Delaunay三角剖分：如果点集V的一个三角剖分T只包含Delaunay边，那么该三角剖分称为Delaunay三角剖分。1.3 Delaunay三角剖分的准则要满足Delaunay三角剖分的定义，必须符合两个重要的准则：1、空圆特性：Delaunay三角网是唯一的（任意四点不能共圆），在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在。如下图所示：2、最大化最小角特性：在散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大。从这个意义上讲，Delaunay三角网是“最接近于规则化的”三角网。具体的说是在两个相邻的三角形构成凸四边形的对角线，在相互交换后，两个内角的最小角不再增大。如下图所示1.4 Delaunay三角剖分的特性以下是Delaunay剖分所具备的优异特性：1、最接近：以最接近的三点形成三角形，且各线段（三角行的边）皆不相交。2、唯一性：不论从区域何处开始构建，最终都将得到一致的结果。3、最优性：任意两个相邻三角形构成的凸四边形的对角线如何可以互换的话，那么两个三角形六个内角中最小角度不会变化。4、最规则：如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大。5、区域性：新增、删除、移动某一个顶点只会影响邻近的三角形。6、具有凸边形的外壳：三角网最外层的边界形成一个凸多边形的外壳。1.5局部最优化处理理论上为了构造Delaunay三角网，Lawson提出的局部优化过程LOP(Local Optimization Procedure),一般三角网经过LOP处理，即可确保为Delaunay三角网，其基本做法如下所示：1、将两个具有共同边的三角形合成一个多边形。2、以最大空圆准则作检查，看其第四个顶点是否在三角形的外接圆内。3、如果在，修正对角线即将对角线对调，即完成局部优化过程的处理。LOP处理过程如下图所示：2 Delaunay剖分的算法Delaunay剖分是一种三角剖分的标准，实现它有多种算法。由表可以看出，三角网生成法的时间效率最低，分治算法的时间效率最高，逐点插入法效率居中。由于区域生长法本质的缺陷，导致其效率受限，这种方法在80年代中期以后已经很少使用。分治算法时间效率相对较高，但是由于其递归执行，所以需要较大的内存空间，导致其空间效率较低。此外，分治法的数据处理及结果的优化需要的工作量也比较大。逐点插入算法实现简单，时间效率比较高，而运行占用的空间也较小，从时间效率和空间效率综合考虑，性价比最高，因而应用广泛。1977年，Lawson提出了逐点插入法建立Delaunay三角网的算法思想。之后Lee和Schachlter,Bowyer,Watson,Sloan,先后进行了发展和完善。他们的算法在初始化三角网的建立方法、定位点所在三角形的过程、以及插入的过程方面各具特点。2.1 Lawson算法逐点插入的Lawson算法是Lawson在1977提出的，该算法思路简单，易于编程实现。基本原理为：首先建立一个大的三角形或多边形，把所有数据点包围起来，向其中插入一点，该点与包含它的三角形三个顶点相连，形成三个新的三角形，然后逐个对它们进行空外接圆检测，同时用Lawson设计的局部最优化过程LOP进行优化，集通过交换对角线的方法来保证所形成的三角网为Delaunay三角网。上述基于散点的构网算法理论严密、唯一性好，网格满足空圆特性，较为理想。由其逐点插入的构网过程可知，遇到非Delaunay边时，通过删除调整，可以构造形成新的Delaunay边。在完成构网后，增加新点时，无需对所有点进行重新构网，只需要对新点的影响三角形范围进行局部联网，且局部联网的方法简单易行。同样，点的删除、移动也可快速动态地进行。但在实际应用当中，这种构网算法当点集较大时网速度也较慢，如果点集范围是非凸区域或者存在内环，则会产生非法三角形。2.2 Bowyer-Watson算法目前采用逐点插入方式生成的Delaunay三角网的算法主要基于Bowyer-Watson算法，Bowyer-Watson算法的主要步骤如下：1）建立初始三角网格：针对给定的点集V,找到一个包含该点集的矩形R,我们称R为辅助窗口。连接R的任意一条对角线，形成两个三角形，作为初始Delaunay三角网格。2）逐点插入：假设目前已经有一个Delaunay三角网格T,现在在它里面再插入一个点P,需要找到该点P所在的三角形。从P所在的三角形开始，搜索该三角形的邻近三角形，并进行空外接圆检测。找到外接圆包含点P的所有的三角形并删除这些三角形，形成一个包含P的多边形空腔，我们称之为Delaunay空腔。然后连接P与Delaunay腔的每一个顶点，形成新的Delaunay三角网格。3）删除辅助窗口R:重复步骤2）,当点集V中所有点都已经插入到三角形网格中后，将顶点包含辅助窗口R的三角形全部删除。在这些步骤中，快速定位点所在的三角形、确定点的影响并构建Delaunay腔的过程是每插入一个点都会进行的。随着点数的增加，三角形数目增加很快，因此缩短这两个过程的计算时间，是提高算法效率的关键。算法执行图示如下：、3 OpenCV中的Delaunay和Voronoi细分学习这部分，也是一个头疼的问题，要理解需要较好的数据结构作为基础。由于自己对数据结构也是敬畏三分，所以下面理解不免有误。OpenCV中现实的Delaunay三角剖分应该是Bowyer-Watson算法。3.1创建一个Delaunay或Voronoi细分我们需要存储Delaunay的内存空间和一个外接矩形（该矩形盒子用来确定虚拟三角形）1. / STORAGE AND STRUCTURE FOR DELAUNAY SUBDIVISION /存储和结构 for三角剖分 2. / 3. CvRect rect = 0, 0, 600, 600 ; /Our outer bounding box /我们的外接边界盒子 4. CvMemStorage* storage; /Storage for the Delaunay subdivsion /用来存储三角剖分 5. storage = cvCreateMemStorage(0); /Initialize the storage /初始化存储器 6. CvSubdiv2D* subdiv; /The subdivision itself / 细分 7. subdiv = init_delaunay( storage, rect); /See this function below /函数返回CvSubdiv类型指针 init_delaunay函数如下，它是一个OpenCV函数，是一个包含一些OpenCV函数的函数包。1. /INITIALIZATION CONVENIENCE FUNCTION FOR DELAUNAY SUBDIVISION /为三角剖分初始化便利函数 2. / 3. CvSubdiv2D* init_delaunay(CvMemStorage* storage,CvRect rect) 4. CvSubdiv2D* subdiv = cvCreateSubdiv2D(CV_SEQ_KINDSUBDIV2D,sizeof(*subdiv),sizeof(CvSubdiv2DPoint),sizeof(CvQuadEdge2D),storage);/为数据申请空间 5. cvInitSubdivDelaunay2D( subdiv, rect ); /rect sets the bounds 6. return subdiv;/返回申请空间的指针 7. 我们知道三角剖分是对散点集进行处理的，我们知道了散点集就可以获得点集的三角剖分。如何传入（插入）散点集呢？这些点必须是32位浮点型，并通过下面的方式插入点：1. CvPoint2D32f fp; /This is our point holder/这是我们点的持有者（容器） 2. for( i = 0; i as_many_points_as_you_want; i+ ) 3. / However you want to set points /如果我们的点集不是32位的，在这里我们将其转为CvPoint2D32f，如下两种方法。 4. / 5. fp = your_32f_point_listi; 6. cvSubdivDelaunay2DInsert( subdiv, fp ); 7. 转换为CvPoint2D32f的两种方法：1）通过宏cvPoint2D32f(double x,double y)2）通过cxtype.h下的cvPointTo32f(CvPoint point)函数将整形点方便的转换为32位浮点型。当可以通过输入点（散点集）得到Delaunay三角剖分后，接下来，我们用一下两个函数设置和清除相关的Voronoi划分：1. cvCalcSubdivVoronoi2D( subdiv ); / Fill out Voronoi data in subdiv /在subdiv中填充Vornoi的数据 2. cvClearSubdivVoronoi2D( subdiv ); / Clear the Voronoi from subdiv/从subdiv中清除Voronoi的数据 CvSubdiv2D结构如下：1. #define CV_SUBDIV2D_FIELDS() 2. CV_GRAPH_FIELDS() 3. int quad_edges; 4. int is_geometry_valid; 5. CvSubdiv2DEdge recent_edge; 6. CvPoint2D32f topleft; 7. CvPoint2D32f bottomright; 8. typedef struct CvSubdiv2D 9. 10. CV_SUBDIV2D_FIELDS() 11. 12. CvSubdiv2D; cpp view plaincopy1. #define CV_GRAPH_FIELDS() 2. CV_SET_FIELDS() /* set of vertices */ 3. CvSet *edges; /* set of edges */ cpp view plaincopy1. #define CV_SET_FIELDS() 2. CV_SEQUENCE_FIELDS() /*inherits from #CvSeq CvSeq */ 3. struct CvSetElem* free_elems; /*list of free nodes */ 平面划分是将一个平面分割为一组不重叠的、能够覆盖整个平面的区域。结构CvSubdiv2D描述了建立在二维点集上的划分结构，其中点集互相连接且构成平面图形，该图形通过结合一些无线连接外部划分点（称为凸形点）的边缘，将一个平面用按照其边缘划分成很多小区域。对于每一个划分操作，都有一个对偶划分与之对应。对偶的意思是小区域与点（划分的顶点）变换角色，即在对偶划分中，小区域被当做一个顶点（以下称为虚拟点）而原始的划分顶点被当做小区域。如下图所示，原始的划分用实线表示，而对偶划分用虚线表示。OpenCV使用Delaunay算法将平面分割成小的三角形区域（该三角形确保包括所有的分割点）开始不断迭代完成。在这种情况下，对偶划分就是输入的二维点集的Voronoi图表。这种划分可以用于对一个平面进行三维分段变换形态变换、平面点的快速 定位以及建立特定的图结构（如NNG,RNG）。CvQuadEdge2DCvQuadEdge2D结构包含了两个Delaunay点和两个Vorionoi点以及连接它们的边缘（假设Voronoi点和边缘已经由函数计算出来，通过上面的函数：cvCalSubdivVoronoi2D(subdiv)）。CvQuadEdge2D定义平面划分中的Quad-edge（四方边缘结构），其结构如下：1. /* one of edges within quad-edge, lower 2 bits is index (0.3) and upper bits are quad-edge pointer */2. typedef long CvSubdiv2DEdge; /四方边缘结构中的一条边缘，低两位表示该边缘的索引号，其他高位表示边缘指针。 3. 4. /*quad-edge structure fields*/四方边缘的结构场 5. #define CV_QUADEDGE2D_FIELDS() 6. int flags; 7. struct CvSubdiv2DPoint*pt4; 8. CvSubdiv2DEdge next4; 9. typedef struct CvQuadEdge2D 10. 11. CV_QUADEDGE2D_FIELDS() 12. 13. CvQuadEdge2D; 四方边缘结构是平面划分的基元，其中包括四条边缘（e,eRot,以及他们的反方向），如下图所示：CvSubdiv2DPointCvSubdiv2DPoint结构包含Delaunay边缘及其相连的顶点。其结构定义如下：cpp view plaincopy1. #define CV_SUBDIV2D_POINT_FIELDS() 2. int flags; 3. CvSubdiv2DEdge first; 4. CvPoint2D32f pt; 5. int id; /This integer can be used to index auxillary data asscoiated with each vertex of the planar subdivison. 6. #define CV_SUBDIV2D_VIRTUAL_POINT_FLAG (1 pt; / to 32f point 3. CvPoint pt = cvPointFrom32f(pt32f); / to an integer point 如何从Delaunay/Voronoi细分得到第一条边或点？有两种方法：1）使用一个外部点定位边缘或顶点 2）遍历一系列点或边缘方法一：使用外部点定位边缘或顶点第一种方法是任取一点，然后在细分中定位该点。该点不一定是三角剖分中点，而可以为任意点。cvSubdiv2DLocate()函数填充三角形的边缘和顶点（如果必要）或者填充该点所在的Voronoi面，函数的声明如下：1. CvSubdiv2DPointLocation cvSubdiv2DLocate( 2. CvSubdiv2D* subdiv, 3. CvPoint2D32f pt, 4. CvSubdiv2DEdge* edge,/要填充的边缘 5. CvSubdiv2DPoint* vertex = NULL/如果需要，则填充顶点 6. ); 请注意，这些不必是最接近的边缘或顶点，它们只需要在三角形上。此函数的返回值按下列的方式说明点的位置：1）CV_PTLOC_INSIDE 点落入某些面；*edge将包含该面的一个边缘。2）CV_PTLOC_ON_ENCODE 点落于边缘；*edge含有这个边缘。3）CV_PTLOC_VERTEX 该点与一个细分顶点重合；*vertex将包含该顶点的指针。4）CV_PTLOC_OUTSIDE_RECT 该点处于细分参考矩形之外；该函数返回后不填充指针。5）CV_PTLOC_ERROR 输入变量无效。方法二：遍历一系列点或边缘外接三角形（虚拟）的三个顶点和三个边的获取：1）首先我们要有一个建立的Delaunay细分。2）我们还需要调用cvCalcSubdivVoronoi2D(subdiv)函数计算相关的Voronoi划分。3）然后我们就能用CvSubdiv2DPoint *outer_vtx3和CvQuadEdge2D*outer_gedges3来存储三个顶点和三条边。如下：cpp view plaincopy1. CvSubdiv2DPoint* outer_vtx3; 2. for( i = 0; i 3; i+ ) 3. outer_vtxi = 4. (CvSubdiv2DPoint*)cvGetSeqElem( (CvSeq*)subdiv, I ); 5. cpp view plaincopy1. CvQuadEdge2D* outer_ qedges3; 2. for( i = 0; i edges), I ); 5. 确定凸包的外接三角形（Bounding triangle）或边缘并遍历凸包我们回想一下，我们通过调用cvInitSubdivDelaunay2D(subdiv,rect)来初始化Delaunay三角剖分。在这种情况下，下面的论述成立：1、如果边缘的起点和终点都在矩形之外，那么此边缘也在细分的外接三角形上。（即如果一个边缘的两个点出了矩形边界，那么该边缘为虚拟三角形的边缘）2、如果边缘的一端在矩形内，一段在矩形外，那么矩形边界内的点落在凸集上，凸集上的每个点与虚拟外接三角形的两顶点相连，这两边相继出现。（英语原文：If you are on an edge with one point inside and one point outside the rect bounds,then the point in bounds is on the convex hull of the set。）（注释：Learning OpenCV中第343页，将In bounds翻译为在边界上，这里个人感觉应该是在翻译为在边界内，为了方便理解，个人猜测如下图，可能不对：）从上述第二个条件可知，我们可以使用cvSubdiv2DNextEdge()，移动到第一条边上，这条边的dst端点在边界内。两个端点都在边界内的这条边在点集的凸包上，因此记下那个点和边。一旦它在凸包上，那么你可以如下遍历所有顶点。（切记参考英文原版，我感觉下面的步骤应该参考具体操作理解。）1、将凸包遍历一周后，通过cvSubdiv2DRotateEdge(CvSubdiv2DEdge edge,0)函数移动到凸包的下一条边。2、接着，两次调用宏cvSubdiv2DNextEdge()就到了凸包的下一条边。跳转到第一步。使用实例：1、使用函数cvSubdiv2DLoacate()遍历Delaunay三角剖分的边。2、使用函数cvFindNearestPoint2D()函数找到距离输入最近的点。3、使用函数draw_subdiv_facet()函数逐步遍历Voronoi面，这个函数是组合函数，实现参考OpenCV的344页。4、使用CvSeqReader逐步遍历Delaunay或者Voronoi边。5、找到了三角剖分的所有顶点，我们就可以利用内敛宏cvTriangleArea()函数来计算剖分的面积。三角剖分例程序：在OpenCV的sample的c目录下：delaunay.c文件代码理解参考图：程序一：1. #include 2. #include 3. #include opencv2/highgui/highgui.hpp 4. #include 5. #include 6. #include 7. using namespace std; 8. using namespace cv; 9. static void help( void ) 10. 11. printf(nThis program demostrates iterative construction ofn/这个程序阐述了delaunay剖分和voronoi细分的迭代构造 12. delaunay triangulation and voronoi tesselation.n 13. It draws a random set of points in an image and then delaunay triangulates them.n/在图像上画出一些随机点，然后进行delaunay三角剖分 14. Usage: n 15. ./delaunay n 16. nThis program builds the traingulation interactively, you may stop this process byn 17. hitting any key.n);/迭代构造三角剖分，如果像停止，则按任意键 18. 19. 20. static CvSubdiv2D* init_delaunay( CvMemStorage* storage,/初始化三角剖分结构，为其分配单元 21. CvRect rect ) 22. 23. CvSubdiv2D* subdiv;/三角剖分的数据单元 24. 25. subdiv = cvCreateSubdiv2D( CV_SEQ_KIND_SUBDIV2D, sizeof(*subdiv), 26. sizeof(CvSubdiv2DPoint), 27. sizeof(CvQuadEdge2D), 28. storage ); 29. cvInitSubdivDelaunay2D( subdiv, rect ); 30. 31. return subdiv; 32. 33. 34. 35. static void draw_subdiv_point( IplImage* img, CvPoint2D32f fp, CvScalar color )/画出三角剖分的顶点 36. 37. cvCircle( img, cvPoint(cvRound(fp.x), cvRound(fp.y),5, color, CV_FILLED, 8, 0 ); 38. 39. 40. 41. static void draw_subdiv_edge( IplImage* img, CvSubdiv2DEdge edge, CvScalar color )/画出三角剖分的边 42. 43. CvSubdiv2DPoint* org_pt;/源顶点 44. CvSubdiv2DPoint* dst_pt;/目地顶点 45. CvPoint2D32f org 46. CvPoint2D32f dst; 47. CvPoint iorg, idst; 48. 49. org_pt = cvSubdiv2DEdgeOrg(edge);/通过边获取顶点 50. dst_pt = cvSubdiv2DEdgeDst(edge); 51. 52. if( org_pt & dst_pt )/如果两个端点不为空 53. 54. org = org_pt-pt; 55. dst = dst_pt-pt; 56. 57. iorg = cvPoint( cvRound( org.x ), cvRound( org.y ); 58. idst = cvPoint cvRound( dst.x ), cvRound( dst.y ); 59. 60. cvLine( img, iorg, idst, color, 1, CV_AA, 0 ); 61. 62. 63. 64. 65. static void draw_subdiv( IplImage* img, CvSubdiv2D* subdiv, 66. CvScalar delaunay_color, CvScalar voronoi_color )/画出剖分和细分 67. 68. CvSeqReader reader; 69. int i, total = subdiv-edges-total/边的数量 70. int elem_size = subdiv-edges-elem_size;/边的大小 71. coutedges).name()edges), &reader, 0 );/使用CvSeqReader遍历Delaunay或者Voronoi边 74. 75. for( i = 0; i total; i+ ) 76. 77. CvQuadEdge2D* edge = (CvQuadEdge2D*)(reader.ptr); 78. 79. if( CV_IS_SET_ELEM( edge ) 80. 81. / draw_subdiv_edge( img, (CvSubdiv2DEdge)edge + 1, voronoi_color ); 82. draw_subdiv_edge( img, (CvSubdiv2DEdge)edge, delaunay_color ); 83. 84. 85. CV_NEXT_SEQ_ELEM( elem_size, reader ); 86. 87. 88. 89. 90. 91. static void locate_point( CvSubdiv2D* subdiv, CvPoint2D32f fp, IplImage* img,/遍历三角剖分的边 92. CvScalar active_color ) 93. 94. CvSubdiv2DEdge e; 95. CvSubdiv2DEdge e0 = 0; 96. CvSubdiv2DPoint* p = 0; 97. 98. cvSubdiv2DLocate( subdiv, fp, &e0, &p ); 99. 100. if( e0 ) 101. 102. e = e0; 103. do 104. 105. draw_subdiv_edge( img, e, active_color ); 106. e = cvSubdiv2DGetEdge(e,CV_NEXT_AROUND_LEF