Voronoi图及其应用简介
本文是博主撰写的学校课程论文(有删改),主题是几何处理领域的常用数据结构 Voronoi 图。
在计算几何领域,Voronoi 图是最有代表性的区域划分数据结构之一,在诸多领域如计算机图形学、物理化学、分子生物学乃至社会科学中均有广泛应用。在不同科学领域,Voronoi 图被发现和重视已有百年历史,以至于随发现者的不同而被赋予了不同的名称,如 Thiessen 多边形、Dirichrit 网络、plesiohedra 等。
作为基础的数据结构,Voronoi 图与 Delaunay 三角化之间的密切联系以及构建 Voronoi 图的基本算法已经得到充分研究。随着时代的发展,这一数据结构的重要性并未丧失,各领域对 Voronoi 图的应用仍可见于近年来的若干前沿研究。本文将给出 Voronoi 图与 Delaunay 图的基本概念和构建算法,并介绍在算法和数据结构和计算机图形学领域中的应用,重点强调计算机图形学领域最前沿(state-of-the-art)的应用成果,以期启发这一领域的后续研究。
预备知识
Voronoi图
在 \(d\) 维实空间 \(\mathbb R^d\) 中,给定点集 \(S = \{\mathbf p_1, \dots, \mathbf p_n\}\),称 \(S\) 中的点为Voronoi图的生成点。对每一生成点 \(\mathbf p_i\),定义
\[ R(\mathbf p_i) = \{ \mathbf p \in \mathbb R^d : || \mathbf p - \mathbf p_i || < || \mathbf p - \mathbf p_j ||, \forall j \neq i \} \]
为 \(\mathbf p_i\) 的Voronoi区域,其中 \(||\mathbf p - \mathbf q||\) 衡量了两点 \(\mathbf p, \mathbf q\) 之间的距离,我们在此采用常见的欧氏距离。Voronoi区域的几何意义相当明显:该区域恰好包含了所有到 \(\mathbf p_i\) 的距离近于到其他 \(\mathbf p_j \ (j \neq i)\) 的距离的点。
注意到,当上式的距离表达式取等时,点 \(\mathbf p\) 恰好处于到 \(\mathbf p_i, \mathbf p_j\) 距离相等的一个 \(d-1\) 维的超平面上。定义
\[ R^=(\mathbf p_i) = \{ \mathbf p \in \mathbb R^d : || \mathbf p - \mathbf p_i || \leq || \mathbf p - \mathbf p_j ||, \forall j \neq i \} \]
为 \(\mathbf p_i\) 的含边界的 Voronoi 区域,以方便我们未来的讨论。
注意到 \(\mathbb R^d = \bigcup_{i=1}^n R^=(\mathbf p_i)\),即整个空间可以被分解为各点的 Voronoi 区域和它们的边界。一个二维的欧氏距离 Voronoi 图的例子如上图所示,我们称图中各Voronoi区域边界为 Voronoi边,图中Voronoi边的相交点为 Voronoi点。一般地,考虑这一空间划分得到的各 Voronoi 区域为面,Voronoi 边为边,Voronoi 点为顶点,由此我们得到了 \(S\) 的 Voronoi图,记作 \(\text{VD}(S)\)。
Delaunay图
给定一个二维 Voronoi 图 \(\text{VD}(S)\),对每对生成点 \((\mathbf p_i, \mathbf p_j)\),我们用线段连接之当且仅当 \(R(\mathbf p_i)\) 与 \(R(\mathbf p_j)\) 相邻,即恰好有一条Voronoi 边将它们相隔。对二维情形,可以证明,该 Voronoi 边在两点连线的中垂线上。如此产生的图如上图所示,视新添加的线段为边(称之为Delaunay边),所有生成点为顶点,那么我们称新的图为Delaunay图。Delaunay 图亦可拓展到高维情况,此处不再赘述。
值得注意的是,如果对给定的一个二维 Delaunay 图和一个点集 \(S\),我们容易通过完全类似的手段还原 Voronoi 图:只需对 \(S\) 中的每对 Voronoi点,用线段连接之当且仅当恰好有一条 Delaunay 边将它们相隔。可以证明,在二维情况下,Voronoi点和生成点的数目规模都是 \(O(n)\) 的,因此,Delaunay 图及其对应 Voronoi 图之间可以用 \(O(n)\) 的时间互相还原。在数学上,Voronoi 图及相应的 Delaunay 图是互为对偶图的关系,指的就是这一性质。
Delaunay图与Delaunay三角化密不可分------一般情况下,Voronoi图的生成点集 \(S\) 的Delaunay三角化结果正是该图对应的Delaunay图,但对于退化情况则需要特殊讨论。如上图所示,一般地,一个 \(d\) 维Voronoi图中的Voronoi点与 \(d+1\) 条Voronoi边相连。而如果该数目超过 \(d+1\),则称这一点为退化点,该Voronoi图出现了退化情况。
可见,退化的二维Voronoi图对应Delaunay图产生了边数大于 \(3\) 的多边形,但处理方法也很简单------我们只需依次连接该多边形的对角线(无所谓连接哪两个顶点),将其分割成若干三角形即可,如上图所示。因此,可以认为Voronoi图的退化情况无伤大雅,Voronoi图与Delaunay三角化有相当密切的一对一关系。
Delaunay三角化是计算机图形学的基础算法之一,能够产生高质量的三角化结果,具有许多优秀的性质,如最大空外接圆,最大化最小角等,被频繁地应用于大量的图形学算法之中。可见,对偶Delaunay图的应用亦十分广泛,鉴于本文的重点是Voronoi图,故此处不再赘述。
构建Voronoi图的算法
本小节将介绍从点集 \(S = {\mathbf p_1, \dots, \mathbf p_n}\) 构建Voronoi图的算法。由于Voronoi图和Delaunay图的对偶关系,两者的构建是等价的,因此所有构建Delaunay图(Delaunay三角化)的算法均可用于构建Voronoi图。据我们所知,大部分的算法构建的是Delaunay图,本小节介绍的也均是构建Delaunay图的算法。为使算法易于描述,我们省略高维情形的拓展,只讨论二维情形的算法。
圆扩张算法
圆扩张算法的思想来源于Delaunay三角化满足的最大空外接圆性质:(非退化情况下)任意Delaunay三角形的外接圆不会包含其他生成点。注意到任意三角形的外接圆唯一,而线段的外接圆不唯一,于是,圆扩张算法可以如下描述:
- 从任意未进行过圆扩张的生成点 \(\mathbf p_i \in S\) 出发,寻找其最近邻点 \(\mathbf p_j\),如上图(a)(b)。
- 确定 \(\mathbf p_i\) 与 \(\mathbf p_j\) 为一条Delaunay边,如上图(c)。
- 若该Delaunay边不是 \(S\) 的凸包的一条边,则向两侧扩张该边的外接圆,分别直到找到一个三角形 \(\mathbf p_i, \mathbf p_j, \mathbf p_k\) 的外接圆;若该Delaunay边是 \(S\) 的凸包的一条边,则只向 \(S\) 的凸包内侧扩张该边的外接圆,直到找到一个三角形 \(\mathbf p_i, \mathbf p_j, \mathbf p_k\) 的外接圆,如上图(d)(e)。
- 确定三角形 \(\mathbf p_i, \mathbf p_j, \mathbf p_k\) 为一个Delaunay三角形,如上图3。回到步骤(1)。
存在平均时间复杂度为 \(O(n \log n)\) 的计算二维凸包的算法;对每个生成点 \(\mathbf p_i\),该算法均执行一次圆扩张,而查询最近点、扩张边外接圆的时间复杂度均为 \(O(n)\),故圆扩张算法的整体时间复杂度为 \(O(n^2)\)。该算法属于朴素算法,思路与求最小生成树的Prim算法(加点法)相当类似,有很大的优化空间,实际工程中不会采用该算法。
抛物面提升算法
抛物面提升算法为将\(d\)维的Delaunay图问题归约到了\(d+1\)维的凸包问题,思想相当巧妙。算法基于接下来给出的一个重要引理。
引理 考虑二维平面上的点\(\mathbf p_i, \mathbf p_j, \mathbf p_k, \mathbf p\),记抛物面提升映射\(f: \mathbf p \in \mathbb R^2 \to \mathbf p' \in \mathbb R^3\),满足\(f(x,y) = (x,y,x^2+y^2)\),如上图所示。记\(\mathbf p_i', \mathbf p_j', \mathbf p_k', \mathbf p'\)分别为这四点经过映射的点。那么,\(\mathbf p\)在\(\mathbf p_i, \mathbf p_j, \mathbf p_k\)确定的圆内/外,当且仅当\(\mathbf p'\)在\(\mathbf p_i', \mathbf p_j', \mathbf p_k'\)确定的平面的下半/上半空间。
引理的证明可参见附录。根据这一引理,我们容易推导出二维Delaunay图问题与三维凸包问题之间的关系。记\(S'\)为点集\(S\)经过抛物面提升映射得到的三维点集。为了满足最大空外接圆性质,\(S\)中的\(\mathbf p_i, \mathbf p_j, \mathbf p_k\)是Delaunay三角形,当且仅当其提升后得到的\(S'\)三角面\(\mathbf p_i', \mathbf p_j', \mathbf p_k'\)恰好是\(S'\)的凸包的一个低面,如此才能保证三角面\(\mathbf p_i', \mathbf p_j', \mathbf p_k'\)的下半空间不包含其他\(S'\)中的点。如上图所示,\(S'\)的凸包的任意低面(lower face)的下方不再包含\(S'\)的其他点,将所有低面投影至\(z = 0\)的\(\mathbb R^2\)平面,即可得到原\(S\)点集的Delaunay图。于是,抛物面提升算法的流程可以如下描述:
- 将\(S\)中的所有点进行抛物面提升映射,得到三维点集\(S'\)。
- 求解\(S'\)的凸包。
- 将\(S'\)的凸包投影至\(S\)所在平面(步骤(1)的逆映射),得到Delaunay图。
一般地,\(d\)维的Delaunay图问题均可如此归约到\(d+1\)维的凸包问题。\(d=2\)时,由于存在平均时间复杂度为\(O(n \log n)\)的计算三维凸包的算法,而映射每个点仅需\(O(n)\)时间,故抛物面提升算法的整体时间复杂度为\(O(n \log n)\),这也是已知的求解二维Delaunay图的算法的时间复杂度的下确界。
其他算法
求解Delaunay图的算法相当丰富,除了时间复杂度之外,数值稳定性也是重要考量之一。由于抛物面提升算法需要计算\(x^2 + y^2\)这一较大的浮点数值,导致一定的精度误差,实际工程中一般不会采用该算法。其他常用的算法包括Divide and Conquer、Fortune's Algorithm等,此处不再赘述。
距离扩展的Voronoi图
通常的Voronoi图的定义采用欧氏距离,可以用其他形式的距离替换之,以满足应用场景的需要。一般地,记\(h(\mathbf p, \mathbf p_i)\)为两点\(\mathbf p, \mathbf p_i\)的广义距离函数,那么在此定义下的Voronoi区域为
\[ R_h(\mathbf p_i) = \bigcap_{\mathbf p_j \in S, j \neq i} \{ \mathbf p \in \mathbb R^d: h(\mathbf p, \mathbf p_i) < h(\mathbf p, \mathbf p_j) \}. \]
一般地,构建距离扩展的Voronoi图的算法与欧氏距离定义下的Voronoi图存在显著不同。对某些特定的距离函数,构建Voronoi图的算法远不及欧氏距离下的经典情形成熟,仍存在较大的研究空间。限于篇幅,本小节仅讨论应用较广的距离函数:球面距离、加法加权距离、power distance、乘法加权距离,以及一般的\(L_p\)距离。
球面距离
若点集\(S\)分布在一个球面上,则Voronoi图的距离函数应当是球面距离函数,如Haversine公式。如此生成的Voronoi图也位于球面上,可用于地理信息系统的建设。例如,在地球上建立Voronoi图时自然应当应用球面距离函数,如上图所示。
加法加权距离
\[ h(\mathbf p, \mathbf p_i) = || \mathbf p - \mathbf p_i ||_2 - w_i, \]
其中\(||\cdot||_2\)表示欧氏距离范数。这就是加法加权距离(additively weighted distance)。使用这种距离函数定义的Voronoi图如上图所示。
使用加法加权距离函数,可以量化各生成点对其Voronoi区域的控制程度。举例说明,当\(\mathbf p\)与\(\mathbf p_i\)的距离小于\(w_i\)时,\(h(\mathbf p, \mathbf p_i) < 0\),若我们对距离函数取绝对值处理,即记\(h(\mathbf p, \mathbf p_i) = \text{abs}(|| \mathbf p - \mathbf p_i ||_2 - w_i)\),这时可以保证\(\mathbf p\)被\(\mathbf p_i\)"绝对控制",即\(\mathbf p\)一定在\(R^=(\mathbf p_i)\)中。加法加权距离下的Voronoi图具有较广泛的应用,如量化不同规模的购物中心对附近居民的吸引程度等。
Power distance
Power distance(又称Laguerre distance)定义如下:
\[ h(\mathbf p, \mathbf p_i) = || \mathbf p - \mathbf p_i ||^2_2 - w_i, \]
称使用这一距离函数的Voronoi图为power diagram。可见,power distance与加法加权距离的区别仅在于power distance将距离取了平方,但这一变动使得power diagram具有许多巧妙的性质。power distance在二维情形下的定义可以由勾股定理解释:考虑在点\(\mathbf p_i\)处作半径为\(\sqrt w_i\)的圆,那么在圆外的任意点\(\mathbf p\)到\(\mathbf p_i\)的power distance的平方恰好是从\(\mathbf p\)向该圆引的切线段的长度。由圆幂定理,容易推导两个生成点对应圆之间的交点所在直线恰好包含Voronoi边,如上图所示。power distance下的Voronoi图可用于高效计算空间中三维球的并集的体积,检测点是否在圆盘的并集内等,应用广泛。
乘法加权距离
与加法加权距离类似,乘法加权距离的定义式为
\[ h(\mathbf p, \mathbf p_i) = \frac{1}{w_i}|| \mathbf p - \mathbf p_i ||_2. \]
使用这种距离函数定义的Voronoi图如上图所示。
乘法加权距离则量化了各生成点对其Voronoi区域的扩张程度,\(w_i\)越大,扩张程度越高。我们取极限说明:当\(w_i \to 0\)时,距离函数\(h(\mathbf p, \mathbf p_i) \to +\infty\),该生成点的Voronoi区域的面积趋向于\(0\);当\(w_i\)充分大时,距离函数\(h(\mathbf p, \mathbf p_i) \to 0\),有利于该生成点扩张其Voronoi区域的大小。乘法加权距离下的Voronoi图也具有较广泛的应用,如量化生态系统下各竞争物种扩张能力等。
\(L_p\)距离
记\(\mathbf p = (a_1, \dots, a_d)\),\(\mathbf q = (b_1, \dots, b_d)\)为\(\mathbb R^d\)中任意两点,定义\(L_p\)距离为
\[ h_p(\mathbf p, \mathbf q) = \left(\sum_{i=1}^d |a_i - b_i|^p \right)^{1/p}, \]
可见\(L_2\)距离即为欧氏距离。\(L_1\)距离又称曼哈顿距离,可衡量规整的城市道路联通的两点之间的距离,上图展示了\(L_1\)距离下的二维Voronoi图的一个例子。定义\(L_\infty\)距离为
\[ h_{\infty}(\mathbf p, \mathbf q) = \lim_{p \to \infty} h_p(\mathbf p, \mathbf q) = \max_{i=1,\dots, d} \{|a_i - b_i|\}, \]
可见\(L_\infty\)距离衡量了两点之间在某个坐标轴的投影距离上的最大差值。在二维情形,\(L_\infty\)距离的应用包括衡量在水平轴和竖直轴上做等速率匀速运动的数控机床从点\(\mathbf p\)运动点\(\mathbf q\)所需的时间。上图展示了\(L_\infty\)距离下的二维Voronoi图的一个例子。
除了\(p=1,2,\infty\)的情况,其他\(L_p\)距离也有相应的应用场景,此处不再赘述。
在算法与数据结构中的应用
Voronoi图在算法和数据结构上的应用相当广泛,降低了求解许多非图形学的计算几何问题的时间复杂度。
最邻近搜索问题
一个最简单也是最典型的例子,就是最邻近搜索(Nearest Neighbor Search)问题:给定\(\mathbb R^d\)空间中的\(n\)个定点构成的点集\(S = \{\mathbf p_1, \dots, \mathbf p_n\}\);接下来\(m\)次查询,第\(j\)次查询将询问点\(\mathbf q_j \in \mathbb R^d\)的最近邻点\(\mathbf p_i \in S\)。
显然,该问题有一个朴素的\(O(mn)\)算法:对每次查询,对所有\(n\)个点线性计算\(\mathbf q_j\)与\(\mathbf p_i\)的距离,取最小距离即可。但若考虑将Voronoi图作为数据结构,可以设计如下的算法:
- 为\(S\)建立Voronoi图(以及相应的用于查询点的位置的数据结构,如Kirkpatrick)。
- 对每次查询\(\mathbf q_j\),根据步骤(1)建立的数据结构,取满足\(\mathbf q_j \in R^=(\mathbf p_i)\)的\(\mathbf p_i\)。
对于二维情形(\(d=2\)),步骤(1)耗时\(O(n \log n)\),步骤(2)对每次查询可做到耗时\(O(\log n)\),故算法的整体时间复杂度为\(O(n \log n + m \log n)\),优于朴素方法。其他空间数据结构如k-d 树、多叉树等亦可解决这一问题,并具有相同的时间复杂度。
容易想到,最邻近搜索问题在地理信息系统中应用广泛,如查询当前坐标距离最近的公共设施的位置等。在计算机图形学中,也存在大量最邻近搜索问题,通过建立Voronoi区域划分的层次结构也是解决方案之一。
最小圆覆盖问题
详见我的另一篇文章。
其他计算几何问题
Voronoi图在计算几何建模的图像分割和路径规划等问题中均得到了应用,限于篇幅,本文不再展开介绍。
在计算机图形学中的应用
众所周知,Delaunay三角化是图形学最重要的基础算法之一,而我们已经知道Delaunay图就是Voronoi图的对偶。除此之外,Voronoi图在计算机图形学中也得到了许多直接的应用。本节将介绍在几何处理和渲染领域中应用Voronoi图的研究成果。
在几何处理领域中的应用
Voronoi图在几何处理领域的经典应用包括计算形体的中轴,以及计算偏移曲线和偏移曲面等。在网格生成与重建领域,有关Voronoi图和Delaunay图的相关工作层出不穷。经典的Tetgen算法利用 Constrained Voronoi 图及其对偶 Constrained Delaunay 三角化,在三维曲面上实现了高效的四面体网格生成算法,后续工作TetWild 等也基于相关概念实现了更稳健的四面体网格生成算法。
在脚手架实体建模的相关研究中,球面上的Voronoi图可用于平滑混合(blending)不同方向的脚手架多面体的伸展。后续工作进一步将该方法扩展到了稳健的管状实体建模领域,这是已知的三维CAD建模前沿应用Voronoi图的研究工作。
从三维Voronoi图重建网格也是一个有趣的研究课题。研究工作给出了一个从球上的Voronoi图还原网格拓扑的GPU并行算法,其效率显著高于过去的串行算法,并在任意分布的Voronoi生成点上取得了稳健的还原结果。
在有向距离场(SDF)定义的距离下,构建Voronoi图可以有效提升后续网格重建操作的效率,因此如何构建这样的预处理数据结构也具有重要研究价值。提出了在SDF下构建受限制的Voronoi图(Restricted Voronoi Diagram)的高效方法, 在此基础上应用Centroidal Voronoi Tessellation算法从SDF表达的曲面生成了高质量的三角网格。
作为基础的数据结构,优化现有的Voronoi图构建算法同样具有重大意义。Inria研究机构的工作提出了一种新颖的构建三维形体的power diagram的GPU并行算法。同样,近年来,出现了一种构建三维\(L_{\infty}\)距离下的Voronoi图的通用算法,高效地完成了对具体积的模型(volumetirc model)上的空间划分任务。可以预见,这些算法对现有利用Voronoi图的许多图形学应用而言具有重大价值。将这些新颖的高效算法应用到图形学研究和工程后,可以更高效地执行Centroidal Voronoi Tessellations等曲面细分算法,推进前沿网格生成和重建领域的研究。
在渲染领域中的应用
在渲染领域,近年来也有部分工作围绕Voronoi图展开。与传统单色光渲染任务相比,全频率渲染(all frequency rendering)任务需要考虑不同频率(或波长)的光线在模型表面上的双向反射分布函数(Bidirectional Reflectance Distribution Function, BRDF)随频率变化而产生不同的结果,进而引入了更多的挑战。在这一领域的光线可见性(或称遮挡性)判定问题上,工作利用球面上的Voronoi图,高效地判定了给定的渲染表面是否能够接收到各频率的光线,其思路来源于球面Voronoi图在各Voronoi区域上的最近距离性质。在Voronoi图对可见性判定的加速计算下,全频率渲染结果相对于传统方法取得了更高的效率,并在存在阴影和遮蔽关系的渲染中取得了更接近于ground truth的效果。
用于VR头盔的注视点渲染(foveated rendering)是虚拟现实技术的研究热点之一。在执行常规体渲染(volume rendering)后,在图像重建(image reconstruction)的后处理环节执行了基于Voronoi区域的最近邻插值,以生成平滑并带有突显注视点效果的渲染结果。
总结
本文对Voronoi图的概念、算法和应用做了综述,详细描述了Voronoi图的计算几何数学背景,讨论了相关的若干计算几何方面的算法和数据结构的问题,并介绍了前沿的图形学学术界对Voronoi图的应用成果。限于篇幅和作者能力,本文无法涵盖学术界全部的前沿Voronoi图研究工作,但所选工作均具有一定的代表性,可供研究人员参考。目前,对Voronoi图的应用研究方兴未艾,可以相信,这一经典的数据结构将持续为各方面的研究贡献更高效的解决方案。
附录
对引理的证明
引理 考虑二维平面上的点\(\mathbf p_i, \mathbf p_j, \mathbf p_k, \mathbf p\),记抛物面提升映射\(f: \mathbf p \in \mathbb R^2 \to \mathbf p' \in \mathbb R^3\),满足\(f(x,y) = (x,y,x^2+y^2)\)。记\(\mathbf p_i', \mathbf p_j', \mathbf p_k', \mathbf p'\)分别为这四点经过映射的点。那么,\(\mathbf p\)在\(\mathbf p_i, \mathbf p_j, \mathbf p_k\)确定的圆内/外,当且仅当\(\mathbf p'\)在\(\mathbf p_i', \mathbf p_j', \mathbf p_k'\)确定的平面的下半/上半空间。
证明
考虑抛物面\(z = x^2 + y^2\)上的点\((a,b,c) \in \mathbb R^3\),过该点的切平面方程为
\[ n_x(x - a) + n_y(y - b) + n_z(z - (a^2 + b^2)) = 0, \]
其中\((n_x,n_y,n_z)\)为该点处抛物面的法向量。通过计算\(f = z - x^2 - y^2\)的梯度可得到\(n_x = 2a, n_y = 2b, n_z = -1\),于是上式的平面方程化为
\[ z = 2ax + 2by - (a^2 + b^2). \]
我们将该平面向
\(z\)轴正方向平移\(\rho^2\)单位长度,得到对抛物面的割平面
\[ z = 2ax + 2by - (a^2 + b^2) + \rho^2. \]
在割平面与抛物面的交线处有\(z = x^2 + y^2\),代入\(z\)化简得到
\[ (x-a)^2 + (y-b)^2 = \rho^2, \]
这意味着割平面与抛物面的交线在\(z=0\)上的投影恰好是一个圆心为\((a, b)\),半径为\(\rho\)的圆,该圆可以由三点\(\mathbf p_i, \mathbf p_j, \mathbf p_k\)确定。注意到对任意二维点\(\mathbf p = (x,y)\),其抛物面提升映射的输出\(\mathbf p' = (x, y, x^2 + y^2)\)与割平面的竖直距离为
\[ 2ax + 2by - (a^2 + b^2) + \rho^2 - (x^2 + y^2) = \rho^2 - ( (x-a)^2 + (y-b)^2 ), \]
于是,当\((x,y)\)在圆内时,\((x-a)^2 + (y-b)^2 < \rho^2\),上式\(>0\),意味着\(\mathbf p'\)在割平面的下半空间;同理,当\((x,y)\)在圆外时,上式\(<0\),意味着\(\mathbf p'\)在割平面的上半空间。证毕。