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] 聯繫我。