凸包を見つけるためのクイックハルアルゴリズム

定義が言うように、ある集合の凸包 最小の凸集合です 多くを含む 。 ペアごとに異なる点の有限集合の凸包は多面体です。
Quickhullアルゴリズムの1次元のケースを実装するには、関数std :: minmax_elementが適しています。 ネット上では、フラットケースのQuickhullアルゴリズムの多くの実装を見つけることができます。 ただし、任意のディメンションの場合、 qhull.orgからの1つのヘビーウェイト実装のみがすぐに見つかります。

用語の横の括弧内に、英語への翻訳を示します（Tuґboのスタイルで、すみません）。 これは、用語に精通しておらず、リンクで英語のテキストを読んだり、ソースコードを扱うことを決めたりする人に役立ちます。
上記の「標準的な」実装はCで記述され（C ++バインダーがあります）、長い間（GNU OctaveおよびMATLABパッケージの一部として）広く使用されているため、十分にテストされています。 ただし、1つのQuickhullアルゴリズムのみに関連するコードを分離することにした場合、問題が発生します。 最も単純なアプリケーションの場合、もっとコンパクトなものでできるようにしたいと思います。 少し試してから、実装を作成しました（ https://github.com/tomilov/quickhull ）：これは依存関係のない唯一のクラスです。 合計約750行。

幾何学的な概念。

多次元エンティティの操作を開始すると、通常の2次元および3次元のオブジェクトとの類似性が有効であることがわかります。 知識をいくらか合理化する必要があります。いくつかの新しい用語を導入し、それらの関係を明確にします。
プログラマーにとって、幾何学的オブジェクトは主にデータ構造です。 最初に、幾何学的定義（厳密ではありません）を示します。 次に、データを格納し、それらを効果的に操作するのに便利な対応するデータ構造の説明を行います。
空間の次元を 。

• シンプレックス（英語シンプレックス）が設定されている アフィン独立点。 親和性の独立性についてさらに説明します。 これらのポイントは、頂点（Eng.units頂点、pl。頂点）と呼ばれます。
• 多面体（多面体、多面体、多面体、英語の多面体、多面体）は最小値として定義されます アフィン独立点（頂点）; シンプレックス（英語の単体ポリトープ）-最も単純なケース 次元のボディ、より少ない頂点を持つ多面体は必然的にゼロになります 次元のボリューム。
• 平行四辺形（英語の平行四辺形）-平らな平行四辺形と平行六面体の一般化。 シンプレックスの場合、構築できます 対応する平行四辺形（それらはすべて同じ体積を持ちます シンプレックスボリューム）、平行四辺形のジェネレーター（ベクトル）として、1つの固定頂点から残りのベクトルに出現するベクトルを取る トップス。
• 凸多面体の概念（Eng。Convex polytope、convex polyhedron）ここでは説明しませんが、それについての直観的な考えはうまくいきます。 多面体の特殊なケースとしてのシンプレックスは、常に凸です。
• 飛行機の概念も説明しません。 囲まれているスペースよりも次元が1つ少ないことに気づくでしょう。
• 単体ファセットが定義されています アフィン独立点（頂点）。 さらに、単体オブジェクト（凸多面体を除く）についてのみ説明するため、「単体」の定義は省略します。 続編の多くのステートメントは、単純な幾何学的オブジェクトの非縮退（すべてのポイントがペアワイズで異なる）の場合にのみ有効です。
• エッジ（英語の尾根）は2つの面の交差として定義され、 トップ。 2つの面に共通のエッジは1つしかありません。 したがって、1つの顔に rib骨を通して隣人。
• ピークが定義されています ドット。 ここでは、ピークと頂点の概念が一致しないため、3次元空間との直感的な類似性が浮かびますが、このようなオブジェクトでは操作しません。 顔は、ピークを介して任意の数の隣接する顔を持つことができます。したがって、メモリと変換中のデータ構造の関連性を維持するコストの両方で、ピークを介した顔の隣接グラフを保存するのはコストがかかると結論付けることができます。

多面体の表面、その境界、その境界の境界など 英語では顔と呼ばれます。 良いリソースへのリンクを提供します。ここでは、 次元空間、可視または観測可能なオブジェクトは 次元オブジェクト、つまり境界 次元オブジェクト。 英語のファセットは正確に観察された顔です。 直線は常に1次元のエンティティです。 このグラデーションのポイント（頂点）はゼロ次元と見なすことができます。
いくつかのセットを扱う ポイント、我々は対応するセットの検討に進むことができます 他のすべてのポイントから1つのポイントを減算することにより、ベクトルを作成します。 したがって、この点を原点（Eng。Origin）として割り当てます。 すべてのポイントが同じ平面にある場合、平面は原点を通過し始めるようにシフトされます。 アフィニティ独立 線形独立性が実行されるときにポイントが実行されます 対応するベクトルの。 ベクトルの線形独立性の定義は、よく知っていることでしょう。 だから -次元空間はこれ以上選択できない それらはすべてアフィン的に独立しています。 これを説明しましょう。3次元空間では、三角形（3次元シンプレックスの面-四面体）は3つの点（もちろん、1つの直線上にない）によって定義されます。 単一の平面がこれらの3点を通過します。 三角形の3つの頂点のセットに追加されたこの平面上の任意のポイントは、3つのアフィン独立ポイントの可能性を4つのアフィン依存セットに変換します。 他の次元についても同様です。
とりわけ、ハイパーボリューム（英語のハイパーボリューム）の概念を導入する必要があります。

アルゴリズム

任意の次元の場合のQuickhullアルゴリズム自体は、 Barber、C. Bradfordによって提案されました。 ドブキン、デビッドP。; Huhdanpaa、Hannu（1996年12月1日）。 「凸包のクイックハルアルゴリズム。」 数学ソフトウェア上のACMトランザクション22（4）：469–483 著者の「標準的な」C実装は、すでに言及したサイトhttp://qhull.org/にあります。これは、C ++インターフェイスhttps://gitorious.org/qhull/qhullを備えたリポジトリです。 元のソースからアルゴリズムを引用します。

 d + 1ポイントのシンプレックスを作成します
各ファセットF
	未割り当ての各点p
		 pがFを超える場合
			 pをFの外部セットに割り当てます
空でない外部セットを持つ各ファセットFについて
	セットの外側のFの最も遠い点pを選択する
	可視セットVをFに初期化します
	 Vのファセットのすべての未訪問のネイバー
		 pがNを超える場合
			 NをVに追加
	地平線の尾根の集合HはVの境界です
	 Hの各リッジR
		 RとPから新しいファセットを作成します
		新しいファセットをその近隣にリンクします
	新しいファセットごとにF '
		 Vのファセットの外部セット内の各未割り当て点q
			 qがF 'より大きい場合
				 qをF 'の外部セットに割り当てる
	 Vのファセットを削除します

凸包は面のリストです。

シンプレックスの開始。

ご覧のとおり、 開始シンプレックスから開始する必要があります。 どれでも選択できます アフィン独立したポイントですが、何らかの安価なヒューリスティックを使用することをお勧めします。 このステップの主な要件は、シンプレックスが可能な限り「狭く」なることがないようにすることです。 これは、浮動小数点演算が使用される場合に重要です。 元の実装では、最大および最小の座標を持つポイントが使用されます（明らかに、それらは凸包の頂点のセットにも含まれます）。 開始シンプレックスから開始して、後続の各ステップで、より大きな多面体を取得します。その不変のプロパティは常にその凸性になります。 一時的な多面体と呼びます。


次に、顔のいわゆる外部ポイントのリスト （英語の外部セット）に残りのポイントを配布する必要があります-これらは、各顔に対して維持され、まだ考慮されていない（以下-未割り当て、英語未割り当て）ポイントのリストですシンプレックスの内部（または将来の一時的な多面体）の反対側の面に直面します。 ここでは、凸面多面体（特にシンプレックス）のプロパティを使用します。これは、空間全体を各面の平面で割ることによって得られる2つの半空間のうちの1つのみに存在します。 このプロパティは、ある点が凸包の面の外部点のリストのいずれにも該当しない場合、その点の内部に確実に配置されるという事実を意味します。 将来的には、アルゴリズムは、外部ポイントのリストに含まれるポイントのどれが顔の平面から最も離れているかについての情報を必要とします。したがって、外部ポイントのリストにポイントを追加するには、どれが最も遠いかに関する情報を保存する必要があります。 任意のポイントは、外部ポイントの1つのリストのみに分類されます。 漸近的な複雑さは、ポイントがリスト全体にどの程度正確に分布しているかには影響されません。
この段階で、顔を設定する方法、顔までの距離を計算する方法、顔に対する点の向きを決定する方法について説明します。


これで、面の平面を説明することができますが、あいまいさが1つあります。開始シンプレックスの面の頂点のリストで、ポイントの列挙の順序を選択する方法です。著者は、頂点のリストを順序付けしません。そうしません。任意であるが固定された次数の頂点に対して平面の方程式を設定すると、この平面に対する内部ポイント（英語の内部ポイント）の方向がわかります。つまり、指向距離が計算されます。距離が負の場合、超平面の法線の座標と原点からのオフセットで符号が変わります。または、フラグの値を設定して、顔が上下逆になっていることを示します（英語を反転）。凸多面体の内部点は、少なくともその頂点（すべての頂点を取得した場合、いわゆる重心）の算術平均になります。


面ごとに、頂点のリスト、外部ポイントのリストを指定することに加えて、隣接する（英語の隣接リスト）面のリストを指定する必要があります。顔ごとに正確に含まれます 要素。開始シンプレックスには面が含まれているため、その面のいずれについても、他のすべての面（ピース）が隣接する面のリストに含まれることは明らかです。

メインループ。

これで、メインループで作業できる一時的なポリゴンができました。メインループでは、この一時的なポリゴンが完成し、多面体の外側に1つのポイントがなくなる瞬間が来るまで、ますます多くのスペースと元のセットのポイントを「キャプチャ」します。次に、このプロセスについて詳しく説明します。つまり、それが私の実装でどのように機能するかです。
メインループはフェイスループです。各反復で、外部ポイントの空でないリストを持つ面の中から、最良の面が選択されます（英語の最良の面）。外部ポイントのリストから最も遠いポイントが他のフェースよりもフェースから遠いという意味で最高です。外部ポイントの空でないリストを持つ面がなくなると、サイクルが終了します。
メインサイクルのボディでは、選択した面は凸包の意図的に非面として削除されます。しかし、最初に、彼女は「コンポーネントを理解します」。外部ポイントのリストから最も遠いポイントが抽出されます-おそらく凸包の上部です。残りのポイントは、割り当てられていないポイントの一時的なリストに移動されます。 。さらに、Facesのリストから始まる隣接するFacesのリストでは、一時的な多面体のFacesの隣接グラフがバイパスされ、各Facesが既に訪れた（Visit）Facesのリストに追加され、ポイントからの可視性がテストされます 。面がポイントから見えない場合、その隣接する面はそれ以上進みません。隣接する面の中で次の可視面に不可視ポイントがある場合、境界の可視面のリストになります（水平線の前ではなく、水平線の前-「水平線を超えて」）。ただし、すべての隣接する顔が表示されている場合、そのような顔は削除用の顔のリストに追加されます。面の隣接グラフを走査した後、境界に表示される面のリストと削除する面のリストには、ポイントから見えるすべてが含まれますファセット。これで、リストからすべての顔が削除されました。さらに、削除された各顔は、その隣接する隣接顔のリストから削除されません。これは安全です。これは、隣接する面が将来的に隣接する面のリストから見える境界であっても、ポイントから見えないエッジへの近接に関する情報のみが使用されるためです。削除する面のリストに続いて、境界に表示される面のリストがスキャンされます。境界に表示される各面は「解析」され（その後、隣接する面のリストと頂点のリストが必要になります）、削除されます。境界の可視面ごとに、新しい面が作成されます（1から新しい顔のリストに追加されます。この場合、隣接する顔のリストがスキャンされ、このリストから見えない各顔について、特定の可視境界を持つ共通のエッジが検索されます。非表示エッジと境界表示エッジに共通の頂点（つまり、共通エッジの頂点）は、さらに、境界表示面の頂点のリストに書き込まれた順序で、新しい面の頂点のリストに書き込まれます。ちなみに、エッジ（Eng。Horizo​​n ridge）は、ポイントから見える一時的な多面体の表面のその部分の地平線の一部です 。 このエッジに属さない境界可視面の唯一の頂点がポイントに置き換えられます 新しい面の頂点のリストでは、残りの頂点は、表示されている面の境界の頂点の（順序付けられた）リストと同じで、同じ場所にあります。したがって、構築された新しい面の頂点のリストは正しい向きを保証します。つまり、正規化された平面方程式を計算するとき、法線は一時的な多面体の外側に向けられ、オフセット値は正しい符号を取得します。すべての境界表示面が処理された後、新しい面のリストが作成されます。これらの面の超平面方程式は、それらの頂点のリストを作成した直後に計算できます。また、新しい面を作成する場合、隣接する面のリストに1つの面、つまり新しい面の頂点のリストを作成するときに境界の可視面に隣接していた対応する非表示面を追加できます。明らかに残り残りの新しい顔の中から隣接する顔を探す必要があります。
判明したように（Google Performance Toolsからプロファイラーでサンプリングするとき）、アルゴリズムは、大きな次元および/または悪い入力データの場合（つまり、ほとんどのポイントが凸包の頂点である場合）に新しい顔の隣接する顔の検索に費やします。私の実装では、最初はアルゴリズムは非常に簡単でした：新しい面のリストからの面のすべての可能なペア（2つのネストされたループ）がソートされ、1つの不一致を持つ頂点リストの順序（ポイントの入力リストからの番号による）の比較が実行されました（修正されたアルゴリズムstd :: set_differenceまたはstd :: set_symmetric_difference）。それは 比較、ここでは新しい顔の数です。しかし、将来、このボトルネックで順序付けられた連想配列（std :: set）を使用して、隣接するものを検索する速度で（バイナリ=）のペアを獲得することができました。一般に、このタスクに関連する知識の分野全体があります（隣人の決定）-これは局所性に敏感なハッシュです。元のqhull実装でハッシュのみ（LSHではなく）が使用されます。その中には、各エッジ（2つの面）に個別のデータ構造（エッジのエッジのリストと、エッジの「上」と「下」に関する情報）が保存され、各フェイスにはエッジのリストが含まれます。面のすべてのエッジの面（新しい面のリストからの面）の近傍を決定するために、それらの頂点の配列のハッシュの関連付けられた順序なし配列が作成され（1回スキップして）、1つの破棄された頂点 ） それは 節約ハッシュ。ハッシュの衝突（およびそれに続くリストの完全な一致）がある場合、（明らかに）これらの2つの面には共通のエッジがあり（ハッシュは、ハッシュテーブルを削減するために削除できます）、それらのそれぞれは、他の隣接する面のリストに追加できます。標準ライブラリもSTLもハッシュコンバイナなどのツールを持っていないため、この方法は使用しません（更新。最良の選択はXORのみであるという情報を受け取りました）。私自身は、std :: vector <std :: size_t>のstd :: hashを、正当性に関する理論的な情報の正当化なしに特化したくありませんでした。それが私が連想順序付けを使用した理由です隣接する面を検索するためのリブコンテナー。その結果（明らかに）、アルゴリズムの私の実装は、速度でqhullを漸近的に失う 前に 回（さまざまな入力データ用）。

比較テスト。

例は、パッケージマネージャー（Ubuntu 14.04 LTS）を介してインストールされたQuickhull qhullアルゴリズムの「標準的な」実装です。 同じqhull-binパッケージから入力を生成するために、 rboxユーティリティが使用されました 。 rbox tn D4 30000 s> /tmp/qhuckhull.txtコマンドは、4次元球体上の30,000ポイントの座標を持つファイルを作成します。 rbox tn D10 30 s> /tmp/quickhull.txtコマンドは、10次元の球体上の30点の座標を持つファイルを作成します。 プログラムが消費するメモリ量は、 -vスイッチを使用した/ usr / bin / timeユーティリティの出力で確認できます。 / usr / bin / time -v bin / quickhull /tmp/quickhull.txtの出力| head -7は、私の実装のメモリとプロセッサ時間の両方の消費量（ファイルの読み取りと端末への出力のクリア）、およびoutput / usr / bin / time -v qconvex s Qt Tv TI /tmp/quickhull.txt- 「標準的な」 qhull実装用。
凸包の面の数の一致により、実装の正確さを判断します。 しかし、（2次元および3次元モードで）デバッグするために、段階的に（ 一時停止コマンドを使用して） gnuplot形式のアニメーション出力を実装しました 。 コミットのどこかです。 プログラムの出力は、凸包とgnuplot形式で表される入力ポイントの番号付きセットです。
さらに、私はrboxの類似物であるrandomboxユーティリティを何とか書き始めました（しかし、完成しませんでした）。 randomboxは 、均一に分布していないポイントを生成するrboxとは異なり、空間に均一に分布しているポイント（英語の均一な空間分布）を生成するユーティリティとして考案されました。 randomboxは、単一の「ダイヤモンド形」多面体の表面に、単位球（英語の単位球）、ボール（英語の球）、シンプレックス（ポイントのある空間に埋め込むことができる任意の次元）で囲まれたポイントのセットを生成できます。ユニットダイヤモンド）、ユニットキューブ（Eng。ユニットキューブ）内、平行四辺形内、標準ユニットシンプレックス（Eng。ユニットシンプレックス）内、およびコーン（Eng。コーン）またはシリンダー（Eng。シリンダー）のボリュームに多くのポイントを投影します） （空間に囲まれた）任意のシンプレックス内に均一に分布する座標を生成するために、インターネット上で（Dirichlet分布を介して）数値的に安定したアルゴリズムを見つけました。 その結果、最初のものを選択しました。
プロジェクトルートからgnuplotを使用して、凸包を「感じ」させるには、コマンドrbox n D3 40 t |を入力します。 bin / quickhull> /tmp/quickhull.plt && gnuplot -p -e "load '/tmp/quickhull.plt'" 。 rbox n D3 40 tは 、境界キューブ内に40ポイントを生成します （英語の境界ボックス）。 キーtは、PRNGの初期値として現在の時間（秒単位）の使用を指定します。 球上の点を取得するには、キーsを追加する必要があります。 非単純な面がどのように壊れているかを見るのも興味深いです： rbox n D3 c -cube、 rbox n D3 729 M1,0,1- ノード全体にあるポイント。 らせんの凸包は美しく見えます（英語のらせん）： rbox n D3 50 l 。 面と同一平面上にある多数のポイントを同時に配置する場合、アルゴリズムオブジェクトの構築時に設定されるeps定数が、結果の正確性に大きな役割を果たし始めます。 たとえば、 std :: numeric_limits <value_type> :: epsilon（）を小さな値として使用する場合、64個の整数ポイントrbox n D3 64 M3,4の 「回転した」立方体の場合、アルゴリズムは正しくない幾何学的オブジェクトになります。凸多角形であること。

結果または遅延TL; DR。

私の仕事の結果は、C ++でのQuickhullアルゴリズムの実装でした。これは非常にコンパクトで、 qhull.orgでの実装よりもそれほど遅くありません。 このアルゴリズムは、STLで慣例となっているように、入力で空間の次元の値、小さな定数epsの値、およびイテレータのペアで指定された範囲として表されるポイントのセットを受け取ります。 最初の段階で、 create_simplexは開始シンプレックスを構築し、入力ポイントを含むアフィン（サブ）空間のポイントベースを返します。 基底の点の数が点を含むユークリッド空間の次元よりも大きい場合、凸包を完成させるアルゴリズムを実行する必要があります。 出力では、アルゴリズムは入力セットの凸包の面を記述するデータ構造の配列を提供します。これが答えです。 失敗した場合、すべてのポイントを含む低次元の部分空間のポイントベースがあります。 Householderアルゴリズムを使用すると、ポイントの最高座標を無効にしながら、何らかの方法で入力セットを回転させることができます。 そのような座標は破棄でき、Quickhullアルゴリズムは低次元の空間に既に適用できます。

用途

このアルゴリズムには多くの用途があります。 凸包自体を見つけることに加えて、凸包、 ドロネー三角形分割とボロノイ線図の間に関連があるという事実のため、「間接」アルゴリズムの適用を見つけることは難しくありません。
この実装はすでにいくつかの人にとって便利になっており 、感謝の気持ちで（謝辞で）支払われました。これは非常に素晴らしいことです。

アップデート。

2015年6月24日： Kurt Mehlhorn、StefanNäher、Thomas Schilz、Stefan Schirra、Michael Seel、Raimund Seidel、およびChristian Uhrigのアルゴリズムに従って、結果の幾何学的構造の凸性の検証を追加しました。 幾何学的プログラムのチェックまたは幾何学的構造の検証。 Procで。 12周年。 ACMシンポジウム。 計算。 Geom。、ページ159〜165、1996年 。
2015年6月24日：結果の多面体の面には、対応する頂点が反対側（英語の反対側）になるように配置された隣接する面の頂点とインデックスのセットが含まれるようになりました。 これにより、構造が使用しやすくなる場合があります。
06/24/2015：ハッシュを使用して、新しく受け取った顔の近傍を決定します。