光怪陆离的世界之Delaunay三角剖分和Voronoi图

时间:2022-07-25
本文章向大家介绍光怪陆离的世界之Delaunay三角剖分和Voronoi图，主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项，具有一定的参考价值，需要的朋友可以参考一下。

缘起

封面图是不是很酷炫? 该图的核心算法就是 Delaunay三角剖分. 这种低多边形的成像效果在现代游戏设计中越来越被喜欢，其中的低多边形都是由三角形组成的。于是我们来学习一下.

分析

首先，先来明确我们的问题

给出很多离散点，我们希望基于这个点集，将点集连接成一定大小的三角形，且分配要相对合理，力求能呈现出漂亮的三角化效果，就像上图那样.

为什么要是三角形，而不是四边形、五边形呢? 因为站在审美的角度，这种低多边形的成像效果在现代设计中越来越被喜欢. 其中的低多边形大都是由三角形组成的.

所以，自然就涉及到三角剖分，所以自然要给出三角剖分的定义（以下定义引自百度百科）.

【定义】三角剖分：假设V是 上的有限点集，称 V 的完全图的子图 T=(V,E) 是 V 的一个三角剖分，如果T是一个可平面图，而且满足 T 中的所有面都是三角面，且所有三角面的合集恰好是V的凸包

ps: 一些图论的概念

  • 完全图是一个无向图，其中每对不同的顶点之间都恰连有一条边相连
  • 可平面图是指 能将图在平面画出且不相交，缘起于电路板布线的设计.

讲清楚了什么叫三角剖分，那么自然就要面对剩下的一个问题：

什么叫做 分配要相对合理?

这就涉及到了 Delaunay 三角剖分，由俄国数学家 B.Delaunay于 1934 年提出.

【定义】Delaunay边：假设 E中的一条边e（两个端点为a,b），如果存在一个圆经过a,b两点，使得圆的严格内部不含V中任何其他的点，那么称 e 是一条 Delaunay边.

【定义】Delaunay三角剖分：如果 T 只包含Delaunay边，那么T被称为Delaunay三角剖分.

来张图直观体会一下三角剖分

上图左边的离散点集 V 的 三角剖分 就是右边.

可以证明 三角剖分 具备以下两个优秀的性质

  1. 空圆特性：三角剖分中的每个三角面的外接圆的严格内部不包含任何 V 中其他的点. 特别的，如果 V 中出现多点共圆（点的个数>3）, 则这种情况被称之为退化情况.
  2. 最大化最小角特性
  3. 唯一性：不论从区域何处开始构建，最终都将得到一致的结果。
  4. 区域性：新增、删除、移动某一个顶点时只会影响临近的三角形。
  5. 具有凸包的外壳：三角网最外层的边界形成一个凸多边形的外壳。

具体画图解释前两个性质.

大家可以看一下上面两幅图. 便知道什么叫做空圆特性，然后第二幅图中的两个三角形的最小的那个内角的度数一定大于第一幅图中的两个三角形的最小的内角的度数. 这就是最大化最小角特性.

所以 Delaunay 三角剖分给我们的直观感觉就是

  1. 空圆特性给我们的感觉是：每个三角面（或者叫三角网格）相对独立，因为任何一个三角网格的外接圆都不会接纳除了这三个点之外的点. 如果你将三角形的外接圆视作是 AT 力场（A.T Field）的话，则 Delaunay 三角剖分给我们展示的是一个有隔阂的世界.
  2. 最大化最小角特性给我们的感觉是：某种程度上保证了每个三角网格的丰满, 而避免了狭长三角形的产生，因为从美学角度，狭长的三角形并不是很惹人喜爱.

所以Delaunay三角剖分其实并不是一种算法，它只是给出了一个好的三角剖分的定义

为了方便，除非特别声明，否则下文提及的三角剖分指的就是 Delaunay三角剖分

三角剖分和其他问题的联系. 其中最著名的问题就是 Voronoi 图（也有文献称之为Thiessen 多边形，即泰森多边形），Voronoi 图是一种将平面分裂成许许多多的多边形区域（称之为瓦片），每块瓦片内部有一个点称之为该瓦片的生成点，所有处于瓦片的严格内部的点距离该瓦片的生成点的距离将严格小于它到其他瓦片的生成点的距离，而两块相邻瓦片的接壤边上的点到两块瓦片的生成点的距离是相等的（所以两块瓦片的生成点的连线将垂直平分两块瓦片的接壤边）.

ps: 这里提一嘴，水立方的外围就是Voronoi图.

基于泰森多边形的这些特点，泰森多边形是特别适合于做5G基站规划的. 为什么这么说呢? 因为实际生产生活中是有如下需求的——在网络整体规划中，我们常遇到的是点状分布的基站，缺乏整体连续的面状性，这使得我们无法直观的估算单站覆盖的面积

如上图所示，红点是现有基站的位置. 然后我们想估算每个基站需要覆盖的面积，因为这样的话，我们就能合理买入合适功率的基站设备，你总不能一个基站我们只需要它覆盖很小的城市面积，然后你买一个非常大功率的基站设备，那岂不是资金的浪费? 那么怎么估算一个红点需要覆盖的城市面积呢? 我们只需要让这些红点成为生成点，然后生成上图对应的泰森多边形就行了. 就像下图这样

一旦泰森多边形，或者说 Voronoi 图被构建，则估算多边形的面积就是轻而易举的事情. 事实上，现在已经有了现成的软件来做这件事情，那就是 MapInfo（美国MapInfo公司的桌面地理信息系统软件，又是老美的公司，唉！）

MapInfo 软件能使得基站覆盖面积的估算变得就像操作简单的 excel 那样简单. 而且，泰森多边形给出了一种估算，或者说衡量点集的分布类型——或者说聚集类型.

只需要计算泰森多边形面积的变异系数（CV）即可. 变异系数在统计学中的定义是标准差除以期望. 如果 CV 很大，则表明点集分布是一小撮一小撮这种，如果 CV 很小，表示点集的分布是均匀的.

当然，靠近边界的泰森多边形的面积很大程度上受到边界的影响. 而这个边界是人为划定的. 例如，你要考察南昌市的基站的覆盖面积情况，你的边界最好是划定在南昌市，如果划定到江西省的话，则 CV 值将变得很大. 这显然是不合理的.

说了这么多，Voronoi图 和三角剖分的关系是什么呢? 你看如下一张图就会顿悟它俩是对偶关系.

A、B、C、D、E、F、G、H、I 这 9个点（也就是前面说的点集 V）你可以视作是基站. 然后我们希望得到这9个离散点的 Voronoi图，也就是上图中的粗黑实线多边形. 如果我们能得到这9个点的 三角剖分, 也就是上图中的细黑实线三角形的话，则还记得我们之前垂直平分线的说法吗? 所以立马我们就知道 a 是三角网格 AEC 的外心（因为 aF=aA=aC）. 因此从一个点集的 三角剖分获取该点集的 Voronoi图 是轻而易举的事情. 只需要获取一个顶点，例如A，以A 为其中一个顶点的所有相邻三角网格（按照顺时针或者逆时针获取），例如上图就是 ACF、AFG、AGH、AHB、ABC ，然后计算这些个三角形的外心坐标，例如上图就是 a、b、c、d、e，则abcde就是 Voronoi图的一个多边形. 所以我们只需要遍历 V 中所有点集，对每个点执行一次上面的程序，得到一个Voronoi图 的多边形即可.

这里顺便说一下如何从A顺时针或者逆时针获取相邻的三角形. 我们任意枚举一个三角形，例如 ACF, 然后选择 AC, 然后以 AC 为边的三角形只有 ACF 和另一个，然后我们选取另一个，那就是 AFG， 然后再选取 AG, 然后选取以 AG 为边的另一个三角形，那就是AGH，然后再选取另一个以 AH 为边的三角形，那就是 AHB， 然后再选取另一个以 AB 为边的三角形 ，那就是ABC，然后发现另一个以 AC 为边的三角形是 ACF, 你看，又回来了，说明我们已经得到了 A 为顶点的一圈的相邻的三角形. 然后进一步就得到了 Voronoi图 的一个多边形.

最后，我们来研究一下 三角剖分 的具体算法. 因为前面说了，三角剖分 并不是一个实际的算法，而仅仅是一个较美的三角剖分的定义而已. 所以历史上有很多 三角剖分 的实现算法. 有分治算法、Lawson算法、Bowyer-Watson算法等，本文介绍最为常用的 Bowyer-Watson算法.

该算法的思想是极为简明易懂的——和凸包类似，也是逐点插入, 不断更新的算法（也就是增量算法）. 以四个点A、B、C、D举例，我们想建立该四点的 三角剖分，首先建立一个超级三角形PQR，这个三角形要把所有4点都包含进去。

接着分析A点，因为A点在三角形PQR的外接圆内部，所以利用A点将三角形PQR分拆成三个子三角形。

再考虑Ｂ点，它只在三角形AQR的外接圆的内部（而不在APQ、APR的外接圆的内部），再将三角形AQR分拆成三个子三角形。

接着是C点，此时我们已经有5个三角形，对这5个三角形每一个检查C点在不在它的外接圆内。经过检测，发现它在三角形APR和三角形ABR的外接圆内。

所以删除三角形APR和三角形ABR的公共边AR，得到四边形ABPR，并将C点与每一个边组成一个三角形。

最后对D点进行分析，它在三角形ABC和三角形BCR的外接圆内，所以应删除公共边BC，再用D点与这两个三角形的其他边形成子三角形。

最后把含有超级三角形的顶点的三角形全部删除，就得到这四个点的三角剖分

如果用一张简短的图表示上述算法中加入一个新的点的核心过程的话，那就是

但是有一个特例就是如果参与构建三角剖分的如果仅仅是三个点的话，则可能最后一步删除炒鸡三角形的时候会出不来结果，例如最后一步的时候是下图这样子的

如果删除超级三角形相关顶点就把所有三角形都删掉了，解决办法是对于这种特殊情况简单的连接三个顶点成为三角形就行了. 对于四个点之后就不会遇到这种情形了. 上面描述的过程其实就是论文【1】中Bowyer-Watson算法的伪代码.

但是论文中也提及了一个比较好的优化（能将算法的复杂度优化到 ）, 那就是事先将点集按照一个维度坐标进行排序. 如果点集的 x 坐标变化范围大的话，则选择 x 这个维度进行排序，否则选择 y 这个维度进行排序. 其实还有一个比较聪明的优化. 就是利用已经排好的序，可以不用遍历整个三角形列表. 因为对于已经确定的 Delaunay 三角形，对于新加入的点是不需要去遍历该三角形的. 因为你是按照排好的序曲遍历这些点的嘛~

我们的优化算法涉及到

  • 一个数组ps：用于存放顶点集
  • 一根链表temp_triangle_list：用于维护可能要分裂的三角形们
  • 一根链表triangle_list：Delaunay 三角形，表示已经确定了的最后一定在 三角剖分中的三角网格.
  • edge_buffer：边缓冲区, 用于存储遍历 temp_triangle_list 的过程中发现要被分裂的三角形的三条边

优化后的算法伪代码如下

struct Point
{
    double x, y;
} ps[maxn];

struct Triangle
{
    Point a, b, c;
};

list<Triangle> triangle_list, temp_triangle_list;

struct Edge
{
    Point a, b;
} edge_buffer[maxn];

struct Circle
{
    Point o;
    double r;
};

int tot; 

ilv delaunay() 
{
    sort(ps); 
    Triangle superTriangle = generateSuperTriangle();
    temp_triangle_list.push(superTriangle);
    for (re i = 0; i < n; i++) 
    {
        tot = 0; 
        for (Triangle triangle in temp_triangle_list)
        {
            Circle outter = getOutterCircle(triangle); 
            if (如果ps[i] 在 outter 右侧) 
            {
                triangle_list.push(triangle); 
            }
            else if (如果ps[i] 在outter 外部) 
            {
                continue; 
            }
            else 
            {
                则triangle一定不是delaunay三角形, 应该将此三角形从链表temp_triangle_list 中移除, 然后拆出来的三条边放进edge_buffer中,因为涉及到移除, 所以才选择链表这种数据结构;
            }
        }
        edge_buffer中去除出现次数>1的边; 
        for (Edge edge in edge_buffer)
        {
            temp_triangle_list.push(new Triangle(edge.a, edge.b, ps[i]));
        }
    }
    将temp_triangle_list并入triangle_list中;
    for (Triangle triangle in triangle_list)
    {
        if (triangle 包含炒鸡三角形的三点中的任何一点)
        {
            triangle_list.remove(triangle);
        }
    }
}

和论文【1】的Bowyer-Watson算法伪代码相比，这里最重要的优化是

我们只需要遍历 temp_triangle_list 中的三角形，而不需要每次遍历 triangle_lisit + temp_triangle_list 中的全部三角形. 而能够这样做的前提在于我们事先排了序.

有了delaunay伪代码做铺垫，再写出 Voronoi 图的伪代码就很简单了.

struct Poly
{
    vector<Point> ps;
} poly[maxn];

int poly_tot;

ilv voronoi ()
{
    for (re i = 0; i < n; i++)
    {
        getTriangles(ps[i]);
        则这一圈的delaunay三角形的外心就构成voronoi图的一个多边形, 将其放入 poly 数组中去.
    }
}

纵观上面的过程，显然我们需要写一个计算三角形外接圆的函数. 以及如何产生炒鸡三角形. 产生炒鸡三角形的思路非常简单，一图胜千言

如上图所示，我们获取点集的 x 坐标的范围 [xmin, xmax] 和 y 坐标的范围 [ymin, ymax]， 则我们可以画出一个矩形，该矩形的四个顶点分别是 (xmin, ymin)、(xmin, ymax)、(xmax, ymin)、(xmax, ymax).

那么炒鸡三角形的三个顶点的坐标显然是

P((xmin+xmax)/2, ymin-(ymax-ymin))， Q(xmin-(xmax-xmin)/2,ymax)、R(xmax+(xmax-xmin)/2, ymax)

所以得到炒鸡三角形的算法是 O(1) 的，我们很满意. 但是实际调试的过程发现一个坑, 就是其实不能这样定炒鸡三角形，而要稍微大一圈才行，不然会遇到问题. 一图胜千言.

也就是我们实际上需要的炒鸡三角形是 P1Q1R1, 而不是 PQR. 为什么呢? 实际程序跑的时候如果你选取 PQR 的话会有问题的. 还是画图.

如果炒鸡三角形是 PQR， 四个点是 A、B、C、D, 则考虑加入 A、D 之后（因为A 和D的 x 坐标小于 B、C 所以肯定先加入的是A、D两个点），temp_triangle_list 中只有 ARP、AQD、ADC 三个三角形. 然后再加入B点之后，因为 B 在这三个三角形的外接圆内部（可能AQD不太像哈，可能是画的比较烂），则 temp_triangle_list 变成 PBQ、BQR、BRP 三个三角形. 注意，这三个三角形都是包含炒鸡三角形的顶点的. 所以最后考虑 C 点的加入时（C的x坐标最大，所以最后被考虑加入），C 会在 BQR、PBR 这两个三角形内部，则最后你会发现，所有temp_triangle_list 和triangle_list 中的三角形都是包含炒鸡三角形的顶点的，则最后所有的三角形都会被移除掉的. 这就坏事儿了. 稍微分析一下上面的过程就会知道，坏事儿的关键在于 B 在 AQD 的外接圆中. 这完全就是因为炒鸡三角形不够大导致的.

最后，为了显示效果，我基于 Windows GDI 实现了基本的图形界面.

代码在 vs 2019 下调试通过. gitee 地址 https://gitee.com/yfscfs/delaunay.git

启动项目之后

输入 10，点击确定之后， 就生成了随机生成的 10个点的 delaunay 三角剖分

然后点击 转Voronni图按钮

再点击 转Delaunay图 按钮之后又会回到上一张图. 点击重置，又会回到项目启动的样子

随机文章
本站知识点必读
最近更新