5.5 应用细节
(Applications in Detail)
接下来我们将讨论维诺图的五个应用,不过详细程度会有所不同: 最近邻、"胖(fat)"三角分割、最大空圆、最小生成树和旅行商路径。
5.5.1. 最近邻(Nearest Neighbors)
我们已经在5.1 节中提及了维诺图在最近邻聚类问题中的应用。 该问题可视为一个查询(query)问题: 离查询点最近的邻居是哪个? 另一种变体是全最近邻 (all nearest neighbors)问题: 找到给定集合中每个点的最近邻。 这在许多领域都有应用,包括生物学、生态学、地理学和物理学。
给定点集 \(P\),当且仅当 \(|a - b| \le \min_{c \neq a} |a - c|\) 时,我们说 \(b\) 是 \(a\) 的最近邻,其中 \(c \in P\)。 我们可以将此关系写为 \(a \rightarrow b\)。注意,该定义关于 \(a\) 和 \(b\) 所扮演的角色是不对称的,这表明该关系本身也不是对称的。 事实上确实如此,如果 \(a \rightarrow b\),并不一定有 \(b \rightarrow a\),如图 5.12 所示。还要注意,一个点可能有多个距离相等的最近邻,例如图中的点 \(d\)。
查询最近邻
给定一个固定的点集 \(P\),在 \(O(n \log n)\) 时间内构建维诺多边形。现在对于查询点 \(q\),其最近邻问题,可以归结为找出它落在哪个(或哪些)维诺多边形区域内, 因为这些维诺多边形区域的站点正是它的最近邻。在一个区域内部定位一个点的问题被称为点定位 point location。这个问题已得到了广泛的研究, 我们将在第 7.11 节中详细讨论它。我们将看到,在这种情况下,查询的事件复杂度是 \(O(\log n)\)。
全最近邻
对于点集 \(P\) 中的每一个点,我们都用一个节点与之关联,并以最近邻的关系在节点之间建立连边,如此构成的图被称为最近邻图(Nearest Neighbor Graph, NNG)。 虽然最近邻得关系是不对称的,我们仍将其定义为无向图。它完全可以被定义为有向图,但在这里我们不需要。
下面的引理给出了一种简洁的最近邻算法,我将证明留给练习 5.5.6[2] 和 [3]。
引理 5.5.1. \(NNG \subseteq D(P)\)
在一个集合中为每个点寻找最近邻的暴力算法需要 \(O(n^2)\) 的时间,但上述引理让我们只需搜索 Delaunay 三角分割的 \(O(n)\) 条边,从而达到 \(O(n \log n)\) 的时间复杂度。
5.5.2 最大化最小角的三角分割 (Triangulation Maximizing the Minimum Angle)
人们通常使用有限元(finite element analysis)来分析复杂形状的结构特性。例如,汽车制造商使用这种方法来对车身进行建模(Field 1986)。 将复杂形状被划分为一个有限元网格,然后在这些网格划分上离散化,来求解动力学的相关微分方程。所用数值过程的稳定性取决于划分的质量,而事实证明, Delaunay 三角剖分是特别好的划分。我们现在将讨论 Delaunay 三角分割在何种意义上是好的。
点集 \(S\) 的一个三角分割是该对象的泛化,而 Delaunay 三角分割是它的一个特例,它的分割线段的端点都在 \(S\) 上,这些线段仅在端点处相交, 并且将 \(S\) 的凸包划分为三角形。就有限元分析而言,"胖(fat)"三角形的分割是最好的。更精确的说法是要,避免含有小角度的三角形。 自然,我们要寻找一个拥有最大最小角的三角分割,也就是说,在所有三角分割中最大化最小角。事实证明,这正是 Delaunay 三角分割。 我们引入一些符号后,再来描述它,可以给出一个更精确的陈述。
设 \(T\) 是点集 \(S\) 的一个三角分割,角度序列 \((\alpha_1, \alpha_2, \dots, \alpha_{3t})\) 是其三角形角度的列表,按从小到大排序,其中 \(t\) 是 \(T\) 中三角形的数量。 数量 \(t\) 对于每个 \(S\) 都是一个常数(练习 5.5.6[4])。我们可以定义同一点集的两个三角分割 \(T\) 和 \(T'\) 之间的关系,来描述三角形有多胖。 如果说 \(T \ge T'\),即\(T\) 比 \(T'\) 更胖,是指 \(T\) 的角度序列在字典序(lexicographically)上大于 \(T'\) 的角度序列。 即要么 \(\alpha_1 > \alpha'_1\),要么 \(\alpha_1 = \alpha'_1\) 且 \(\alpha_2 > \alpha'_2\),要么 \(\alpha_1 = \alpha'_1\) 且 \(\alpha_2 = \alpha'_2\) 且 \(\alpha_3 > \alpha'_3\),依此类推。
Edelsbrunner(1987,p. 302)证明了如下定理。这表明 Delaunay 三角分割最大化了最小角度。
定理 5.5.2. Delaunay 三角分割 \(T = D(P)\) 是最胖的,对于 \(P\) 的任何其他三角分割 \(T'\),都有 \(T \ge T'\)。
5.5.3. 最大空圆 (Largest Empty Circle)
我们在 5.1 节中提到了在一组站点 \(S\) 中寻找最大空圆的问题。 这样一个圆的圆心是开设新商店的好位置。Toussaint (1983a) 提到了另一个应用,尽可能远离一组城市站点来定位核反应堆。我们现在详细研究最大空圆问题。 除非对圆心的位置施加某些约束,否则这个问题几乎没有意义,因为在任何有限点集之外,总是存在任意大的空圆。因此,我们将这一问题表述成如下的形式:
最大空圆(Largest Empty Circle)。寻找一个最大的空圆,其圆心位于包含 \(n\) 个站点 \(S\) 的闭合凸包内, 所谓"空"是指其内部不包含任何站点,所谓"最大"是指不存在其他半径严格更大的此类圆。
设 \(f(p)\) 为以点 \(p\) 为圆心的最大空圆的半径,我们要寻找的是该函数在 \(S\) 的凸包 \(H = \mathcal{H}(S)\) 内所有 \(p\) 上的最大值。 但是,这些最大值的候选点似乎有无限多个。计算几何中的一个常见操作是将无限的候选集缩减为一个小的有限列表,然后高效地找到它们。 本节中我们遵循这一思路,首先非正式地论证只有某些特定的点 \(p\) 才是 \(f\) 最大值的真正候选点。
凸包内部的圆心 (Centers Inside the Hull)
想象一下从凸包 \(H\) 内的一点 \(p\) 开始膨胀画圆。这个圆第一次碰到站点集 \(S = \{s_1, \dots, s_n\}\) 中的某个站点时的半径,就是 \(f(p)\) 的值。 在本小节中,我们暂时假设 \(p\) 严格位于 \(H\) 的内部。如果半径为 \(f(p)\) 时,圆内只包含一个站点 \(s_1\),那么很明显 \(f(p)\) 不可能是半径函数的最大值。 因为如果 \(p\) 沿着从 \(s_1\) 穿过 \(p\) 的射线 \(s_1p\),向远离 \(s_1\) 的方向移动到 \(p'\),如图 5.13 中上方的圆所示,\(f(p')\) 会更大。 因此,\(p\) 不可能是 \(f\) 的局部最大值,因为在 \(p\) 的任何邻域内都存在一个点 \(p'\) 使得 \(f\) 的值更大。 注意,假设 \(p\) 严格位于凸包内部,可以保证存在如上所述的同样位于 \(H\) 内的点 \(p'\)。
现在我们假设在半径为 \(f(p)\) 时,圆恰好包含两个站点 \(s_1\) 和 \(s_2\)。此时,\(f(p)\) 不可能取到最大值。如果 \(p\) 沿着 \(s_1s_2\) 的垂直平分线,向远离线段 \(s_1s_2\) 的方向移动, 如图 5.13 下方的圆所示,\(f(p')\) 会再次变大。理解这一点的另一种方式是通过第 5.4.4 节讨论的以站点为中心的圆锥的交集。 两个这样的圆锥的交线(图 5.9)代表了垂直平分线上的点到站点的距离。由于这条曲线是向上的双曲线分支,垂直平分线的任何内点都不是局部最大值,在两个可能的方向之一,距离总是增加的。
只有当圆包含三个站点时,\(f(p)\) 才可能达到最大值。如果这三个站点"跨越(straddle)"圆心 \(p\),即它们张开的范围超过半圆(如图 5.3 所示), 那么 \(p\) 向任何方向的移动都会导致 \(p\) 离某个站点更近,从而减小 \(f(p)\)。现在我们已经确认了这一事实:
引理 5.5.3. 如果最大空圆的圆心 \(p\) 严格位于站点凸包 \(\mathcal{H}(S)\) 的内部,那么 \(p\) 必定与某个维诺顶点重合。
注意,并不一定每个维诺顶点都代表 \(f(p)\) 的局部最大值(练习 5.5.6[5])。
凸包上的圆心 (Centers on the Hull)
现在我们来看一下圆心 \(p\) 直接位于凸包 \(H = \mathcal{H}(S)\) 上的情况。我们之前的论证之所以不适用,是因为将 \(p\) 移动到 \(p'\) 可能会使其移出凸包, 而我们需要将圆心约束在了凸包内。我们现在论证,一个最大圆必须包含两个站点。
假设 \(f(p)\) 是最大值,其中 \(p\) 在 \(H\) 上,且圆只包含一个站点 \(s_1\)。首先,\(p\) 不可能位于 \(H\) 的顶点上,因为 \(H\) 的顶点都是站点本身,这将意味着 \(f(p) = 0\)。 因此,\(p\) 位于 \(H\) 的某条边 \(h\) 的内部。那么,沿着 \(h\) 向某一方向移动 \(p\) 必然会增大其与 \(s_1\) 的距离,如图 5.14 所示。 我们可以通过想象以 \(s_1\) 为顶点的圆锥被一个垂直平面(图 5.9)切割,可以直观地看到这一点。
然而,如果以 \(p\) 为圆心的圆包含两个站点 \(s_1\) 和 \(s_2\),那么有可能沿着站点角平分线增加距离的方向是指向凸包外部的。因此,\(f(p)\) 完全有可能处于局部最大值。我们已经证明了如下事实:
引理 5.5.4. 如果最大空圆的圆心 \(p\) 位于站点的凸包 \(\mathcal{H}(S)\) 上,那么 \(p\) 必定位于一条维诺边上。
算法
现在我们找到了一组有限的点,它们可能是最大空圆的圆心。即维诺顶点,以及维诺边与站点凸包的交点。这引出了Algorithm 5.1 中的算法,该算法由 Toussaint(1983a) 提出。
注意,并非每个维诺顶点都一定在凸包内部,如图 5.14所示,这就需要在算法中进行 \(v \in H\) 的检查。该算法的一个朴素实现需要 \(O(n^2)\) 的时间, 但在 \(H\) 中定位一个维诺顶点以及求维诺边与 \(e\) 的交点都可以在 \(O(\log n)\) 时间内完成,这些优化使得整体算法的时间复杂度达到 \(O(n \log n)\)。我们将细节留给练习 5.5.6[6]。
5.5.4. 最小生成树 (Minimum Spanning Tree)
一个点集的的最小生成树(Minimum Spanning Tree, MST) 是能够连接所有点的并且边长最小的树,或者说是,一棵最短的树其节点恰好是该集合中的点。 如果我们按照欧几里得长度来计算线段两端之间的距离,那么这棵树被称为欧几里得最小生成树,缩写为 EMST。在这里,我们将只考虑欧几里得长度,因此将省略这个冗余的修饰词。 图 5.15 给出了一个欧几里得最小生成树的例子。MST 有许多应用。例如,许多局域网采用树形结构连接主机节点。MST 是使总导线长度最小化的网络拓扑,这通常也能使成本和时间延迟最小化。
Kruskal 算法
在这里,我们考虑计算平面上一组点的 MST 的问题。我们先来看看为图 \(G\) 计算 MST 这一更一般的问题。虽然不是那么直观,但却可以用一种简单的贪心策略可以找到 MST。 这种算法基于一个简单的直觉,最短的树应该由最短的边组成。我们可以通过不断添加尚未探索过的最短边来逐步构建这样一棵树,这个过程中需要保持树的性质——无环。 该算法被称为 Kruskal 算法,最早可追溯到 1956 年。
设 \(T\) 为逐步构建的树,并用符号 \(T + e\) 表示树 \(T\) 与边 \(e\) 的并集。Kruskal 算法如算法 5.2 所示。我们不会停下来证明该算法的正确性, 而是断言其复杂度由其一开始的排序主导的。这需要 \(O(E \log E)\) 的时间,其中 \(E\) 是图中的边数。
\(\text{MST} \subseteq \mathcal{D}(P)\)
平面上一组点,共有 \(\binom{n}{2}\) 条可能的边,因此如果在完全图上进行排序,排序的复杂度为 \(O(n^2 \log n)\)。但回想一下,Delaunay 三角分割的边在某种意义上记录了邻近信息, 因此我们希望只需要使用 Delaunay 边就可以构建 MST。幸运的是这是可以的,正如以下定理所示。
定理 5.5.5. 最小生成树是 Delaunay 三角分割的子集,即,\(\text{MST} \subseteq \mathcal{D}(P)\)。
证明:我们想证明如果 \(ab \in \text{MST}\),那么 \(ab \in \mathcal{D}\)。假设 \(ab \in \text{MST}\) 且假设相反的情况 \(ab \notin \mathcal{D}\)。 然后我们采用反证法,通过证明所谓的假设 MST 不是最小的来推出矛盾。
如果 \(ab \in \mathcal{D}\),那么存在经过 \(a\) 和 \(b\) 的空圆(根据性质 P7 和定理 5.3.1)。 所以如果 \(ab \notin \mathcal{D}\),那么经过 \(a\) 和 \(b\) 的圆就不是空圆。因此在 \(ab\) 为直径的圆上或圆内,一定有一个点 \(c\)。
所以假设 \(c\) 在以 \(ab\) 为直径的圆上或圆内,如图 5.16 所示。那么 \(|ac| < |ab|\) 且 \(|bc| < |ab|\),即使 \(c\) 在圆上这些不等式也成立。 从所谓的 MST 中移除 \(ab\) 会导致分裂成两棵树,\(a\) 在其中一棵 \(T_a\) 中,\(b\) 在另一棵 \(T_b\) 中。不失一般性的,假设 \(c\) 在 \(T_a\) 中。 移除 \(ab\) 并添加 \(bc\) 会得到一棵新树 \(T' = T - ab + bc\)。这棵树更短,因此使用 \(ab\) 的树不能是最小的,这就矛盾了,假设\(ab \notin \mathcal{D}\)是不成立的。 所以必须有 \(ab \in \mathcal{D}\)。□
这就产生了一个改进的 Kruskal 算法。首先找到 Delaunay 三角分割,时间复杂度为 \(O(n \log n)\),然后对其中的 \(O(n)\) 条边进行排序。 事实证明,Kruskal 算法的其余部分可以在 \(O(n \log n)\) 的时间内运行,因此对于平面上的 \(n\) 个点,找到 MST 的总复杂度是 \(O(n \log n)\)。
5.5.5. 旅行商问题 (Travelling Salesperson Problem)
旅行商问题是计算机科学中研究最广泛的问题之一。在给定的点集中找到一条遍历每个点的最短闭合路径,这样的路径被称为旅行商路径(traveling salesperson path, TSP)。 可以想象这些点是销售员在回家之前必须按任意顺序访问的城市。这个问题具有巨大的实际意义,不仅因为这一个应用,还因为许多其他问题可以归约为它。 不幸的是,该问题已被证明是 NP-hard 的,这是一个专业术语,意味着目前没有已知的多项式算法可以解决它(Garey & Johnson 1979)。而且在撰写本文时,似乎也不太可能找到这样的算法。 实际意义和难解性的结合促使人们寻找有效的启发式和近似算法。基于 Delaunay 三角分割通过最小生成树来实现,是最简单的近似算法之一。
这个想法虽然很简单,但它能干呀。找到点集的最小生成树(MST),然后简单地按照图 5.17 所示的方式往返跟随它。这样构建的回路长度正好是 MST 长度的两倍,因为树的每条边都在每个方向上遍历了一次。
现在我们来计算一个界,看看这个 double-MST 回路有多差。设 \(M\) 为最小生成树的长度,\(M_2\) 为 double-MST 的长度。当然 \(M_2 = 2M\)。 设 \(T\) 为旅行商路径的长度,\(T_1\) 为去掉一条边后的 TSP 路径的长度。注意,\(T_1\) 是一棵生成树。显然有以下不等式成立:
$$ T_1 \lt T $$ $$ M \le T_1 $$ $$ M \lt T $$ $$ M_2 \lt 2T $$如此我们得到了一个上界: double-MST 的长度不会超过最优 TSP 长度的两倍。
这个结果可以通过各种启发式方法改进。我只概述最简单的一种,不重复访问同一个地点。从起点开始遍历 double-MST 路径,并加以修改。 如果下一个地点已经被访问过,则跳过该地点,并考虑连接到 double-MST 回路上的下一个地点。这具有抄近道去某些地点的效果。 如果我们按照在 double-MST 回路中被访问的顺序对地点进行索引,那么某个地点 \(s_i\) 可能会通过一条直线段连接到 \(s_j\),而在 double-MST 回路中, 将是一条弯曲路径 \(s_i, s_{i+1}, \dots, s_{j-1}, s_j\)。根据三角不等式,直线路径总是比弯曲路径短,因此这种启发式方法只能缩短路径。 如图 5.18 所示。注意,缩短后的路径可能会自相交。
不幸的是,这种启发式方法并不能保证性能有所提升。但一种 Christofides 的启发式规则可以。它使用一组被称为"最小欧几里得匹配 minimum Euclidean matching"的线段作为抄近道的指导, 能保证路径长度不超过最优解的 \((3/2)T\),也就是说,误差不会超过 50%。更复杂的启发式方法通常能在最优解的很小百分比内找到路径(Bentley 1992), 尽管对于算法本身而言,这一点并没有得到保证。最近的一个令人兴奋的理论突破是一种"多项式时间近似方案 polynomial-time approximation scheme",用于 TSP。 在这种方案中,人们所能期望的最好结果是 NP-complete 问题的近似解。这是一种在 \((1+\epsilon)\) 倍最优解范围内进行拟合的方法,时间复杂度为 \(O(n^k)\), 其中 \(p\) 与 \(1/\epsilon\) 成正比(参见 Arora (1996) 和 Mitchell (1996))。
5.5.6. 练习
- NNG 的度数 Degree of NNG。有向最近邻图(directed Nearest Neighbor Graph, NNG)中,\(n\) 个点的节点在二维空间中的最大出度是多少? 最大入度是多少? 请给出达到你答案的示例,并尝试证明它们确实是最大值。
- [easy] NNG 和 \(\mathcal{D}\)。给出一个例子说明 NNG 可以是 \(\mathcal{D}(P)\) 的真子集。
- \(\text{NNG} \subseteq D \)。 证明引理 5.5.1。如果 \(b\) 是 \(a\) 的最近邻,那么 \(ab \in \mathcal{D}(P)\)。
- 三角剖分中的三角形数量 Number of triangles in a triangulation。证明对于某些固定的点集 \(S\),其任何三角分割的三角形数量 \(t\)都是一个常数,即, \(S\) 的所有三角剖分都有相同的 \(t\)。
- 维诺顶点不是局部最大值 Voronoi vertex not a local max。构造一个点集,使其拥有一个严格位于凸包内部的维诺顶点 \(p\),使得 \(f(p)\) 不是局部最大值, 其中 \(f\) 是 5.5.3 节中定义的半径函数。
- 空圆算法 Empty circle algorithm。用伪代码详述如何实现空圆算法(算法 5.1),其时间复杂度为 \(O(n \log n)\)。
- 相对邻域图 Relative Neighborhood Graph。点集 \(p_1, \dots, p_n\) 的相对邻域图(RNG)是一个图。其节点就是这些点,
并且当且仅当 \(p_i\) 和 \(p_j\) 彼此之间的距离至少与它们到任何其他点的距离一样近时,两个节点 \(p_i\) 和 \(p_j\) 才由一条弧连接,即,如果
$$
|p_i - p_j| \le \min_{m \neq i, j} \max(|p_i - p_m|, |p_j - p_m|)
$$
参见 Jaromczyk & Toussaint (1992)。该方程确定了 \(p_i\) 和 \(p_j\) 相邻时的一个"禁止 forbidden"区域,点 \(p_m\) 不可能位于该区域内。
这与定理 5.3.1 不无相似之处。该区域称为 \(Lune(p_i, p_j)\),它是以 \(p_i\) 和 \(p_j\) 为中心、半径为 \(|p_i - p_j|\) 的两个开圆盘的交集。
- 设计一个"暴力 brute-force"算法来构建 RNG。不必担心效率。它的时间复杂度是多少?
- 证明 \(RNG \subseteq \mathcal{D}(P)\)。RNG 的每一条边也是 Delaunay 三角分割的一条边。(与定理 5.5.5 对比。)
- 使用 b 来设计一个更快的算法。
- 三维 Delaunay 三角分割的大小 Size of Delaunay triangulation in three dimensions。我们已经证明过,二维 Delaunay 三角分割的大小是线性的,即 \(O(n)\)。
证明在三维空间中并非如此。三维空间中 \(\mathcal{D}(P)\) 的大小可以是二次的。参考二维版本定义三维空间的 \(\mathcal{D}(P)\)。它是 \(\mathcal{V}(P)\) 的对偶,
是没有唯一最近邻的点的轨迹(locus)。
设 \(P\) 是一个由两部分组成的点集:
- \(n/2\) 个点均匀分布在 \(xy\) 平面上以原点为中心的一个圆上
- \(n/2\) 个点均匀分布 \(z\) 轴上,关于原点对称。
- 三维空间中相对近邻图的大小 Size of Relative Neighborhood Graph in three dimensions。练习[7] 已经证明二维空间中 \(RNG \subseteq \mathcal{D}(P)\) 成立。 这种包含关系在任意维度上都成立。Jaromczyk & Toussaint 1992 已经证明三维 RNG 的大小是 \(O(n^{4/3})\),因此它比 Delaunay 三角分割更小。 但它的上界看起来很松,Jaromczyk & Kowaluk (1991) 猜想它的大小是 \(O(n)\)。确认这一猜想是一个开放问题。 尝试确定练习 [8] 中例子的 RNG 其 \(\mathcal{D}(P)\) 是二次的。
- \(MST \subseteq RNG\)。证明 MST 的每一条边也是 RNG 的一条边。(与定理 5.5.5 对比。)
- 最远点维诺图 Furthest-point Voronoi diagram。定义点集 \(P\) 的最远点维诺图 \(\mathcal{F}(P)\),将平面的每个点关联到它的"最远邻域 (furthest neighbor)",即,离它最远的站点。
指向一个最远邻域的点的区域形成最远点维诺区域,这些区域的边界形成 \(\mathcal{F}(P)\) 的边。参见图 5.19。
- 两个站点的 \(\mathcal{F}(P)\) 是什么?
- 三个站点的 \(\mathcal{F}(P)\) 是什么?
- 推导最远点维诺图的一些结构性质,类似于 5.3.1 节和 5.3.2 节中 Delaunay 和维诺图的性质。可以借助图 5.19 来提出假设。
- 最小生成圆 Minimum spanning circle。证明最远点维诺图可用于计算包围给定点集的"最小半径圆"。假设 \(\mathcal{F}(P)\) 是可用的。
