CS Geometry

前言

本文基于原书英文版《Computational Geometry: Algorithms and Applications(Mark de Berg, Otfried Cheong, Marc van Kreveld, Mark Overmars)》撰写，部分中文词汇参考清华大学出版社邓俊辉译本《计算几何：算法与应用》。

我尽量用自己的语言来总结归纳，以及给出部分我认为比较好的题目的答案，可能还会包含一些核心算法的具体 C++ 实现代码。

撰写本文的目的是深入学习理论计算几何，尤其是有深度的、有科研价值的课题。重点不在于完全掌握背诵所有算法，而是学习其中真正的思想内涵。

本书作者前言

每个章节都由一个问题开始，用于解决这个几何问题的概念和技巧就是这个章节真正的主题。给出应用的目的是给予读者启发。

目录

  • 计算几何-导言
  • 线段相交-专题图叠合
  • 多边形三角剖分-画廊看守
  • 线性规划-铸模制造
  • 正交区域查找-数据库查询
  • 点定位-找到自己的位置
  • Voronoi 图-邮局问题
  • 排列与对偶-光线追踪超采样
  • Delaunay 三角剖分-高度插值
  • 更多几何数据结构-截窗
  • 凸包-杂项
  • 空间二分-画家算法
  • 机器人运动规划-随意所之
  • 四叉树-非均匀网格生成
  • 可见性图-求最短路径
  • 单纯形区域查找-再论截窗

第一章：计算几何 - 导言 Computational Geometry - Introduction

三个例子：

  1. 走在大学里，有很多电话亭，你想找到最近的一个，拿到一张平面地图，可以利用 Voronoi 图将其划分；
  2. 知道要去哪个电话亭了，但是中间有建筑物障碍，拿到一张地图，但这时候需要让一个机器人利用程序去绕着障碍物走最短路；
  3. 获得了两张地图，一个是建筑物图，一个是道路图，你需要将他们叠在一起（Overlay），来展示两个图之间的组合信息。

1.1 简单例子：凸包 Convex Hulls

对于几何背景的算法问题，好的方案通常有两个要素：其几何性质和算法技巧，二者缺一不可。

关于二维凸包问题，至少我还算了解，其定义就是给定平面上的点集 $P = \{p_1, p_2, \cdots, p_n\}$，计算一个点集，能包含所有点，以顺时针输出，称其为凸包。

作者给出的第一个 SlowConvexHull(P) 确实有点蠢（甚至比我想的最笨的方法还要笨），先枚举任意两点 $p,q$ 形成的边 $\overrightarrow{pq}$，然后对于每个边，如果剩下的任意点 $r$ 都在这个边的右边，那么这个边就是凸包的一个边。很明显，这个算法是 $O(n^3)$ 的。

而且，我们在这个算法里还需要考虑，如果 $r$ 在边 $\overrightarrow{pq}$ 之上（lies on），我们该怎么做呢？这就是一个经典的退化情况。我们一般来说在一开始思考几何问题的时候会忽略这些问题，但是在实际情况下一定要引起重视（打 ACM 天天遇到 corner case）。针对上面这个算法，我们把边的判定条件改成“$r$ 在边的右边或在边上”即可。

同时，我们还忽略了精度问题，简单地认为我们可以 somehow 准确判断一个点在一条边的左边或者右边，但如果点是以浮点坐标展现，实际上用浮 点计算的时候会出现误差导致误判。

尽管我们已经证明了算法的正确性，以及处理特殊情况，但是它仍然不具有鲁棒性，小的计算错误可能会导致它以完全意想不到的方式失败。

鲁棒性 Robustness：健全的，耐用的。

当输入包含退化情况、坐标误差或浮点舍入误差时，算法仍能给出几何上和拓扑上一致、合法的结果，而不是因为一次符号误判就崩掉或输出矛盾结构。

然后作者又提出了一种求凸包的算法（Graham’s Scan），说实话这个方法是我当时第一次听到凸包这个概念时，想的类似方法。就是先把所有点按 x 轴排序，然后从最左边开始加点，每次加一个点，然后判断当前最尾端的三个点是否形成一个“右转”的折线，没有形成的话就删掉三个点中间的点。这样下来能找到整个凸包的上凸壳，再通过类似的算法找个下凸壳就好了。

这个算法里，每个点最多入队一次，出队一次，所以后面的部分复杂度是 $O(n)$
的，但是第一步排序的复杂度是 $O(n \log n)$ 的，所以还是一个 $O(n \log n)$ 的算法。

这个跟传统教的分治思想求凸包还不太一样，但其算法时间复杂度都是 $O(n \log n)$。

思考：能否用数学严谨证明求凸包算法的时间复杂度下限就是 $O(n \log n)$？

答案：可以，通过构造抛物线，将排序问题归约到凸包问题，已知排序问题下限是 $\Omega(n \log n)$。

退化和鲁棒性 Degeneracies and Robustness

一个计算几何算法通常有三个步骤：

  1. 忽略会影响几何概念的任何事情，尝试设计或理解一个算法；
  2. 调整算法设计来应对退化情况，例如用字典序等来处理某些情况，还有一个常用技巧是符号扰动；
  3. 实际实现，需要考虑到鲁棒性。

符号扰动 symbolic perturbation schemes：想象把输入点移动一个无限小的距离，使它们不再退化，但实际上不真正修改坐标，例如将 $x_i$ 变成 $x_i’ = x_i + \varepsilon_i$

1.3 应用领域

  • 计算机图形学：用来给计算机屏幕创造图片 / 模型，有二维 / 三维的不同应用，线段、多边形、物体、光线、材质等等；
  • 机器人学：主要是三维空间的问题，如何操纵机械臂，控制小车移动，很多时候会用到计算几何；
  • 地理信息系统：记录真实地理信息，模拟地形变化，城市轨道问题，管道系统问题等等，真实世界的地理是一个非常好的计算几何温床；
  • CAD / CAM：计算机辅助设计一类，常用于画图什么的，感觉是图形学的高级应用；
  • 其他应用：分子模型、模式认知，甚至数据库、数据结构；

练习

1.1 一个集合 $S$ 的凸包被定义为所有包含 $S$ 的凸集的交集。对于一个点集的凸包，它就是具有最小周长的凸集。我们想证明这两个条件是等价的。

a. 证明两个凸集的交集还是凸的。这能引申出有限数量凸集的交集都是凸的。
b. 证明包含一些点的最小周长多边形 $\mathcal{P}$ 是凸的。
c. 证明任意包含点集 $P$ 的凸集，包含最小周长多边形 $\mathcal{P}$。

答案

(a) 凸集 $S$ 的一个定义是：

$$
\forall x, y \in S, \forall \lambda \in [0,1]
$$

$$
\lambda x + (1-\lambda) y \in S
$$

其中，$x, y$ 都是点。

那么对于凸集 $A, B$，我们任取 $x, y \in A \cap B$，有对于 $\forall \lambda \in [0,1]$：

$$
\lambda x + (1-\lambda) y \in A
$$

$$
\lambda x + (1-\lambda) y \in B
$$

因此这个点同时属于 $A$ 和 $B$，也就是：

$$
\lambda x + (1-\lambda) y \in A \cap B
$$

所以 $A \cap B$ 也是凸集。

(b) 假设最小周长多边形 $\mathcal{P}$ 不是凸的，如下图实线所示：

那么一定会存在三个点，例如图中 $A, B, C$，使得凹顶点 $C$ 在内部形成三角形不等式

$$
|AC| + |BC| > |AB|
$$

此时我们发现，这个多边形并不是我们能围成的最小周长多边形，边折线 $ACB$ 可以替换成边 $AB$ ,与前提条件相矛盾，因此假设不成立。

所以最小周长多边形是凸的。

(c) 关键是要证明 $\mathcal{P}$ 的每一个顶点都属于 $P$。

如上图所示，假设 $\mathcal{P}$ 存在点 $A \notin P$，其隔壁两点为 $B$ 和 $C$，则我们在 $AB$ 和 $AC$ 上分别取足够靠近 $A$ 点的两点 $D$ 和 $E$，连接 $DE$。

因为我们可以使 $AD$ 和 $AE$ 非常小，让 $\triangle ADE$ 不包含 $P$ 的点，这样以来，根据三角形不等式

$$
|AD| + |AE| > |DE|
$$

我们可以把 $\mathcal{P}$ 的 $A$ 点换成 $DE$，从而得到更小的“最小周长多边形”，但这与原本所谓的“最小”冲突，因此假设不成立。

所以 $\mathcal{P}$ 上的顶点都属于 $P$，那么包含点集 $P$ 的凸包肯定也包含最小多边形 $\mathcal{P}$。

那么回到一开始的问题，证明最上面那两个条件是等价的，设凸集的交集是 $H$，那么 (a) 证明了 $H$ 是凸集，(b) 证明了 $H \subseteq \mathcal{P}$，(c) 证明了 $\mathcal{P} \subseteq H$，所以 $H=\mathcal{P}$，即两条件等价。

1.3 $E$ 是一个凸多边形未排序的 $n$ 条边的集合，描述一个 $O(n \log n)$ 的算法能顺时针输出这个凸多边形的顶点。

答案

如上图所示，我们只需要想象把所有边向量都移动到同一个点出发的线段，按照角度排序 $[0,2\pi)$，其沿着角度角度顺时针走的边，即在凸多边形上顺时针走的边，时间复杂度瓶颈在排序的 $O(n \log n)$。

实际操作可以通过边的两点的坐标算出这个角度。

1.4 对于前文提到的凸包算法，我们有能力检验一个点 $r$ 是否在一条线段 $pq$ 的左边或者右边。我们使 $p=(p_x, p_y), q=(q_x, q_y), r = (r_x, r_y)$。

a. 证明行列式

$$
D =
\begin{vmatrix}
1 & p_x & p_y \\
1 & q_x & q_y \\
1 & r_x & r_y
\end{vmatrix}
$$

的符号决定了 $r$ 在线段 $pq$ 的左边还右边。
b. 证明 $|D|$ 实际上是 $p,q,r$ 围成三角形面积的两倍。

答案

(a) 如上图所示，我们令

$$
P = (1, p_x, p_y) \\
Q = (1, q_x, q_y) \\
R = (1, r_x, r_y)
$$

然后有

$$
\det(D) = \det\left(D^T\right)
$$

且

$$
\det(A) = \text{volumn of the polytope formed by } a_1, a_2 \cdots, a_n \\
A = (a_1, a_2, \cdots, a_n)
$$

所以，我们可以得到 $\det(D)$ 就是由 $\overrightarrow{OP}, \overrightarrow{OQ}, \overrightarrow{OR}$ 张成的有向平行六面体体积，这里的有向指的是，$R$ 在 $PQ$ 的左右能决定这个体积的符号正负。

(b) 由平行六面体的性质得知，图中的三棱锥体积有

$$
V_{OPQR} = \frac{|D|}6
$$

且

$$
V_{OPQR} = \frac13 S_{\triangle PQR} h \\
h = 1
$$

所以：

$$
|D| = 2 S_{\triangle PQR}
$$

上面这个做法是我第一次接触这个题想到的做法，比较巧妙，也符合直观感觉，但是我跟 GPT 聊了一下，发现其实还有直接从行列式代数角度出发的做法。

行列式的两种展开方式：

(1) 通过行变换

$$
\begin{aligned}
D &=
\begin{vmatrix}
1 & p_x & p_y \\
1 & q_x & q_y \\
1 & r_x & r_y
\end{vmatrix} =
\begin{vmatrix}
1 & p_x & p_y \\
0 & q_x-p_x & q_y-p_y \\
0 & r_x-p_x & r_y-p_y
\end{vmatrix} =
\begin{vmatrix}
q_x-p_x & q_y-p_y \\
r_x-p_x & r_y-p_y
\end{vmatrix} \\
&= (q_x-p_x)(r_y-p_y) - (q_y-p_y)(r_x-p_x) \\
&= \overrightarrow{pq} \times \overrightarrow{pr}
\end{aligned}
$$

这里成功地将这个行列式化简成了两个二维向量的叉积，这个东西在计算几何里非常重要，其得到的数值也是具有方向性的，$\overrightarrow{pq} \times \overrightarrow{pr} > 0$ 意味着从 $\overrightarrow{pq}$ 逆时针旋转到 $\overrightarrow{pr}$，表示 $r$ 在 $\overrightarrow{pq}$ 的左边，反之同理。

其绝对值，也是二维行列式的值，即 $\overrightarrow{pq}$ 和 $\overrightarrow{pr}$ 围成的平行四边形的面积，可以很好的解决 (b) 问。

(2) 第二种方式是直接通过拉普拉斯展开

$$
D = q_x r_y - q_y r_x - p_x r_y + p_x q_y + p_y r_x - p_y q_x
$$

但这种方式在这里会容易卡住（我第一次就卡死了），因为需要补两项 $+p_x p_y - p_x p_y$ 才可以继续化简，但是并不是那么容易观察出来，接下来步骤与 (1) 相同，不多赘述。

第二章：线段相交 - 专题图叠合 Line Segment Intersection - Thematic Map Overlay

在地图应用里，地图一般都有所谓的“分层图”，比如河流、植被、公路、铁路等不同层次的网格叠加在地图上，然后用户也可以选择只看一个层次的内容。

在这种时候，不同层次的网格之间的交点，比如河流与道路的交点等，通常能引起我们的注意。

2.1 线段相交 Line Segement Intersaction

关于两张网格图叠在一起的线段相交问题，我们想知道有会相交在哪些点，经过一系列定义，我们可以形式化问题为：给定由 $n$ 个封闭线段组成的集合 $S$，求出这些线段的所有交点。

最容易想到的就是 $O(n^2)$ 暴力算法：对于每两个边，都求一下它们的交点。

但考虑到如果边没有排布那么密集，有些边离得特别远，这时候如果全都两两计算的话，会浪费非常多的时间，那么有没有一种办法，能在实际应用中，将时间复杂度降到比 $O(n^2)$ 要小很多，以至于说，我们可以将答案交点数量 $k$ 列入时间复杂度中，这种算法我们一般称为输出敏感算法 Output-Sensitive Algorithm。

关于输出敏感算法，我之前打的算法竞赛中，一般是不涉及相关内容的，这也是我想提到的，竞赛算法与科研算法的一些区别，科研中常用到的剪枝技巧、np-hard 问题的处理、输出敏感算法等都是竞赛学不到的内容。随着算法科研不断发展，传统多项式算法可能会越来越难挖掘，那么很多时候新的论文就会去研究非多项式问题里能不能将次方复杂度的实际运行时间压缩到一个可接受的范围。

这时候观察到线段相交的几何性质：相比于距离很远的两条线段，距离相近的线段更容易产生交点，或者说才是线段相交的“候选人”。

所以，作者在这里提出了一种扫描线算法，与其说是算法，不如说是一种在计算几何中常见的思想。虽然在二维平面中，实数坐标是无限的，看上去用一条直线扫描过去需要有无数的状态，但是实际上很多时候，我们只需要关心一些关键节点，比如线段的头和尾、线段交点等，这些关键节点的数量在大部分实际情况中是有限的。

想象一根从上往下扫描的直线，在中间的某个状态，它会与某些线段相交，而且我们可以试图维护“正在被直线相交的线段”里，哪些线段是所谓“相邻的线段”。

再往下考虑，我们发现，对于这种状态的变化，只在扫描线扫到线段头，线段尾，和线段交点，这三种情况时出现，作者将这三个点统称为 event point 事件点。

如上图所示，本来 $s_k$ 和 $s_m$ 是邻居，但是经过了 $s_k$ 和 $s_l$ 的交点后，发生了“邻居交换”。这个就是我们想要维护的一个核心问题。

还有一个事情是，刚开始的时候我们就获得了线段头尾的坐标，但我们还需要知道交点的具体坐标，这个东西并不能在扫描线扫到的时候自动检测（毕竟我们只扫事件点），但是可以在每次扫到一个线段头的时候，检测这个新的线段是否会与已知的左右两线段相交，然后立刻算出交点，将其加入事件点，如下图所示。

同样的，在扫描线到了线段尾，即要离开某个线段的时候，这个邻居状态也会发生变化。

到目前为止，我们已经对算法有了个初步的框架，接下来就是要处理细节和数据结构问题。

在计算几何中，与数据结构的联系实际上相当紧密，很多时候会出现利用数据结构优化计算几何算法。

首先，最容易想到的当然是维护一个 event queue $Q$，能按照顺序给我们吐出下一个事件点。这里有个细节就是如果一条边是水平的（头尾节点 $x$ 坐标相等），我们认为所有 $x$ 坐标小的排在前面，这个也是一个常用技巧，大部分时候是不影响结果的。

具体来讲，$Q$ 需要支持：

  1. 取出下一个事件点
  2. 插入新交点
  3. 检查点 $p$ 是否存在
  4. 按顺序储存

所以这里可以直接采用平衡二叉树，在 C++ 里使用 set 即可实现（？）。

其次，我们还需要维护一个状态结构 $T$，用于存储当前与扫描线相交的线段，按交点从左到右排序，也就是我们想要维护的“邻居”。这个结构需要支持：

  1. 插入线段
  2. 删除线段
  3. 找前驱 / 后继
  4. 改变顺序
  5. 找到经过事件点的线段

这里一看就知道用平衡二叉树了，非常适合。

快来参与讨论吧~