计算几何课程
04Voronoi图(Voronoi Diagram)
本章主题介绍平面如何划分。通过对平面划分成距离点最近的区域，与Voronoi图相关联。Voronoi图是一个经典问题，同时介绍了计算几何中的一个经典数据结构DCEL。
- Introduction
- 介绍Voronoi图的来源，应用场景。有一篇很不错的背景介绍文章。
- Terminologies
- 介绍Voronoi图相关概念术语。
- Properties
- 介绍Voronoi的特性。主要是等距性和非空性。
- Complexity
- 简单介绍Voronoi的存储复杂性。
- Representation
- 介绍对Voronoi的存储，表示。主要方式为主区域刨分(Subdivision)，同时也适用于其他一些表示。同时为了快速的查找，引入了DCEL。
- DCEL
- 介绍了DECL数据结构来存储平面图结构。
- Hardness
- 规约技术界定了Voronoi图的复杂度，使用的是排序方式来规约操作。
- Sorted Sets
- 关注在某些情况下Voronoi图的求解问题是否会更容易。例如在某个方向排好顺序，类似凸包的Graham Scan方式，是否可以有效降低复杂度。
- $VD_{sorted}$
- $VD_{sorted}$即排序的Voronoi图复杂度。在这里引入了$\epsilon{}-Closeness$问题并将其归约到$VD_{sorted}$来说明排序对于降低Voronoi图的复杂度并没有帮助。值得一说的是，即便知道Voronoi图的问题也不能帮助排序，看似两个问题结构毫无关联。
- Naive Construction
- 简单说明一下最朴素的构造方式。
- Incremental Construction
- 增量式构造。即每次增加一个点的方式来构造Voronoi图
- Divide-And-Conquer
- 基于分而治之的思路，将平面划分，分成多个区域然后合并起来操作。
- Plane-Sweep
- 基于扫描线的方式来生成Voronoi图结构。基本想法是，根据已经扫描过的点来确定一些潜在的生成Voronoi点的事件位置。高层想法是基于抛物线结构。
Voronoi图介绍
Voronoi图第一次由G. F. Voronoi(1868 - 1908)第一次正式的提出并且界定，于是后来这种图结构，一般都称之为Voronoi图。直观来看Voronoi图由一系列点和边构成，这些边会把点隔开，各自围成一个区域，在该区域内的点距离该点的距离是最近的，如下图：
- 这些中心点，我们称之为Voronoi Site，简称Site。
- 这些区域，我们称之为Voronoi Cell，简成为Cell。从几何上来看，该Site之所以拥有这个Cell，是因为这个Cell中的点距离这个Site是最近的。
用数学的方式定义来说如下：
Voronoi图由一些d维空间中点$S=\Set{p_1,\cdots,p_n}$构成，每个点$p_i$称之为Site。对于每一个Site来说，存在一个区域Cell，该Cell包围目标$p_i$，该区域可表示为$Cell(p_i)=\set{q\in{}\mathbb{R}^n|d(q,p_i) 根据这些定义我们可以得到一些特性： 每个Cell都是凸的。 因为对于每个Cell可以通过$p_i$跟其余$n-1$个点求所在半平面的交来来确定，因为半平面可以认为是一个凸多边形。而凸多边形在交运算下保持，所以Cell也是凸的。 只当所有点都在一条线上的时候会出现直线。因为存在直线，那么所有点间平分线都必须平行。有一个相交则都有交点。所以所有Site都在一条直线上。所以直线数目最哦多不超过$n-1$条。 Non-empty Cells 考察每个Site中，以$p_i$为中心的一个半径为$s$的$Disk(s,p_i)$。可以看到这个Disk必然会落到自己的Cell中。这是因为这里点集是一个有限集合，所以每个点之间的距离一定有一个最小值，取$s$为最小值的一半，这个$Disk(s,p_i)$里的点必然离这个$p_i$最近，所以会落到在Cell中。 Empty Disks 这个地方可以简单讨论一下。当$x$落在Cell中的时候。可以发现，这个Disk的最大半径就是$x$到Cell所对应Site的距离。当$x$落在Voronoi Edge上时，其对应的最大半径就是$x$到该边相邻两个Cell上对应Site的距离。这两个距离时相等的。当$x$落在Voronoi Vertex上时，同理。这些都可以通过对Voronoi图的定义来得到。 Nearest = Concyclic 这里并没有强调，这个最近邻Site有几个，就算有多个，由上面讨论，加上圆的定义可知，这个是显然的。实际上也可以反过来说，如果有两个最近邻Site，那可以认为$x$落在一条边上。如果有不少于3个，那么这个$x$一定是一个Voronoi Vertex。 注意这里没有说一定要是最近邻的点了。所以实际上取一个圆上4个点，使其构成的四边形不相互平行就可以了。 最后我们可以考察一个点$x$在边上的运动。当一个点在边上运动的时候，考察其最大空圆，可以看到这个空圆伴随变化。等到了Voronoi Vertex的时候，就会看到这个最大空圆犹如发生了分裂，沿着其余各个线的方向传递过去。 下面就是考虑Voronoi图的存储描述问题，如果其存储过大，那显然其结构不会简单，时间复杂度不会低。所以这个时候就需要对其规模进行一个界定。其实方法就是从平面图出发，可以注意到Voronoi图一定是一个平面图结构，就可以掏出我们无敌的欧拉公式来处理他，处理方式可以如下图。 可以添加一个无穷远点，并且将所有射线连接到这个无穷远点。这样就会成为一个标准的平面连通图。对于欧拉公式来说我们有：
其中$v$是顶点个数，$e$是边个数，$f$是面个数，而$c$则是连通分支个数。通常来说这个$c=1$，即只有一个连通分支。也是熟知连通平面图的欧拉公式。这里面我们有$n$个Site，所以$f=n$。 此时我们还可以注意到对于Voronoi图来说，每个Vertex其度数一定大于等于3，即$3v 所以我们可以得出上界 可以看到对于$n$个Site的Voronoi图产生的Vertex和Edge都是线性的。 DCEL结构 接下来就是讲如何存储Voronoi图，这里用的是一种通用的结构，子区域刨分结构。该结构可以用来存储所有类似的平面结构，即讲一个平面细分成一个一个的区域结构。而这些区域互相没有重叠，同时没有任何覆盖。Voronoi图则相当于这种形式的一个特例。 现在来看这样一个数据结构的设计，首先有个HighLevel思想，设计需要达成什么目标。这里通俗来讲需要： 而DCEL结构切可以满足上述要求。DCEL是一种基于边的数据结构，结构大体来说由点，边，面三部分来组成。其最主要的技术就是边分解技术，即将一条有向边，分解成两个有向边，一个正向，一个逆向。而形成的数据结构可由下面三类表组成： Half-Edge Half-Edge数据结构如图。其相当于一条边中的那半条有向边。作为一条边的两部分，其要记录其孪生兄弟边(twin)。然后作为有向边要记录其起点(ori)。同时因为有向，可以定义其左右方向，就可以区分其两侧面的位置。这里约定以有向边的左侧定义其关联面，称之为Incidence Face，记录为inc。因为然后为了遍历需求，其还要记录其前后的Hals-Edge。记录在pred和succ。 注意的是有向边的终点可以通过其twin的ori来确定。同时这个前驱和后继是唯一的，这个唯一性可以由其Incidence Face来确定。 Vertex Vertex数据结构如图。对于一个顶点来说，要记录其坐标，即x和y。同时要记录这个顶点上关联的边。这里的记录方式是借助Half-Edge，通过类似于链表的方式来实现。准确来说，其随机记录一个出去的Incidence Edge来作为找到围绕顶点遍历的第一条边。随后通过后面的方式来遍历。 Face Face数据结构如图。Face结构更少，其存储的只有一个inc字段。这个存储着一条Hals-Edge 简述几个遍历操作思路 通过f的inc找到其边上的一个Half-Edge，然后通过其前驱后继遍历即可。 通过v的inc找到第一条Half-Edge，我们记录为e。然后注意到e->pred这条有向边，这个有向边，以v为终点。这相当于找到了另一条跟v关联的边。然后现在关注其孪生边e->pred->twin。这个边将是一个以v为源点的边，同时其所属在另一个平面中，这说明其pred边将不会是原先的e。这个可以通过简单画图看出来。所以可以通过这样的方式不断翻越，遍历v所关联的所有边。 现在考虑给定点集合，构造出上述DCEL的时间复杂度。那么我们首先需要界定一下这个问题大致的难度。 这里直接的想法就是规约。找到一个方法来界定下限。这里用到的规约方法就是如下的： 其中$1D-VD$即1维的Voronoi图问题，对于1维来说，这即给定一系列点，确定其每个点间的最近距离平分线位置，所以规约如下进行。 现在剩下到2维的规约。这里直接给出了一个$Sorting\leq_N{}2D-VD$的规约。 其思路就是将排序点集放在圆上，然后通过2维Voronoi图算法来处理。然后此时因为所有点都在圆上，所以此时Voronoi图只有唯一的一个Vertex即那个圆的圆心。并且每个圆上每两个相邻点有一个平分射线，有Vertex发出设想外部。所以规约可以如下进行： 这里还提到了一种跟泛用的规约，即用抛物线映射，而不是圆来提到一样会得到相同的结论。 现在考虑一些特殊结构情况，即拥有一些特殊结构后，其复杂度是否会降低？ 例如给出一个排序集合是否可以降低复杂度？在凸包问题中这个是可以降低到$O(n)$量级的。但是很可惜的是在Voronoi图问题中，即便已经排好序其任然要花费$\Omega(nlogn)$时间。 这个问题实际上说明了凸包和Voronoi图结构本质上的区别。 区分这两个方式就是通过仿射变换来看(affine transform)。在仿射变化下，凸包具有不变性，即原来构成凸包的点，依旧是变化后点集构成凸包的点。但是对于Voronoi图来说，则不是，例如一个压缩变化，那么Voronoi则会发生巨大改变。实际上Voronoi图来说，强烈依赖于距离这个概念，距离这个概念则是几何中相当重要的一部分。 $\mathbf{VD_{sorted}}$问题 现在来说明，排好序的结构，对于Voronoi图来说，依然是$\Omega(nlogn)$。这个问题归结于$\epsilon-Closeness$问题，对于这个问题是Min-Gap的判定问题，其时间复杂度依然是$\Theta(nlog)$的。剩下就是看规约构造这个排序问题。 这个规约可以用如下图表示： 根据规约我们应该是 所以首先是一步骤提升变化对于$\epsilon-Closeness$的输入$x_1,\cdots,x_n$我们构造如下点$(x_i,\frac{i}{n}\epsilon)$作为Voronoi图的输入。这一步类似于提升操作，把所有点升到2维中去。 然后我们通过Voronoi图算法给出Voronoi的DCEL结构。此时我们关注穿过x轴的那些Face和Edge。这个是很容易做到的，因为一定有一个负无穷远点。判断这个点在什么Site里面即可，然后通过Site找到Face不断遍历边判断是否跟X轴相交，即可找到x轴穿过的所有Site。 现在我们关注每个点$x$是否都落在自己提升后点$(x_i,\frac{i}{n}\epsilon)$所张成的Cell中。这两有两种情况 那么存在一个提升后$I$，其对应的原来点$i$落在了$J$张成的Cell中。那么如图所示，我们可以简单得到关系: 最简单来讲，我们可以通过原始暴力算法来做。方法思路就是，每个点原始占用整个平面，然后通过计算每个点跟其他点的平分线来剪切整个平面，最后依次计算出每个Voronoi图Cell。显然这个方法非常笨重，甚至无法有效构造DCEL结构。所以后面来看更高效的思路。 增量方法构造。 思路就是一开始相当于只有一个Site点，此时DCEL即全平面，然后按照顺序逐个去添加Site点。每添加一个Site点都会跟原来已经存在的Site点进行竞争，从原来某些Cell中划去部分面积。 所以核心问题就是，当平面上已经有了$k$个点，如何添加第$k+1$个点，同时操作原来Voronoi图的DCEL结构来适配新的Voronoi图。这步操作如下图： 算法流程如下： 之所以通过交点翻越边的方式可行。可以认为这是因为对于新加入Site点的Cell，其边是一个连续的凸包。所以对于每一个Cell边界上的点，都会有两条平分线构成的边经过。我们只要按顺序找出来即可。 现在来看这个算法的复杂度是多少。根据步骤我们一个个去查看。 所以每次引入一个点其时间复杂度是$O(nlogn)$量级，所以整个算法时间度是$O(n^2logn)$。 这里可能有一些疑惑，其实上面这个算法是不是不一定能达到最坏情况。毕竟平分线与凸多边形求交部分这个凸多边形的边数不一定与$n$持平，因为我们知道整个Voronoi图的所有边不过$O(3n-6)$条，新加入点的Cell边界也不会有$O(n)$量级。 不过可惜的是这种最坏例子确实存在，老师也举例了出来，就是形如抛物线上的一系列点$(x_i,x_i^2)$。 分而治之构造。 思路就是一开始把点切分成两部分，然后两部分都求解出对应的Voronoi图结构。然后剩下就是把两个Voronoi图结构合并起来。 思路很简单，犹如归并排序。但关键点其实就在这个怎么合并上面。 合并的大致方式可通过下图来表示： 我么设左右点集为$S_L,S_R$对应构的Voronoi图为$VD(S_L),VD(S_R)$。 在Merge开始前，左右两个Voronoi图都会铺满整个平面。这时其边界上的Cell一定会出现某些冲突。例如图中表示的一样，左边的Cell(l)和右边的Cell(r)有部分重叠在一起。这样就需要重新划分边界，使得两者不再重叠。方法就是找到两个Site点的平分线(bisector)来重新切割两个Cell，然后逐个去解决每一对冲突的Cell。这思路是显然的，不过如果不仔细处理的化，这个方式将会是低效的。 这里其实有个关键问题，就是需要去操作的点有哪些？如果盲目的去找，那么其实有大量的无用比较。这里有一个直接观察就是，合并的裂缝只在边界上，且只影响很小的一部分区域。而且，可以看到这条裂缝是连续的，可以设想，从一个点出发连续的去探索下一个点的方式是可行的。 我们称这条缝合线线为轮廓线(Contour)，我们可以看到这条轮廓线的一些特性： 为什么不会出现转弯 老师给的提示已经非常正式化了。可以定义轮廓线分别到左右两侧Site点距离带有方向符号。类比向量，例如向X轴正向为正，向负向为负。如果出现转弯，左右两侧距离会出现正负变化。这说明左右凸包中的点在X轴上次序发生改变。而这与之前的Sorting冲突。 所以现在思路就是找到这条轮廓线的第一条。可以说明的是，这条线一定是，左右点集凸包$conv(S_L),conv(S_R)$上下公切线(Common Supporting Lines)的平分线(Bisect)。因为无穷远点也会在轮廓线(Contour)上，所以其到最近邻两侧Site的距离是一样的。而这两侧的Site实际就是$conv(S_L),conv(S_R)$的公切线两端点。 具体构造方式 现在来看整个具体的构造过程。假设此时左右Site点集合分别为$S_L,S_R$。对应构造的Voronoi图为$VD(S_L),VD(S_R)$。整个过程可以如下描述： 复杂度分析 在上面步骤中有一个关键细节。如果只是按照传统的线与多边形求交计算，那么可能需要遍历Cell的所有边界，其实根据上面直线跟凸包相交至少也要$O(logn)$时间。而且存在情况其会反复扫描求交其中某一侧的Cell。如下图： 这有可能导致一次扫描求交是$\Omega(n*m)$量级的。而这有个方式而已规避这个反复扫描过程。将累计扫描求交限制在不超过这个Cell大小量级。 其观察就是这个轮廓线一定是某个Cell的边界，而某个Cell边界是凸的。凸的就说明在这个边界上行进是都沿着统一个方向转弯的。对应而来就可以发现轮廓线跟原来凸包边界上的交点，也是按照统一个方向前进的。例如途中黄色1，2，3点。都是按照一个方向去前进的，这说明，这些交点分布是单调的，进而每次求交都排除了一部分不用再计算的边界，进而控制住时间复杂度。 由此我可以得出整体算法的时间复杂度。对于一个左右大小为$n,m$的凸包来说： 也就是说整体来说，这个合并过程不超过$O(n+m)$量级，进而算法整体构造时间复杂度不超过$O(nlogn)$量级。 扫描线方式。 当用扫描线构造Voronoi图算法首先会碰到一个难题，就是无法预测Voronoi图边的构造位置，必然存在情况，要回退去处理之前某些边的终止位置情况。 而扫描线算法大致的思路基于一件事情，即到目标Site点和到目标直线的距离关系。如果到点和到线的距离相等就是一个抛物线，这个抛物线划分了整个区域，确定了距离关系。 这个算法构造非常精妙，基于原来Fortune的思路，进行了简化，构造了这抛物线的扫描线。但是这抛物线扫描线只是为了直观理解算法的运行方式而构建的几何直观。最后再落到细节实现上。 这个奇怪的扫描下可以用下图表示： 相对于原来的扫描线，其划分成了两部分。一条直线，以及跟在其之后的由无数抛物线构成的曲线，我们称之为海岸线(BeachLine)。 可以看出来，这样扫描线一共把区域划分成了三部分。 Fronzen区域的Voronoi图是固定，其不会再受到扫描线之前那些未处理Site点的影响。这是因为对于这个区域的一个点$x$来说，其到扫描线的距离，大于到区域内任意Site点的距离。这说明，这个点$x$一定是区域内某个Cell内的点，不会归属于未扫描到的点。 所以现在可以看到，如果我么可以确定这个BeachLine。我么就可以有效的确定，未扫描到的点会影响的区域，并且对于Fronzen区域我么可以有效构造Voronoi图。而我们前面知道BeachLine是由无数条抛物线构成的最低包络线。现在就要看怎么通过无数抛物线构造BeachLine，如下图： 现在我么可以直接的感觉到，这个抛物线的相交点是一个关键点。恰恰是这些交点勾勒出了整个Voronoi图形状。直观来看： 这个过程可以通过下面几个分析看出来。而且在上面过程中，会看出来这个过程看似是连续的，但是计算机只能处理离散的东西。所以我么会更加关注这个交点的几个关键生命周期。称之为Significant Events，可以分之为两类： Circle Events 其运动效果如图，与SiteEvent不同的是，CircleEvent是同态生成的。其当两个交点合并成一个的时候触发。所以需要检测这种事件并将它加入到事件队列中去。 而这种相交事件显然只能发生在相邻的三段弧上，所以当我么产生新的弧时，有新的三段弧结构，我么就去检测制造这个Event并将其放在要处理的事件队列中去。 而这个事件位置的计算方式可以通过三段弧来确定。例如图，相邻的三段弧分别为$arc(i),arc(j),arc(k)$。这三段弧可以确定一个圆，而这个圆在扫描线上的最远点$c$，就是这个Circle Event该处理的位置。 首相对于$i,j,k$三点围城的圆，一定是空的。如果有其他点，三段弧一定不会相临。而当扫描线到$c$这个位置时，$i,j,k$三点到扫面线的距离都相等。这意味着原来两个Break Point已经合并成一个点。从Voronoi图角度来说，这个点就是Voronoi Vertex了。此时对于BeachLine的结果就是，两个交点融合，同时3段弧中间的那条消失。此时操作一共有六件事情： Site Events 先看如图情况，当扫描到一个点$s$时。可以想到这个时候会有一条从$s$发出的直线，刺穿BeachLine。这条线是一条退化的抛物线，可以想象后面这条包抛物线会逐渐向两边延展。这个向两边延展的过程，会画出一条新的edge产生。注意的是，一开始我们并不知道这条边时候结束。但是结合circle event，可以知道后面会有一个circle event给出这条edge的终止位置。即途中的$c$。 所以显然这个$s$点位置就是这个Site Event位置。当遇到其时，我么主要更新BeachLine以及加一条新的，未定向的Voronoi Edge。其具体操作如下： 以上就是两种主要的事件。可以确认这两种事件组成了所有的事件队列。对于每一种情况我们都通过Event处理了。现在来看复杂度分析。 复杂度分析 算法主要由以下几部分构成。 所以总体而言算法复杂度不超过$O(nlogn)$。
非空Cell特性：每个Site都会有一个非空的Cell
空圆特性：给定平面上任何一个点$x$。存在一个$Disk(s,x)$为空的。即Disk不包含任何Site。
最近邻Site皆共圆。
DCEL结构
Voronoi图复杂度分析
Voronoi图构造
Incremental Construction
Divide-And-Conqure
Plane-Sweep
