Pythonの計算幾何学：理論から応用まで

私の経験では、人々が計算幾何学を考えるとき、彼らは通常、次の2つのことのいずれかを考えます。

1. うわー、それは複雑に聞こえます。
2. そうそう、凸包。

この投稿では、計算幾何学に光を当てたいと思います。まず、主題の簡単な概要から始めて、私自身の経験に基づいた実践的なアドバイスに移ります（主題をうまく理解している場合はスキップしてください）。

何が大騒ぎですか？

凸包計算幾何学アルゴリズムは通常、入門アルゴリズムコースに含まれていますが、計算幾何学ははるかに豊富な主題であり、平均的な開発者/コンピューター科学者から十分な注目を集めることはめったにありません（ゲームなどを作成している場合を除く）。

理論的に興味をそそる…

理論的な観点から、計算幾何学の質問はしばしば非常に興味深いものです。 答え、説得力のある; そして、彼らが到達する道はさまざまでした。 私の意見では、これらの資質だけでも研究する価値のある分野になっています。

たとえば、アートギャラリーの問題について考えてみます。私たちはアートギャラリーを所有しており、アートワークを保護するために防犯カメラを設置したいと考えています。 しかし、予算が限られているため、できるだけ少ないカメラを使用したいと考えています。 カメラは何台必要ですか？

これを計算による幾何学的表記に変換すると、ギャラリーの「平面図」は単なるポリゴンになります。 また、エルボーグリースを使用すると、 n / 3のカメラが、どんなに乱雑であっても、 n個の頂点のポリゴンに常に十分であることを証明できます。 証明自体は、双対グラフ、いくつかのグラフ理論、三角測量などを使用します。

ここでは、巧妙な証明手法と、それ自体で評価できるほど興味深い結果が見られます。 しかし、理論的な関連性があなたにとって十分でない場合は…

そして重要な実践

前述したように、ゲーム開発は計算幾何学の適用に大きく依存しています（たとえば、衝突検出は、オブジェクトのセットの凸包の計算に依存することがよくあります）。 地理データの保存と計算の実行に使用される地理情報システム（GIS）も同様です。 ロボット工学も（たとえば、可視性と計画の問題のために）。

なんでそんなにタフなの？

かなり単純な計算幾何学の問題を考えてみましょう。点と多角形が与えられた場合、その点は多角形の内側にありますか？ （これは、ポイントインポリゴンまたはPIP問題と呼ばれます。）

PIPは、計算幾何学が（一見）難しい理由を示す素晴らしい仕事をします。 人間の目には、これは難しい質問ではありません。 次の図が表示され、ポイントがポリゴン内にあることがすぐにわかります。

比較的複雑なポリゴンの場合でも、答えは1、2秒以上はわかりません。 しかし、この問題をコンピューターにフィードすると、次のように表示される場合があります。

 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は、引数として3つのポイントA、B、およびCを取り、質問します。これらの3つのポイントは、反時計回りの回転（時計回りの回転に対して）を構成しますか？ 言い換えれば、A-> B-> Cは反時計回りの角度ですか？

たとえば、緑のポイントはCCWですが、赤のポイントはそうではありません。

CCWが重要な理由

CCWは、構築可能な基本的な操作を提供します。 それは私たちに計算幾何学の問題を厳密にしそして解決し始める場所を与えてくれます。

その力を理解するために、2つの例を考えてみましょう。

凸性の決定

1つ目：ポリゴンが与えられた場合、それが凸状であるかどうかを判断できますか？ 凸性は非常に貴重な特性です。ポリゴンが凸状であることを知っていると、パフォーマンスを桁違いに向上させることができます。 具体的な例として、凸多角形ではLog（n）時間で実行されるが、多くの凹多角形では失敗する、かなり単純なPIPアルゴリズムがあります。

直感的には、このギャップは理にかなっています。凸型の形状は「優れた」ものですが、凹型の形状は鋭いエッジが出入りする可能性があります。同じルールに従わないだけです。

凸性を決定するための単純な（しかし自明ではない）計算幾何学アルゴリズムは、連続する頂点のすべてのトリプレットがCCWであることを確認することです。 これには、数行の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からのCCW曲線は180未満の角度に対応することに注意してください。これは、凸面を定義するための広く教えられている方法です。）

線線交叉

2番目の例として、線分交差について考えます。これは、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の重要性を理解したところで、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を考えてみましょう。 それらの外積ABxBCをとると、これはz軸に沿ったベクトルになります。 しかし、どちらの方向（つまり、+ zまたは-z）ですか？ 結局のところ、外積が正の場合、回転は反時計回りです。 それ以外の場合は、時計回りです。

この定義は、線形代数や右手の法則などを十分に理解していない限り、直感的ではないように見えます。しかし、それが抽象化の理由です。CCWを考えるときは、計算ではなく直感的な定義を考えてください。 値はすぐに明らかになります。

Pythonを使用した計算幾何学とプログラミングへの私の飛び込み

この1か月間、Pythonでいくつかの計算幾何学アルゴリズムの実装に取り​​組んできました。 次のいくつかのセクションでそれらを利用するので、GitHubにある計算幾何学アプリケーションについて説明します。

注：私の経験は確かに限られています。 私はこのようなことに何年もではなく何ヶ月も取り組んできたので、一粒の塩でアドバイスを受けてください。 そうは言っても、私はこの数か月で多くのことを学んだので、これらのヒントが役立つことを願っています。

カークパトリックのアルゴリズム

私の仕事の中核は、ポイント位置のカークパトリックのアルゴリズムの実装でした。 問題の説明は次のようになります。平面の細分割（平面内の重なり合わないポリゴンの束）と点Pが与えられた場合、どのポリゴンにPが含まれますか？ ステロイドのポリゴンの点を考えてみてください。単一のポリゴンではなく、平面一杯のステロイドがあります。

ユースケースとして、Webページを考えてみましょう。 ユーザーがマウスをクリックすると、Webページはユーザーが何をクリックしたかをできるだけ早く把握する必要があります。 ボタンAでしたか？ リンクBでしたか？ Webページは重複しないポリゴンで構成されているため、Kirkpatrickのアルゴリズムが役立つ位置にあります。

アルゴリズムについては詳しく説明しませんが、詳細についてはこちらをご覧ください。

最小境界三角形

サブタスクとして、線形時間で最小の囲み/境界三角形を計算する（つまり、凸状の点のセットを囲む最小の三角形を見つける）ためのO'Rourkeのアルゴリズムも実装しました。

注：最小境界三角形の計算は、計算自体が線形時間であるため、カークパトリックのアルゴリズムの漸近的パフォーマンスを助けたり傷つけたりすることはありませんが、美的目的には役立ちます。

実用的なアドバイス、アプリケーション、および懸念事項

前のセクションでは、計算幾何学を厳密に推論することが難しい理由に焦点を当てました。

実際には、まったく新しい多くの懸念に対処する必要があります。

CCWを覚えていますか？ 素晴らしいセグエとして、その優れた性質のもう1つを見てみましょう。それは、浮動小数点エラーの危険から私たちを保護します。

浮動小数点エラー：CCWが王様である理由

私の計算幾何学のコースでは、私が数えきれないほど多くの論文を発表した尊敬されている教授であるバーナード・チャゼルは、アルゴリズムや解決策を説明しようとするときに角度について言及できないという規則を作りました。

なんで？ 角度が乱雑です。 角度は「汚い」です。 角度を計算する必要がある場合は、除算するか、近似（たとえば、円周率を含むもの）または三角関数を使用する必要があります。

コードで角度を計算する必要がある場合、ほとんどの場合、近似します。 浮動小数点の精度がわずかに低下します。これは、同等性をテストするときに重要です。 2つの異なる方法で平面内のある点を解くことができます。もちろん、 p1.x == p2.x and p1.y == p2.yあると期待できます。 しかし、実際には、このチェックは頻繁に失敗します。 さらに（そして明らかに）、これらのポイントは異なるハッシュを持ちます。

さらに悪いことに、計算を通じて小さな違いが伝播するにつれて、エラーの程度が大きくなります。 （いくつかのより科学的な例については、このペーパーでは、凸包またはドロネー三角形分割を計算するときに何がうまくいかない可能性があるかについて説明します。）

それで、これについて何ができるでしょうか？

ほぼ等しい

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モジュールで見られた）1つのアプローチは、2つの数値が小さい値のイプシロンだけ異なる場合、それらを等しいものとして定義することです。 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未満の差がある2つのポイントを計算しますか。 このタイプのオーバーライドを実装することを強くお勧めします。 同様の方法を行に使用できます。次に例を示します。

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

もちろん、より高度なソリューションが提案されています。 たとえば、「正確な幾何学的計算」の考え方（このペーパーで説明）は、プログラム内のすべての決定パスを、正確な数値ではなく、計算の符号のみに依存させ、関連する懸念の多くを取り除くことを目的としています。浮動小数点計算に。 私たちのほぼ平等なアプローチは表面を傷つけるだけですが、実際には十分であることがよくあります。

CCWは王様です

より高いレベルでは、角度や点座標などの正確な計算量の観点からソリューションを定義することさえ（おそらく）問題があります。 症状だけに対処するのではなく（つまり、ほぼ等しい浮動小数点エラーをalmostEqual ）、原因に対処してみませんか？ 解決策：角度の観点から考えるのではなく、CCWの観点から考えてください。これは、浮動小数点計算に関連する懸念を抽象化するのに役立ちます。

具体的な例を次に示します。凸多角形P 、頂点v 、および多角形の外側の点uがあるとします。 線uvがvの上または下でPと交差するか、まったく交差しないかを一定時間でどのように判断できますか？

いくつかの正確な線交点を計算する必要があるため、ブルートフォースソリューション（一定ではなく線形時間であることに加えて）は問題があります。

私が見た1つの一定時間のアプローチには、次のものが含まれます。

• 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分かかる迅速で汚い解決策は、ポリゴンを含むボックス内のランダムなポイントを選択し、ポイント自体がポリゴン内にあるかどうかを確認することです。

たとえば、2回見逃して、3番目のポイントでのみ有効なサンプルを見つける場合があります。

コードは次のとおりです。

 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で言うように、「完了は完璧よりも優れています」。 同じことが計算幾何学の問題にも当てはまります。

視覚的および自動化されたテスト

計算幾何学で取り組む問題の多くは、視覚化が容易な質または量の観点から定義されているため、視覚的なテストは特に重要ですが、それだけでは不十分です。 理想的なテストスイートには、視覚的テストとランダム化された自動テストの組み合わせがあります。

ここでも、例を挙げて進めます。 Kirkpatrickのアルゴリズムの実装をテストすることを検討してください。 1つのステップで、アルゴリズムは指定されたポリゴンを三角形で境界付け、ポリゴンと外側の三角形の間の領域を三角形分割する必要があります。 これは視覚的な例です。緑色の実線は最初のポリゴンを定義し、破線は三角形分割された領域を定義します。

この三角測量が正しく実行されたことを確認することは、コードで確認するのは非常に困難ですが、人間の目にはすぐにわかります。 注：視覚的なテストを支援するために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：ポリゴンの高速三角測量のための優れたライブラリ。 また、穴のあるポリゴンをサポートします（これは多くの場合重要です）。 C ++で記述されたpoly2triにはPythonバインディングもあり、起動と実行は非常に簡単でした。 関数呼び出しの好みについては、上記のtriangulateメソッドを参照してください。
• scipy.spatial：凸包、ドロネー三角形分割などを計算するための関数が含まれています。 高速（いつものように）、信頼性など。注： toNumpyメソッドで独自のPointデータ型を使用すると便利であることがわかりました。defnp 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]までご連絡ください。