Python中的几何计算（从理论到应用） _几何计算

本文概述

有什么大惊小怪的？
为什么这么难？
CCW简介
我深入研究Python的计算几何和编程
实用建议, 应用程序和关注点
总结

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

哇, 听起来很复杂。
哦, 是的, 凸包。

在这篇文章中, 我想对计算几何学有所启发, 首先简要介绍一下该主题, 然后再根据自己的经验提出一些实用的建议(如果你对这个主题有很好的了解, 请跳过)。
有什么大惊小怪的？虽然入门算法课程通常包含凸包计算几何算法, 但计算几何是一个更为丰富的学科, 很少得到普通开发人员/计算机科学家的足够关注(除非你正在制作游戏或类似的东西)。
从理论上讲很有趣……
从理论上讲, 计算几何学中的问题通常非常有趣。答案令人信服；以及达到这些目标的路径各不相同。我认为, 仅凭这些素质, 它就成为值得研究的领域。
例如, 考虑美术馆问题：我们拥有美术馆, 并希望安装安全摄像机来保护我们的艺术品。但是我们的预算很紧张, 因此我们希望使用尽可能少的相机。我们需要多少个相机？
当我们将其转换为计算几何符号时, 画廊的” 平面图” 就是一个简单的多边形。借助一些肘部润滑脂, 我们可以证明, 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)

【Python中的几何计算（从理论到应用）】对人脑来说直观的东西并没有那么容易地翻译成计算机语言。
更抽象地讲(忽略了用代码表示这些东西的需要), 我们在该学科中看到的问题很难在计算几何算法中进行严格化(“ 使严谨” )。在不使用” 点在多边形内, 如果点在多边形内” 这样的重言式语言下, 我们将如何描述多边形中的场景？其中许多属性是如此基础, 如此基础, 以至于很难对其进行具体定义。
在不使用诸如” 如果在多边形内部, 如果它在多边形内部” 之类的重言式语言的情况下, 我们如何描述多边形点场景？
困难, 但并非不可能。例如, 你可以使用以下定义来严格指定多边形中的点：

如果从该点开始的无限光线与奇数个多边形边缘相交(称为奇偶规则), 则该点位于多边形内部。
如果点的绕线数不为零(定义为定义该多边形的曲线绕该点行进的次数), 则该点在多边形内。

除非你对计算几何有一定的经验, 否则这些定义可能不会成为你现有词汇的一部分。也许这象征着计算几何如何推动你以不同的方式思考。
CCW简介既然我们已经对计算几何问题的重要性和难易程度有所了解, 现在该开始动手了。
该主题的骨干是看似强大的原始操作：逆时针或简称” CCW” 。 (我现在警告你：CCW会一次又一次弹出。)
CCW将A, B和C这三个点作为参数并问：这三个点是否构成逆时针旋转(相对于顺时针旋转)？换句话说, A-> B-> C是逆时针方向的角度吗？
例如, 绿色点是CCW, 而红色点不是：

文章图片
为什么CCW很重要
CCW为我们提供了可以构建的原始操作。它为我们提供了一个开始严格解决计算几何问题的地方。
为了让你对其功能有所了解, 我们来看两个示例。
确定凸度第一个：给定多边形, 你能否确定它是否凸出？凸性是无价的属性：知道多边形是凸的通常可以使性能提高几个数量级。举一个具体的例子：有一个相当简单的PIP算法, 对于凸多边形, 该算法以Log(n)的时间运行, 但对于许多凹多边形, 则失败。
从直觉上讲, 这种差距是有道理的：凸形是” 不错的” , 而凹形可以有锐利的边缘凸出和凸出-它们只是不遵循相同的规则。
一种确定凸度的简单(但非显而易见)的计算几何算法是检查连续顶点的每个三元组是否都是CCW。这仅需要几行Python几何代码(假设以逆时针顺序提供点-如果点以顺时针顺序, 则你希望所有三元组都顺时针旋转)：

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的CCW曲线对应的角度小于180o, 这是定义凸度的一种广泛教授的方法。)
线相交作为第二个示例, 请考虑线段相交, 也可以仅使用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 (B.x - A.x) * (C.y - A.y) > (B.y - A.y) * (C.x - A.x)

要了解此定义的来源, 请考虑向量AB和BC。如果取他们的叉积AB x BC, 这将是沿z轴的向量。但是朝哪个方向(即+ z或-z)？事实证明, 如果叉积为正, 则转弯为逆时针方向；反之亦然。否则, 它是顺时针方向。
除非你对线性代数, 右手法则等有足够的了解, 否则此定义似乎并不直观。但这就是我们拥有抽象的原因-当你考虑CCW时, 只需考虑其直观定义而不是其计算即可。该值将立即清除。
我深入研究Python的计算几何和编程在过去的一个月中, 我一直在努力用Python实现几种计算几何算法。在接下来的几节中, 我将继续使用它们, 我将花一点时间描述我的计算几何应用程序, 该应用程序可在GitHub上找到。
注意：我的经验有限。由于我从事这类产品已有几个月而不是几年的时间, 因此请多加盐。就是说, 这几个月我学到了很多东西, 所以希望这些技巧有用。
柯克帕特里克算法
我工作的核心是实现柯克帕特里克(Kirkpatrick)的点定位算法。问题陈述将类似于：给定一个平面细分(平面中一堆非重叠的多边形)和一个点P, 哪个多边形包含P？考虑类固醇上的多边形点-你将拥有一个平面, 而不是单个多边形。
作为一个用例, 考虑一个网页。当用户单击鼠标时, 网页需要尽快找出用户单击的内容。是按钮A吗？是链接B吗？该网页由不重叠的多边形组成, 因此Kirkpatrick的算法将可以很好地为你提供帮助。
虽然我不会深入讨论该算法, 但你可以在此处了解更多信息。
最小边界三角形
作为子任务, 我还实现了O’ Rourke的算法, 用于在线性时间内计算最小的包围/边界三角形(即, 找到包围凸点集的最小三角形)。
注意：由于计算本身是线性时间, 因此计算最小边界三角形不会影响或损害Kirkpatrick算法的渐近性能, 但对美学目的很有用。
实用建议, 应用程序和关注点前面的部分集中讨论了为什么很难严格地推理计算几何。
实际上, 我们必须处理一系列全新的问题。
还记得CCW吗？作为一个不错的话题, 让我们来看看它的另一项高品质：它可以保护我们免受浮点错误的危害。
浮点错误：为什么CCW为王
在我的计算几何学课程中, 一位着名的教授伯纳德·查泽尔(Bernard Chazelle)发表的论文超出了我的想象, 这是一条规则, 即在尝试描述算法或解决方案时, 我们不能提及角度。
一条规则就是我们甚至不能提到角度。为什么？角度是凌乱的-角度是” 脏” 的。
为什么？角度是凌乱的。角度是” 肮脏的” 。当必须计算角度时, 你需要除法或使用一些近似值(例如, 涉及Pi的任何东西)或某些三角函数。
当你必须在代码中计算角度时, 几乎总是会近似。你将获得一些微小的浮点精度, 这在测试相等性时很重要。你可以通过两种不同的方法求解平面中的某个点, 并且当然可以期望p1.x == p2.x和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模块中看到的方法)是将两个数字定义为相等(如果它们之间仅相差很小的值)。在Python中, 你可能具有：

def almostEqual(x, y, EPSILON=1e-5): return abs(x - y) < EPSILONclass 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))

当然, 已经提出了更高级的解决方案。例如, “ 精确几何计算” 学派(本文所述)旨在使程序中的所有决策路径仅取决于某些计算的符号, 而不是其确切的数值, 从而消除了许多相关的问题。到浮点计算。我们接近平等的方法只是摸索, 但在实践中通常就足够了。
CCW为王在更高的层次上, (甚至可以说)问题是, 我们甚至根据角度或点坐标之类的精确计算量来定义解决方案。为何不单独解决这些症状(即以几乎等于的结果冲洗浮点错误), 为什么不解决原因呢？解决方案：不用考虑角度, 而是考虑CCW, 这将有助于抽象化与浮点计算有关的问题。
这是一个具体的示例：假设你在多边形外部有一些凸多边形P, 一个顶点v和一些点u。你如何确定uv线是否在恒定时间内与v高于或低于v相交？
蛮力解决方案(除了是线性时间, 而不是常数)会带来问题, 因为你必须计算一些确切的线相交点。
我见过的一种固定时间方法涉及：

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

幸运的是, 作者使用了上面的几乎相等的技术来平滑浮点错误。
我认为, 最好完全避免浮点错误的问题。如果花几分钟时间在纸上看问题, 则可以找到完全基于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 u.y > verts[v].y else: return u.y < verts[v].y

该解决方案并不是立即可以用肉眼看到的, 而是用计算几何算法的语言表达的：” 同一边” 是该可信赖算法的经典元素。
完成比完美更重要
在计算几何文献中, 看似简单的操作通常涉及大量的巫术。这为你提供了一个选择：你可以按照一些很难定义的高级解决方案来解决困难的问题, 也可以用一些蛮力以简单的方式进行操作。
同样, 我将使用一个示例：从任意多边形中采样随机内部点。换句话说, 我给你一个简单的多边形, 并且你在其中给了我一个随机点(均匀分布在多边形上)。
通常, 内部点是测试所必需的。在这种情况下, 你对生成它们的计算几何算法没有任何特定的运行时要求。快速而又肮脏的解决方案要花大约2分钟的时间才能实现, 方法是在包含多边形的框中选择一个随机点, 然后查看该点本身是否在多边形内。
例如, 我们可能会错过两次, 而仅在第三点上找到有效样本：

文章图片
这是代码：

class Polygon(object): ... def interiorPoint(self): """Returns a random point interior point""" min_x = min([p.x for p in self.points]) max_x = max([p.x for p in self.points]) min_y = min([p.y for p in self.points]) max_y = max([p.y 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 pdef 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

这就是所谓的拒绝抽样：随机抽取分数, 直到满足你的标准为止。尽管可能需要几个样本才能找到符合你标准的点, 但实际上, 对于你的测试套件而言, 差异可以忽略不计。那为什么还要努力工作呢？总结：在需要时不要害怕走肮脏的路线。
顺便说一句：如果你想要一种精确的随机抽样算法, 下面将介绍一个巧妙的算法。要点：

对多边形进行三角剖分(即, 将其分成三角形)。
选择一个与面积成正比的三角形。
从所选三角形内取一个随机点(恒定时间操作)。

请注意, 此算法要求你对多边形进行三角剖分, 这立即在算法上强加了一个不同的运行时绑定, 并且有必要提供一个用于对任意多边形进行三角剖分的库(我将poly2tri与Python绑定一起使用)。

from p2t import CDTclass Triangle(object): ... def area(self): return abs((B.x * A.y - A.x * B.y) + (C.x * B.y - B.x * C.y) + (A.x * C.y - C.x * A.y)) / 2def 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) * Cclass Polygon(object): ... def triangulate(self): # Triangulate poly with hole cdt = CDT(poly.points) triangles = cdt.triangulate()def convert(t): A = Point(t.a.x, t.a.y) B = Point(t.b.x, t.b.y) C = Point(t.c.x, t.c.y) 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上所说, “ 做总比做得好” 。计算几何问题也是如此。
可视化和自动化测试
由于你在计算几何中处理的许多问题都是根据易于可视化的质量或数量来定义的, 因此视觉测试尤其重要-尽管其本身不足。理想的测试套件将结合视觉和随机自动测试。
理想的测试套件将结合视觉和随机自动测试。
同样, 我们以身作则。考虑测试我们的Kirkpatrick算法的实现。第一步, 算法需要用三角形限制给定的多边形, 并对多边形和外部三角形之间的区域进行三角剖分。这是一个视觉示例, 其中实心绿色线定义了初始多边形, 而虚线定义了三角形区域：

文章图片
很难通过代码来确认已正确执行了该三角剖分, 但是肉眼可以立即看到。注意：我强烈建议你使用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绑定, 并且很容易启动和运行。请参见上面的三角剖分方法, 以了解函数调用。
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开发人员, 对计算几何的美感有所了解, 这是一个充满了令人着迷的问题和同样令人着迷的应用程序的主题。
在实践中, 计算几何实现提出了独特的挑战, 这些挑战将促使你锻炼新的, 令人兴奋的问题解决能力。
如果你有兴趣了解更多信息或对我有任何疑问, 可以通过[电子邮件保护]与我联系。