Python 中的计算几何：从理论到应用

当人们想到计算几何时，根据我的经验，他们通常会想到以下两件事之一：

1. 哇，听起来很复杂。
2. 哦，是的，凸包。

在这篇文章中，我想对计算几何进行一些说明，首先是对该主题的简要概述，然后再根据我自己的经验提出一些实用的建议（如果你对这个主题有很好的了解，请跳过）。

有什么大惊小怪的？

虽然凸包计算几何算法通常包含在介绍性算法课程中，但计算几何是一个更丰富的主题，很少得到普通开发人员/计算机科学家的足够关注（除非你正在制作游戏或其他东西）。

理论上很有趣…

从理论的角度来看，计算几何中的问题通常非常有趣。 答案，令人信服； 以及到达它们的路径，各不相同。 在我看来，仅这些品质就使它成为一个值得研究的领域。

例如，考虑美术馆问题：我们拥有一家美术馆，并希望安装安全摄像头来保护我们的艺术品。 但是我们的预算很紧，所以我们希望使用尽可能少的摄像机。 我们需要多少台相机？

当我们将其转换为计算几何符号时，画廊的“平面图”只是一个简单的多边形。 并且用一些肘部油脂，我们可以证明n/3 个相机对于n个顶点上的多边形总是足够的，无论它多么混乱。 证明本身使用对偶图、一些图论、三角剖分等等。

在这里，我们看到了一种巧妙的证明技术和一个足够好奇的结果，可以单独欣赏。 但是，如果理论相关性对您来说还不够……

和重要的实践

正如我之前提到的，游戏开发在很大程度上依赖于计算几何的应用（例如，碰撞检测通常依赖于计算一组对象的凸包）； 地理信息系统 (GIS) 也是如此，用于存储和执行地理数据计算； 和机器人技术（例如，用于可见性和规划问题）。

为什么这么难？

让我们来看一个相当简单的计算几何问题：给定一个点和一个多边形，该点是否位于多边形内部？ （这称为多边形中的点或 PIP 问题。）

PIP 很好地说明了为什么计算几何会（看似）困难。 在人眼看来，这不是一个难题。 我们看到下图，很明显该点在多边形中：

即使对于相对复杂的多边形，答案也不会超过一两秒。 但是当我们将此问题输入计算机时，它可能会看到以下内容：

 poly = Polygon([Point(0, 5), Point(1, 1), Point(3, 0), Point(7, 2), Point(7, 6), Point(2, 7)]) point = Point(5.5, 2.5) poly.contains(point)

对人脑来说直观的东西并不能那么容易地转化为计算机语言。

更抽象地说（并且忽略在代码中表示这些东西的需要），我们在这门学科中看到的问题很难在计算几何算法中严格化（“变得严谨”）。 如果不使用“如果一个点在多边形内部，那么它就在多边形内部”这样的重言式语言，我们将如何描述多边形中的点场景？ 其中许多属性是如此基础和基础，以至于很难具体定义它们。

困难，但并非不可能。 例如，您可以使用以下定义严格化多边形中的点：

• 如果从该点开始的任何无限射线与奇数个多边形边相交，则该点位于多边形内（称为奇偶规则）。
• 如果一个点具有非零绕组数（定义为定义多边形的曲线围绕该点行进的次数），则该点位于多边形内。

除非您对计算几何有一定的经验，否则这些定义可能不会成为您现有词汇表的一部分。 也许这就是计算几何如何推动你以不同方式思考的象征。

介绍 CCW

现在我们已经了解了计算几何问题的重要性和难度，是时候弄湿我们的手了。

该主题的主干是一个看似强大的原始操作：逆时针，或简称为“CCW”。 （我现在警告你：CCW 会一次又一次地弹出。）

CCW 将三个点 A、B 和 C 作为参数并询问：这三个点是否构成逆时针转动（相对于顺时针转动）？ 换句话说，A -> B -> C 是逆时针角度吗？

例如，绿点是 CCW，而红点不是：

为什么 CCW 很重要

CCW 为我们提供了一个可以构建的原始操作。 它为我们提供了一个开始严谨和解决计算几何问题的地方。

为了让您了解它的力量，让我们考虑两个例子。

确定凸度

第一个：给定一个多边形，你能确定它是否是凸的吗？ 凸性是一个无价的属性：知道你的多边形是凸的通常可以让你将性能提高几个数量级。 作为一个具体的例子：有一个相当简单的 PIP 算法，它在 Log(n) 时间内运行凸多边形，但对于许多凹多边形失败。

直观地说，这种差距是有道理的：凸形是“不错的”，而凹形可以有尖锐的边缘突出和突出——它们只是不遵循相同的规则。

用于确定凸性的简单（但非显而易见）计算几何算法是检查连续顶点的每个三元组是否为逆时针方向。 这只需要几行 Python 几何代码（假设这些points是按逆时针顺序提供的——如果points是按顺时针顺序提供的，那么您会希望所有三元组都是顺时针的）：

 class Polygon(object): ... def isConvex(self): for i in range(self.n): # Check every triplet of points A = self.points[i % self.n] B = self.points[(i + 1) % self.n] C = self.points[(i + 2) % self.n] if not ccw(A, B, C): return False return True

用几个例子在纸上试试这个。 你甚至可以使用这个结果来定义凸度。 （为了使事情更直观，请注意从 A -> B -> C 的逆时针曲线对应的角度小于 180，这是一种广泛教授的定义凸性的方法。）

线交点

作为第二个例子，考虑线段相交，也可以单独使用 CCW 解决：

 def intersect(a1, b1, a2, b2): """Returns True if line segments a1b1 and a2b2 intersect.""" return ccw(a1, b1, a2) != ccw(a1, b1, b2) and ccw(a2, b2, a1) != ccw(a2, b2, b1)

为什么会这样？ 线段相交也可以表述为：给定具有端点 A 和 B 的线段，另一段的端点 C 和 D 是否位于 AB 的同一侧？ 换句话说，如果从 A -> B -> C 和 A -> B -> D 的转弯方向相同，则线段不能相交。 当我们使用这种类型的语言时，很明显这样的问题是 CCW 的生计。

严格的定义

现在我们已经了解了 CCW 的重要性，让我们看看它是如何计算的。 给定点 A、B 和 C：

 def ccw(A, B, C): """Tests whether the turn formed by A, B, and C is ccw""" return (Bx - Ax) * (Cy - Ay) > (By - Ay) * (Cx - Ax)

要了解此定义的来源，请考虑向量 AB 和 BC。 如果我们取它们的叉积 AB x BC，这将是一个沿 z 轴的向量。 但在哪个方向（即+z 或-z）？ 事实证明，如果叉积为正，则转为逆时针； 否则，它是顺时针的。

除非你对线性代数、右手定则等有很好的理解，否则这个定义看起来不直观。但这就是我们有抽象的原因——当你想到 CCW 时，只考虑它的直观定义而不是它的计算。 该值将立即清晰。

我潜入计算几何和使用 Python 编程

在过去的一个月里，我一直致力于在 Python 中实现几种计算几何算法。 因为我将在接下来的几节中使用它们，所以我将花一点时间来描述我的计算几何应用程序，可以在 GitHub 上找到。

注意：我的经验当然是有限的。 由于我已经在这些东西上工作了几个月而不是几年，所以请接受我的建议。 也就是说，在那几个月里我学到了很多东西，所以我希望这些技巧对我有用。

柯克帕特里克算法

我工作的核心是实施柯克帕特里克算法的点定位。 问题陈述类似于：给定一个平面细分（平面中的一堆不重叠的多边形）和一个点 P，哪个多边形包含 P？ 想想类固醇上的多边形点——而不是单个多边形，你有一个平面。

作为一个用例，考虑一个网页。 当用户点击鼠标时，网页需要尽快弄清楚用户点击了什么。 是按钮A吗？ 是链接B吗？ 该网页由不重叠的多边形组成，因此柯克帕特里克算法可以很好地提供帮助。

虽然我不会深入讨论该算法，但您可以在此处了解更多信息。

最小边界三角形

作为一个子任务，我还实现了 O'Rourke 算法，用于在线性时间内计算最小包围/边界三角形（即找到包围凸点集的最小三角形）。

注意：计算最小边界三角形不会帮助或损害 Kirkpatrick 算法的渐近性能，因为计算本身是线性时间的，但它对于美学目的很有用。

实用建议、应用和关注点

前面的部分重点讨论了为什么计算几何很难严格推理。

在实践中，我们必须处理一系列全新的问题。

还记得 CCW 吗？ 作为一个不错的选择，让我们看看它的另一个优点：它保护我们免受浮点错误的危险。

浮点错误：为什么 CCW 为王

在我的计算几何课程中，Bernard Chazelle 是一位受人尊敬的教授，他发表的论文数不胜数，他规定我们在尝试描述算法或解决方案时不能提及角度。

为什么？ 角度很乱。 角度是“脏的”。 当你必须计算一个角度时，你需要除法，或者使用一些近似值（例如任何涉及 Pi 的东西）或一些三角函数。

当你必须在代码中计算一个角度时，你几乎总是在逼近。 你会偏离一些微小的浮点精度——这在你测试相等性时很重要。 您可以通过两种不同的方法求解平面中的某个点，当然，期望p1.x == p2.x and p1.y == p2.y 。 但是，实际上，这种检查经常会失败。 此外（而且很明显），这些点将具有不同的哈希值。

更糟糕的是，随着您的微小差异在您的计算中传播，您的错误程度会增加。 （对于一些更科学的例子，本文介绍了计算凸包或 Delaunay 三角剖分时可能出现的问题。）

那么，我们能做些什么呢？

几乎相等

Python 计算几何问题的一部分是，我们需要在一个事物很少精确的世界中精确。 与处理角度相比，这将成为一个更常见的问题。 考虑以下：

 # Define two random points p1 = RandomPoint() p2 = RandomPoint() # Take the line through them l1 = Line(p1, p2) # Shift both points up by sqrt(2) p1.y += sqrt(2) p2.y += sqrt(2) l2 = Line(p1, p2) # Slope 'should' be the same? if abs(l1.slope - l2.slope) > 0: print "Error!" # Error!

事实上，这段代码会打印“错误！” 大约 70% 的时间（根据经验）。 我们可以通过对平等的定义稍微宽容一些来解决这个问题。 也就是说，通过牺牲一定程度的准确性。

我使用过的一种方法（例如在某些 OpenCV 模块中看到的）是将两个数字定义为相等，如果它们仅相差一些小值 epsilon。 在 Python 中，您可能有：

 def almostEqual(x, y, EPSILON=1e-5): return abs(x - y) < EPSILON class Point(object): ... def __eq__(self, that): return (almostEqual(self.x, that.x) and almostEqual(self.y, that.y))

在实践中，这非常有帮助。 很少，如果有的话，你会计算两个相差小于 1e-5 的点实际上是不同的点。 我强烈建议实施这种类型的覆盖。 类似的方法可以用于线条，例如：

 class Line(object): ... def __eq__(self, that): return (almostEqual(self.slope, that.slope) and almostEqual(self.intercept, that.intercept))

当然，已经提出了更先进的解决方案。 例如，“精确几何计算”学派（在本文中描述）旨在让程序中的所有决策路径仅取决于某些计算的符号，而不是其精确数值，从而消除了许多相关的问题到浮点计算。 我们的近似相等方法只是触及表面，但在实践中通常就足够了。

逆时针为王

在更高的层次上，我们甚至根据角度或点坐标等精确计算量来定义我们的解决方案（可以说）是有问题的。 与其单独解决症状（即，用almostEqual浮点错误），为什么不解决原因呢？ 解决方案：不要从角度思考，而是从 CCW的角度思考，这将有助于抽象出与浮点计算相关的问题。

这是一个具体的例子：假设你有一些凸多边形P ，一个顶点v ，以及多边形外部的一些点u 。 您如何确定线uv是否在恒定时间内与P相交于v上方或下方，或者根本不相交？

蛮力解决方案（除了是线性时间，而不是恒定的）将是有问题的，因为您必须计算一些精确的线交点。

我见过的一种恒定时间方法包括：

• 使用arctan2计算一些角度。
• 通过乘以 180/Pi 将这些角度转换为度数。
• 检查这些不同角度之间的关系。

幸运的是，作者使用上面的almostEqual技术来平滑浮点错误。

在我看来，最好完全避免浮点错误的问题。 如果您花几分钟时间在纸上查看问题，您可以获得完全基于 CCW 的解决方案。 直觉：如果与v相邻的顶点在uv的同一侧，则线不相交； 否则，查看u和v是否在相邻顶点之间的直线的同一侧，并根据结果比较它们的高度。

这是用于测试v以上交集的 Python 代码（下面的交集只是反转了比较的方向）：

 def intersectsAbove(verts, v, u): """ Returns True if uv intersects the polygon defined by 'verts' above v. Assumes v is the index of a vertex in 'verts', and u is outside of the polygon. """ n = len(verts) # Test if two adjacent vertices are on same side of line (implies # tangency) if ccw(u, verts[v], verts[(v - 1) % n]) == ccw(u, verts[v], verts[(v + 1) % n]): return False # Test if u and v are on same side of line from adjacent # vertices if ccw(verts[(v - 1) % n], verts[(v + 1) % n], u) == ccw(verts[(v - 1) % n], verts[(v + 1) % n], verts[v]): return uy > verts[v].y else: return uy < verts[v].y

解决方案对肉眼来说并不是很明显，但它采用计算几何算法的语言：“同一侧”是该可靠算法的经典元素。

完成比完美更重要

在计算几何文献中，看似简单的操作往往涉及相当多的魔法。 这给了你一个选择：你可以用困难的方式做事，按照一些论文定义了一个非常先进的解决方案来解决一个不那么先进的问题——或者你可以用一点蛮力用简单的方法做事。

同样，我将使用一个示例：从任意多边形中采样一个随机内部点。 换句话说，我给你一个简单的多边形，你给我它里面的一个随机点（均匀分布在多边形上）。

通常，测试需要内部点。 在这种情况下，您对生成它们的计算几何算法没有任何特定的运行时要求（在合理范围内）。 快速而肮脏的解决方案需要大约 2 分钟的时间来实现，它是在包含多边形的框中选择一个随机点，并查看该点本身是否在多边形内。

例如，我们可能会错过两次，只在第三点找到有效样本：

这是代码：

 class Polygon(object): ... def interiorPoint(self): """Returns a random point interior point""" min_x = min([px for p in self.points]) max_x = max([px for p in self.points]) min_y = min([py for p in self.points]) max_y = max([py for p in self.points]) def x(): return min_x + random() * (max_x - min_x) def y(): return min_y + random() * (max_y - min_y) p = Point(x(), y()) while not self.contains(p): p = Point(x(), y()) return p def contains(self, p): for i in range(self.n): p1 = self.points[i] p2 = self.points[(i + 1) % self.n] p3 = self.points[(i + 2) % self.n] if not ccw(p1, p2, p3): return False return True

这被称为拒绝抽样：随机抽取一个点，直到满足您的标准。 虽然可能需要多个样本才能找到符合您标准的点，但实际上，您的测试套件的差异可以忽略不计。 那为什么还要更努力呢？ 总结：当场合需要时，不要害怕走肮脏的路线。

顺便说一句：如果你想要一个精确的随机抽样算法，我在下面实现了一个聪明的算法。 它的要点：

1. 对多边形进行三角剖分（即，将其分成三角形）。
2. 选择一个概率与其面积成正比的三角形。
3. 从所选三角形中取一个随机点（恒定时间操作）。

请注意，此算法需要您对多边形进行三角剖分，这会立即对算法施加不同的运行时限制，并且您必须有一个用于对任意多边形进行三角剖分的库（我使用带有 Python 绑定的 poly2tri）。

 from p2t import CDT class Triangle(object): ... def area(self): return abs((Bx * Ay - Ax * By) + (Cx * By - Bx * Cy) + (Ax * Cy - Cx * Ay)) / 2 def interiorPoint(self): r1 = random() r2 = random() # From http://www.cs.princeton.edu/~funk/tog02.pdf return (1 - sqrt(r1)) * A + sqrt(r1) * (1 - r2) * B + r2 * sqrt(r1) * C class Polygon(object): ... def triangulate(self): # Triangulate poly with hole cdt = CDT(poly.points) triangles = cdt.triangulate() def convert(t): A = Point(tax, tay) B = Point(tbx, tby) C = Point(tcx, tcy) return Triangle(A, B, C) return map(convert, triangles) def interiorPoint(self): # Triangulate polygon triangles = self.triangulate() areas = [t.area() for t in triangles] total = sum(areas) # Calculate normalized areas probabilities = [area / total for area in areas] weighted_triangles = zip(triangles, probabilities) # Sample triangles according to area r = random() count = 0 for (triangle, prob) in weighted_triangles: count += prob # Take random point from chosen triangle if count > r: return triangle.interiorPoint()

希望从代码中可以看出额外的努力。 请记住：正如 Facebook 所说，“完成胜于完美”。 计算几何问题也是如此。

可视化和自动化测试

由于您在计算几何中处理的许多问题都是根据易于可视化的质量或数量来定义的，因此视觉测试尤其重要——尽管它本身还不够。 理想的测试套件将结合视觉和随机自动化测试。

同样，我们继续举例。 考虑测试我们对柯克帕特里克算法的实现。 在一个步骤中，该算法需要用一个三角形来限制给定的多边形，并对多边形和外三角形之间的区域进行三角剖分。 这是一个视觉示例，其中绿色实线定义了初始多边形，虚线定义了三角区域：

确认这个三角测量已经正确执行是很难通过代码来验证的，但是人眼可以立即看出这一点。 注意：我强烈建议使用 Matplotlib 来帮助您进行视觉测试——这里有一个很好的指南。

稍后，我们将要验证算法是否正确定位点。 一种随机、自动化的方法是为每个多边形生成一堆内部点，并确保我们返回所需的多边形。 在代码中：

 class TestLocator(unittest.TestCase): ... def runLocator(self, polygons): # Pre-process regions l = Locator(polygons) # Ensure correctness for polygon in polygons: # Test 100 random interior points per region for k in range(100): target = polygon.interiorPoint() target_polygon = l.locate(target) self.assertEqual(polygon, target_polygon) self.assertTrue(target_polygon.contains(target))

然后，我们可以在不同的多边形集上使用runLocator方法，从而为我们提供了一个多样化的测试套件。

开源解决方案

无论您选择哪种编程语言，计算几何都有一套不错的开源库和解决方案可供使用（尽管 C++ 库似乎出现了不成比例的数量）。

使用现有开源解决方案（如 Python 中的科学计算）的好处是众所周知的，并且已被广泛讨论，因此我不会在这里继续讨论。 但我想我会提到一些我发现有用的以 Python 为中心的资源：

• poly2tri：一个很棒的多边形快速三角剖分库。 还支持（这通常是至关重要的）带有孔的多边形。 poly2tri 用 C++ 编写，还具有 Python 绑定，并且很容易启动和运行。 请参阅上面的triangulate方法以了解函数调用。
• scipy.spatial：包括用于计算凸包、Delaunay 三角剖分等的函数。 快速（一如既往）、可靠等。注意：我发现使用我自己的Point数据类型和toNumpy方法很有用： def np(self): return [self.x, self.y] 。 然后，我可以轻松调用 scipy.spatial 方法，例如： scipy.spatial.ConvexHull(np.array(map(lambda p: p.np()), points)) 。
• OpenCV：开源计算机视觉库有一些不错的独立计算几何模块。 特别是，在自己实现它之前，我使用了它的最小封闭三角形函数一段时间。

结论

我希望这篇文章能让你作为 Python 开发人员领略计算几何的美妙之处，这是一个充满迷人问题和同样迷人应用的学科。

在实践中，计算几何实现提出了独特的挑战，将推动您锻炼新的和令人兴奋的解决问题的技能。

如果您有兴趣了解更多信息或对我有任何疑问，可以通过 [email protected] 联系我。