计算几何——凸包问题

Posted 2020-07-03

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

 

Knowledge Dependence：

阅读本文前你只需要有基本的几何知识和算法知识即可。代码实现需要一丢丢C++基础。由于作者很懒不想画图，所以你还需要一点脑内小剧场。不过相信我，自己在脑中一步一步地形象化问题将会很有意思。不信？往下看咯。

写在前面

计算几何就是几何版本的「算法与数据结构」，此话不假。她即有算法的抽象性和逻辑性，又有几何的形象与直观。这种兼具两种看似矛盾特性的特点，使得她成为了一门十分有趣味性的学科。

「计算几何」，重在「计算」。即我们如借助计算机，选取高效的算法，来解决几何问题。看，最终依然落脚在了「算法」上，所以要是说「计算几何」是包上了「几何」外衣的「算法设计」，完全可以。

既然研究的是算法，那一定也是很多更高层学科的基础。是的，诸如图形学、CAD，甚至是人工智能的路径规划等问题，都可以看到计算几何的身影。

每个学科都有其重要且基础的基本问题，计算几何也不例外，凸包问题便是在计算几何占据这样地位的问题。

接下来我将介绍凸包问题。邓老师曾说过，他希望他的课高中毕业生也能听懂，我也将尽我所能以精简明晰的语言来介绍。

 

什么是凸包（Convex Hull）

严格的数学定义是乏味难懂的，简单的说，在一面墙上的一块区域钉了很多钉子，你用一条橡皮筋，把它扯大使其能够包括所有的钉子，然后松手，此时橡皮筋形成的多边形就是凸包。

可以简单地看出，凸包至少有这样两个特点：

• 封闭的简单多边形（中间不含孔洞）；
• 所有的内角都小于180°；

所谓的凸包问题，就是在给定这些「钉子」（点集）的情况下，如何快速求出「橡皮筋」（凸包）的问题。

 

极点（Extreme Point）算法

什么是极点？

很简单，参与了固定「橡皮筋」的「钉子」就是极点。

现在我们把注意力放到极点上。为什么呢？因为如果我们求出了所有的极点自然就已经求出凸包了吖。

那如何确定给定点集中哪些点会是极点呢？或者换一个问法是，极点和其他点有什么本质上的不同呢？

很清楚地我们可以发现，极点不会被包含在由任何其他三个点构成的三角形中，而非极点一定可以找到由其他三个点组成的三角形，将它包围起来。

所以，很自然地，我们可以得到这样的找极点算法：

枚举可以由给定点集组成的所有的三角形，对于每一个这样的三角形，再枚举所有点（不包含组成这样三角形的顶点），如果对于当前点，它落在这个三角形内，那么我们可以立即判定它一定不是极点，可以排除在外。最后在所有的枚举完成之后，所有未被我们排除的点即构成极点。我们此时也就得到了凸包。

看似很清楚简明，完了没？

……还有一个问题，我们如何编码实现判断一个点是不是落在某个三角形之类呢？

稍加思考我们就会发现，我们需要把问题进一步转化成跟简单的问题。也就是，判断出点与三角形之间的关系。

我们定义逆时针为正方向，把三角形的三条边都看成有向线段。

此时，如果我们可以判断出点都在三角形三条边的左侧，则该点一定是落在三角形内的。可以一一枚举验证，反之亦反。

是的，问题被更加简化了。

再来，我们如何判断一个点在有向线段的左侧？

学过解析几何的朋友一定对这个问题不陌生，「有向面积」（Directed Area）这个概念就可以拯救我们（如果你还不清楚，see：有向面积_百度百科）。在我们当前定义的正方向的前提下，如果点和有向线段的两个断点之间构成的有向面积是正的，那么这个点就一定在这个有向线段的左侧。注意，这个方法在我们后面的算法中依然会经常用到，所以如果你对此还有疑惑，请一定去了解清楚。

OK，至此，我们的问题已经全部解决了。

以上所有分析问题的过程，无处不体现了数学中的「化归」思想，或者说算法中「分治」思想，即把一个大问题分解成更易解决的小问题，层层化归，最后变成我们可以解决的问题，原问题自然就解决了。这，也就是我们中国人常说的「大事化小，小事化了」。你会发现，后面我们依然会经常运用这种思想。

代码实现：（每一步都非常清晰，且有注释，建议结合代码再梳理一遍思路）

 1 /*******************************
 2  * Extreme_Point_Algorithm for Convex Hull construction
 3  * Time Complexity:    O(n^4)
 4 ********************************/
 5 
 6 void extremePoint(Point s[], int n) {
 7     // 初始化：“无罪推论”，认为所有的点都是极点，再一一排除
 8     for (int s = 0; s < n; s++) S[s].extreme = true;
 9 
10     //枚举所有的三角形
11     for (int p = 0; p < n; p++)
12         for (int q = p + 1; q < n; q++)
13             for (int r = q + 1; r < n; r++) {
14                 // 枚举所有点
15                 for (int s = 0; s < n; s++) {
16                     // 去掉顶点和已经被判定为非极点的点
17                     if (s == p || s == q || s == r || !S[s].extreme)
18                         continue;
19                     // 进行 InTriangle Test
20                     if (InTriangle(S[p], S[q], S[r], S[s]))
21                         S[s].extreme = false;
22                 }
23             }
24 }
25 
26 bool InTriangle(Point p, Point q, Point r, Point s) {
27     // 分别对三角形每条边做 ToLeft Test
28     bool pqLeft = ToLeft(p, q, s);
29     bool qrLeft = ToLeft(q, r, s);
30     bool rpLeft = ToLeft(r, p, s);
31     // 如果 ToLeft Test 结果都一样则说明在此三角形中
32     return (pqLeft == qrLeft) && (qrLeft == rpLeft);
33 }
34 
35 bool ToLeft(Point p, Point q, Point s) {
36     // 将 ToLeft Test 等价为求2倍“方向面积”，避免除法或三角运算带来的精度问题
37     return Area2(p, q, s) > 0;
38 }
39 
40 int Area2(Point p, Point q, Point s) {
41     // 2倍“有向面积”
42     return
43         p.x * q.y - p.y * q.x
44         + q.x * s.y - q.y * s.x
45         + s.x * p.y - s.y * p.x;
46 }

虽然我们看似已经把这个问题给解决了，但是熟悉算法复杂度分析的朋友可能一开始就意识到了，这可是个时间复杂度高达 O(n^4) 的算法，在实际应用中这完全是无法被接受的！

所以，在本系列文章的接下来的几篇中，我们将一步一步地，介绍复杂度越来越低的更多计算几何经典算法。同时我们也会简要分析一下凸包问题的理论时间复杂度下界，并且你会看到，的确存在可以达到理论时间复杂度下界的算法。

 

【To Be Continued】