计算几何——凸包问题

Posted 2020-07-04

注：本文是2016年春季清华大学邓俊辉老师《计算几何》MOOC课程的简要个人总结系列之一，我将同步课程内容更新。不过有可能写的不完全是课程内容，也包含一些个人理解。如果你在看完本文后开始对计算几何感兴趣，请前往相应的MOOC平台完整学习邓老师的课程。如此精心设计和编排的课程，不应该被辜负。在此感谢邓老师！

Knowledge Dependence：

阅读本文前你只需要有基本的几何知识和算法知识（本篇重点要求复杂度分析与选择排序算法）即可。代码实现需要一丢丢C++基础。

由于作者很懒不想画图，所以在本篇中你需要结合文字表述和脑内小剧场，自己动手在纸上画图。如果不亲自动手，这一篇可能会相对难以理解。

 

我们仔细反思上一篇中极边（Extreme Edge）算法效率的低下的原因在于：对于我们所要求得的解，我们花费了太多无用的枚举在对结果不起贡献的有向线段上，也就是说，我们花费了大量时间在枚举非极边。

有没有一种办法让我们在找到的已知极边的基础上（此时所有极边链接起来是一条不间断的「折线段」），直接找到与当前我们找到的「折线段」的某一端点相连的下一条极边，直至我们的「折线段」首尾相接闭合呢？

 

Javis March 算法

我们同样规定逆时针方向为正方向、当前我们已经找到的极边构成「折线段」的端点是 k 、k 的前驱极点是 i（即我们刚刚找到的新极边的另一端点）。我们此时需要做的，就是立足于 k，找到下一条极边的另一个端点 s。

我们可以肯定的是，s 一定会在给定点集中除去 k 和 i 的其他点中产生，那么是什么原因可以让 s 脱颖而出被我们找到呢？

此时我们需要做一条辅助线。我们将我们刚刚找的极边 ik 延长（方向为从 i 指向 k 的方向，依然是逆时针），显然作为极边， ik 右侧一定是空的。再把 k 和候选的 s（candidate s）连接起来（方向从 k 指向 s ，依然逆时针），发现我们真正要找的 s （real s）会有这样一个特点：它是在所有 candidate s 中使得向量 ik 和 ks 夹角最小的 s，这样的 s 就是我们要找的 real s。

OK，到这里我们看似已经找到了算法的关键。但上述描述都还过于数学化，在我们具体实现之前，我们依然还需要回答三个同样关键的问题：

 

1. 我们如何判断计算夹角？或者无需计算具体值，但是一定能够比较夹角大小？

其实很简单，同样也是利用我们已经在前两篇中都用到的方法——ToLeft Test。具体方法为：假设我们现在已经找了某个局部最优的 candidate s，现在我们枚举下一个点 s‘。连接 ks（方向从 k 指向 s），用 s‘ 对其做 ToLeft Test，如果 s‘ 确实在 ks 的左边，那么 s 依然是比 s‘ 更优的点。但如果反之，s‘ 在 ks 的右边，那么 s‘ 将会以更优的身份的身份取代 s 成为新的局部最优的 candidate s。一次次迭代，直到枚举完所有的点，我们就找到了我们要找的 real s。

 

2. 我们如何找到初始情况下的极点 k 呢？只有立足于此点，我们才能开始「亦步亦趋」地构造我们想要找的「折线段」。

同样不难。从人类视觉上来说，我们总能一下子找到平面上这堆点的最「角落」的那个点，计算机同样也可以。不失一般性，我们可以找到 Lowest-Then-Leftmost 的那个点。具体来说，就是 y 值最小的点中 x 值也最小的点。感性地说，就是最左下角落的点。可以发现，Lowest-Then-Leftmost Point 必然是一个极点。我们将这个点简称为 LTL。

 

3. 因为我们确定的第一个极点 k 还并没有前驱 i，那我们如何找到第一条作为角度参考线的 ik 呢？

其实我们并不需要找到这样的 i。因为 k 已经是「最低」的那个点了，其他所有点都至少不比它低，所以我们可以直接把水平参考线，也就是 x 轴的方向（即水平正无穷方向），规定为初始情况下的角度参考线。以此我们可以找到对于初始状态下的 s。

 

至此，我们已经可以给出具体的代码实现：

 1 /*******************************
 2 * Jarvis_March_Algorithm for Convex Hull construction
 3 * Time Complexity:    O(n^2)
 4 ********************************/
 5 
 6 void Jarvis(Point S[], int n)
 7 {
 8     // 初始化：“有罪推论”
 9     for (int k = 0; k < n; k++) S[k].extreme = false;
10     // 找到 Lowest-Then-Leftmost Point
11     int ltl = LTL(S, n); int k = ltl;
12     do {
13 
14     // k 一定是极点，s 待求
15         S[k].extreme = true; int s = -1;
16         // 枚举所有 candidate s
17         for (int t = 0; t < n; t++)
18         {
19             // 除 k、s 外，如果在右侧或者是初始情况（s == -1），则更新 s
20             if (t != k && t != s &&
21                 (s == -1 || !ToLeft(S[k], S[s], S[t])))
22                 s = t;
23         }
24         // 将 k 的后继更新为 s，并用 s 迭代 k
25         S[k].succ = s; k = s;
26     } while (ltl != k); // 如果 k 已经是 ltl 了，则折线段已经闭合，凸包已构成
27 }
28 
29 int LTL(Point S[], int n) {
30     // 初始化 ltl 为下标为 0 的点
31     int ltl = 0;
32     for (int k = 1; k < n; k++)
33     {
34         // 更低 或者 在同样低时更左，则更新 ltl
35         if (S[k].y < S[ltl].y ||
36             (S[k].y == S[ltl].y && S[k].x < S[ltl].x))
37             ltl = k;
38     }
39     return ltl;
40 }
41 
42 bool ToLeft(Point p, Point q, Point s) {
43     // 将 ToLeft Test 等价为求2倍“方向面积”，避免除法或三角运算带来的精度问题
44     return Area2(p, q, s) > 0;
45 }
46 
47 int Area2(Point p, Point q, Point s) {
48     // 2倍“有向面积”
49     return
50         p.x * q.y - p.y * q.x
51         + q.x * s.y - q.y * s.x
52         + s.x * p.y - s.y * p.x;
53 }

在实现之后，我们反观此算法，有两点需要点明：

1. 其实这个算法的核心思想跟我们特别熟悉的另一个算法的核心思想是完全一致的，那就是选择排序。两个算法都是不断维护一个局部解的增量思想，一步一步达到全局解；都是从当前未知解中选取某个标准下的最值加入到局部解的边界。唯一不同的是，一个是一维数据，一个是二维的数据。

2. 这是一个输出敏感型算法（Output Sensitivity）：

• • 如果我们进行时间复杂度分析，可以发现初始化 LTL 为 O(n)，主体部分的 do-while 循环是 O(n^2) 的，也就是该算法在最坏情况下也不过是 O(n^2) 的。但是什么情况下才会产生严格的最坏情况呢？答案是当给定点集中的所有点都是极点的时候，这时候我们需要构造一个凸包就外围 do-while 循环就需要「行进」n 步。

  • 可是，相对应的，我们同样存在最好的情况：也就是即便我们有 n 个点，最后构成的凸包依然是一个三角形时，此时外围 do-while 循环只需要简单地走3步即可闭合起来，求得凸包。最优情况下竟然达到了线性 O(n)！

  • 因此，我们说，这是一个输出敏感型算法，也就是说，我们算法的具体的复杂度是由算法输出的规模来决定的。如果我们定义最后得到的凸包有 h 条极边，那么我们可以准确的说这个算法的时间复杂度就是 O(nh) 的。如果你足够敏感，你可能会觉得这是一个经典的「先有鸡，还是先有蛋」的问题。其实不然，一旦问题一给出，那么对于这个点集的凸包就已经唯一确定了，h 也就已经确定了，只是我们可能还没具体求解出来而已。

回顾一下，我们至此已经了解了三个凸包构造算法，并且每个算法都比上一个算法在时间复杂度上要优一个数量级。

就像我们研究排序算法那样，此时我们不禁要问，凸包构造算法的理论时间复杂度是多少呢？

下一篇我们就简要分析这个问题。

 

【To Be Continued】