Delaunay三角剖分

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的内存空间和一个外接矩形（该矩形盒子用来确定虚拟三角形）
init_delaunay函数如下，它是一个OpenCV函数，是一个包含一些OpenCV函数的函数包。

我们知道三角剖分是对散点集进行处理的，我们知道了散点集就可以获得点集的三角剖分。

如何传入（插入）散点集呢？这些点必须是32位浮点型，并通过下面的方式插入点：
转换为CvPoint2D32f的两种方法：
1）通过宏cvPoint2D32f(double x,double y)
2）通过cxtype.h下的cvPointTo32f(CvPoint point)函数将整形点方便的转换为32位浮点型。

当可以通过输入点（散点集）得到Delaunay三角剖分后，接下来，我们用一下两个函数设置和清除相关的Voronoi划分：
CvSubdiv2D结构如下：
[cpp]view plaincopy
[cpp]view plaincopy
平面划分是将一个平面分割为一组不重叠的、能够覆盖整个平面的区域。

结构CvSubdiv2D描述了建立在二维点集上的划分结构，其中点集互相连接且构成平面图形，该图形通过结合一些无线连接外部划分点（称为凸形点）的边缘，将一个平面用按照其边缘划分成很多小区域。

对于每一个划分操作，都有一个对偶划分与之对应。

对偶的意思是小区域与点（划分的顶点）变换角色，即在对偶划分中，小区域被当做一个顶点（以下称为虚拟点）而原始的划分顶点被当做小区域。

如下图所示，原始的划分用实线表示，而对偶划分用虚线表示。

OpenCV使用Delaunay算法将平面分割成小的三角形区域（该三角形确保包括所有的分割点）开始不断迭代完成。

在这种情况下，对偶划分就是输入的二维点集的Voronoi图表。

这种划分可以用于对一个平面进行三维分段变换、形态变换、平面点的快速定位以及建立特定的图结构（如NNG,RNG）。

CvQuadEdge2D
CvQuadEdge2D结构包含了两个Delaunay点和两个Vorionoi点以及连接它们的边缘（假设Voronoi点和边缘已经由函数计算出来，通过上面的函数：cvCalSubdivVoronoi2D(subdiv)）。

CvQuadEdge2D定义平面划分中的Quad-edge（四方边缘结构），其结构如下：
四方边缘结构是平面划分的基元，其中包括四条边缘（e,eRot,以及他们的反方向），如下图所示：
CvSubdiv2DPoint
CvSubdiv2DPoint结构包含Delaunay边缘及其相连的顶点。

其结构定义如下：
[cpp]view plaincopy
边缘的遍历
Subdiv2DRotateEdge函数
功能：函数Subdiv2DRotateEdge根据输入的边缘返回四方边缘结构中的一条边缘。

格式：
[cpp]view plaincopy
参数：
edge划分的边缘（不是四方边缘结构，即不是CvQuadEdge2D)
rotate 确定函数根据输入的边缘返回同一四方边缘结构中的哪条边缘，为下列值之一：1）0 输入边缘（如果e是输入边缘，则为e）。

2）1 旋转比那样（eRot）
3）2 逆边缘（e的反向边缘）。

4）3旋转比那样的反向边缘（eRot的反向边缘）。

cvSubdiv2DGetEdge函数
使用该函数我们可以遍历Delaunay图。

该函数返回与输入边缘相关的边缘。

格式：
参数：
edge 划分的边缘（不是四方边缘结构）
type 确定函数返回哪条相关边缘，为下列值之一：
CV_NEXT_AROUND_ORG 边缘原点的下一条（eOnext，如果e是输入边）。

CV_NEXT_AROUND_DST 边缘顶点的下一条(eDnext）
CV_PREV_AROUND_ORG 边缘原点的前一条(eRnext的反向）
CV_PREV_AROUND_DST边缘终点的前一条（eLnext的反向）
CV_NEXT_AROUND_LEFT 左区域的下一条(eLnext) 或下一个左平面
CV_NEXT_AROUND_RIGHT 右区域的下一条(eRnext) 或下一个右平面
CV_PREV_AROUND_LEFT 左区域的前一条(eOnext的反向）或前一个左平面
CV_PREV_AROUND_RIGHT 右区域的前一条（eDnext的反向)或前一个右平面
来至边缘的点
我们可以通过下面的两个函数获得Delaunay或者Voronoi边缘的两个实际点：
下面是将CvSubdiv2DPoint点转换为更熟悉的点CvPoint2D32f 或者CvPoint:
如何从Delaunay/Voronoi细分得到第一条边或点？
有两种方法：1）使用一个外部点定位边缘或顶点2）遍历一系列点或边缘
方法一：使用外部点定位边缘或顶点
第一种方法是任取一点，然后在细分中定位该点。

该点不一定是三角剖分中点，而可以为任意点。

cvSubdiv2DLocate()函数填充三角形的边缘和顶点（如果必要）或者填充该点所在的Voronoi面，函数的声明如下：
请注意，这些不必是最接近的边缘或顶点，它们只需要在三角形上。

此函数的返回值按下列的方式说明点的位置：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_vtx[3]和CvQuadEdge2D*outer_gedges[3]来存储三个顶点和三条边。

如下：
[cpp]view plaincopy
[cpp]view plaincopy
确定凸包的外接三角形（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文件代码理解参考图：
程序一：
执行效果如下图所示：黑色线表示三角剖分、绿色线表示相应的Voronoi细分，顶点用红色表示。

程序二:
和程序一类似，只是加了根据一个点，然后获取该点在所在的剖分平面（即包含该点的三角形），然后遍历该三角形：[cpp]view plaincopy
执行效果：
上图中我们选择了图像的中心点，图像大小为600×600，故中心点为300×300，我们通过函数cvSubdiv2DLocate()得知该点落在了面内，然后我们该函数得到该面的一个边缘，我们有了边缘信息，如果我们帮它当做程序1上的参考图中的e，那么我就可以获得参考图中的所有边的信息，我们也可以根据这些边的信息获取的这些边的端点。

并注意一下控制台的输出信息：即first边的方向向下，second为方向向上,third是方向向下。

程序三：
[cpp]view plaincopy
执行结果：
上程序中，通过isGoodTri()函数可以有效地防止重画并且消除虚拟三角形的顶点。

我要说的：
1、程序的关键是通过上面**号强调的利用CvSeqReader遍历所有的Delaunay或者Voronoi.
2、Delaunay或者Voronoi之间的桥梁就是通过函数，cvSubdiv2DRotateEdge()实现边的旋转（切换）。

3、这个Delaunay是个有向图，还是双向的。

具体怎么实现的，没有较深的数据结构为基础是很难理解的，至少现在我还不是很懂具体的实现，关于它的使用还是摸索阶段。

OpenCV水很深，如果这东西只是为我们用，或者了解，而不是研究，我们没有必要纠结的太深。

4、文章中的错误，还请指正，一起学习。