ドロネー三角形分割の期待最速アルゴリズム
本記事はデータ構造とアルゴリズム Advent Calendar 2019 の 19日目の記事です。
18 日目は @Akazawa_Naoki さんの「ハッシュチェーン、それは歴史を抱え込みながら成長していくデータ構造」でした。
20 日目は @flare さんの「ビットコインのデータ構造」です。
本記事では、計算幾何でしばしば扱われるドロネー三角形分割という構造と、それを構築する期待最速アルゴリズムを解説し、その実装を紹介し、検証・実験の結果もお見せします。ただし、本記事ではアルゴリズムの計算量解析は解説しません。また、本記事は幾何用語を知らなくても雰囲気で読めます。
一般の三角形分割は、同アドベントカレンダーの 2 日目の @yurahuna さんによる記事「三角形分割の数え上げとランダムサンプリング」に登場しました。こちらでは、理論的に興味深い問題が扱われています。
参考資料は下記の通りです。
Computational Geometry: Algorithms and Applications さん、PDF で見れる部分があるんですね……（この記事を書いていた 12/18 に知りました）。
実装だけ知りという方はこちらへ。
はじめに
2次元平面上の $n \in \mathbb{N}$ 個の点からなる集合 $P = \{p_1, \ldots, p_n\} \subset \mathbb{R}^2$ があるとします。このとき、$P$ の三角形分割とは、平面に $P$ 中の 2 点を繋ぐ線分を交差がないように加えていったとき、これ以上線分を加えられなくなった時点での図形です。例えば下図のようなもので、図でも示唆されている通り、必ず $P$ の凸包を含みます。
Computational Geometry: Algorithms and Applications から引用
三角形分割を考える動機として、よく例に挙げられるのは地形生成です。三角形分割の構造と各点に割り振られた高さに基づいて 3 次元プロットをすると、凹凸の付いた地形ができあがります。下図にその例を示します。これを応用すると、現実世界の地形や物体を近似した CG が作れそうですね。
Computational Geometry: Algorithms and Applications から引用
ところで、例示した三角形分割はドロネー三角形分割という特殊な構造になっています。$P$ のドロネー三角形分割とは、$P$ の三角形分割であって、すべての三角形の内角を昇順に並べた列が辞書順最大であるものを指します。ぶっちゃけ、肉眼ではわかりませんね。
地形生成において、ドロネー三角形分割を用いると高さの変化がなめらかになるという効果があります。以下の図はその様子を示したものです。辺の中点に着目したとき、左図のドロネー三角形分割の方が周囲の点から見て自然な高さを設定できています。
Computational Geometry: Algorithms and Applications から引用
また、ドロネー三角形分割には次の嬉しい性質があります。
$P$ に対するユークリッド最小全域木は、ドロネー三角形分割に含まれる。
三角形分割に含まれる辺の個数は、凸包上の点の個数を $k$ として $3n - 3 - k$ 個になることが知られています。よって、ドロネー三角形分割があれば、ユークリッド最小全域木を $O(n \log n)$ 時間で求めることができます。もちろん、ドロネー三角形分割の構築も $O(n \log n)$ 時間でなければなりません。
本記事では、点集合 $P$ が与えられたとき、 $P$ のドロネー三角形分割を構築する問題を扱い、期待計算量 $O(n \log n)$ のアルゴリズムを解説します。さらに、実装を紹介し、検証・実験を行った結果もお見せします。
アルゴリズム
まず、アルゴリズムの概要を述べ、重要な操作である「辺フリップ」を解説します。その後、各ステップを愚直になぞったときのボトルネックについて述べ、それを改善するために導入される「ドロネー木」を解説します。今後、三角形分割を面を共有しない三角形の集合として考えます。
概要
アルゴリズムは以下のように記述できます。ここで、各ステップにおける暫定の三角形分割を $T$ とします。
- 三角形分割 $T$ を「$P$ を覆う十分に大きな三角形」で初期化する。
- 点 $p_r \in P$ をランダムに選択する。
- $p_r$ を含む三角形 $\triangle \in T$ を探す。
- 三角形 $\triangle$ を $p_r$ を含む小さな三角形に分割する。
- 「違反辺」がなくなるまで辺フリップを行う。
- $P$ のすべての点が $T$ に追加されるまでステップ 2-5 を繰り返す。
- $P$ の凸包上の辺であって $T$ に張られていないものがあれば追加する。
- ステップ 1 で生成した三角形とそれに付随する辺を削除する。
$P$ を覆う十分に大きな三角形を計算する方法は様々です。私は冒頭で紹介した MIT の講義資料にある方法をオススメしますが、記事の短縮のため解説しません。
ここでは、ステップ 4 のイメージを掴んでもらいます。以下に例を図示します。
Computational Geometry: Algorithms and Applications から引用
左図は三角形 $(p_i, p_j, p_k)$ を点 $p_r$ で分割したときの様子です。3 つの小さな三角形に分かれており、それぞれ $(p_i, p_j, p_r)$ と $(p_j, p_k, p_r)$ および $(p_k, p_i, p_r)$ です。
もしも $p_r$ が辺 $(p_i, p_j)$ の上にある場合は、その辺を持つ隣り合う 2 つの三角形を考えます。それらの三角形の点であって $p_i, p_j$ でないものを $p_k, p_l$ としします。このとき、4 つの小さな三角形が生成され、それぞれ $(p_i, p_k, p_r$) と $(p_i,p_l,p_r)$ と $(p_j,p_k,p_r)$ および $(p_j,p_l,p_r)$ です。この様子は右図に示されています。
ちなみに、以降の解説で MIT の講義資料から採用しなかった処理があったり、上記アルゴリズムのステップ 7 がMITの講義資料にはなかったりしますが、本記事で解説する手順においては簡単な観察から正当性を示すことができます（詳細略）。MIT の講義資料が実装者に不親切だっただけです。
違反辺と辺フリップ
違反辺とは三角形分割がドロネー三角形分割でないことの証拠となる辺のことです。ここでは『Thales の定理』から得られる違反辺の特徴付けを紹介します。
2 つの三角形 $(p_i, p_j, p_r)$ および $(p_i, p_j, p_k)$ があるとき、辺 $(p_i, p_j)$ が違反辺であることと 3 点 $p_i, p_j, p_r$ を通る円が $p_k$ を内包することは等価である。
ドロネー三角形分割の本質は辺フリップにより違反辺をなくすことです。やることはとても簡単で、上記の特徴付けにより辺 $(p_i, p_j)$ が違反辺と判定されたならば、辺 $(p_i, p_j)$ を削除し辺 $(p_k, p_r)$ を追加します。これらの様子を図示すると以下のようになります。
Computational Geometry: Algorithms and Applications から引用
では、ステップ 5 ではどのように辺フリップを適用すればよいでしょうか？実は、$p_r$ を端点に持つ辺は違反辺にならないことが示せます（詳細略）。よって、$p_r$ と向かい合うような辺にだけ違反辺の判定を行い、違反しているならば辺フリップを行えばよいです。
このとき、辺フリップにより新たに $p_r$ と向かい合う辺が発生しえるため、再帰的に辺フリップが行われます。上図で言えば、辺 $(p_i,p_j)$ を辺フリップした後、新たに辺 $(p_i, p_k)$ と $(p_j, p_k)$ について違反辺の判定を行います。疑似コードは以下のようになります。
Computational Geometry: Algorithms and Applications から引用
ボトルネック
このアルゴリズムについて次の 2 つの事実が知られています。
生成される三角形の個数の期待値は $9n+1$ 以下である。
1 点の挿入後に行われる辺フリップの回数の期待値は $O(1)$ である。
このことから、アルゴリズムを愚直になぞった場合、ステップ 5 が実は高速であり、ステップ 3 『 $p_r$ を含む三角形 $\triangle \in T$ を探す。』の計算量 $O(|T|) = O(n)$ が支配的になります。よって、ステップ 3 がアルゴリズムのボトルネックです。
ドロネー木
本記事が目指す期待計算量 $O(n \log n)$ のアルゴリズムのためには、ステップ 3 を改善する必要があります。そこで、ドロネー木と呼ばれるデータ構造が導入されます。
ドロネー木は木ではなく非巡回有向グラフ (DAG) です（半ギレ）。どのような DAG かというと、ステップ 4 およびステップ 5 で削除された三角形と生成された三角形に対して、削除された側を親とし生成された側を子とする親子関係を保存する DAG です。ドロネー木の根にはステップ 1 で生成した三角形が保存されるようにします。
このとき、a. 1 つの三角形を親として 3 つの三角形が生成される場合と、b. 2 つの三角形を親として 4 つの三角形が生成される場合、c. 2 つの三角形を親として 2 つの三角形が生成される場合があります。下図に a, c を含む例を示します（左図が a で右図が c ）。b は c の派生系なので察してください。
Computational Geometry: Algorithms and Applications から引用
左図では、$\triangle_1$ の中に点が追加されたことで、分割が発生し $\triangle_1$ の子として 3 つの三角形が生成されています。右図では、$\triangle_2$ ともうひとつの三角形が違反辺を共有していたため、それらの子として 2 つの三角形 $\triangle_4, \triangle_5$ が生成されています。
ドロネー木があれば、ステップ 3 で点 $p_r$ を包含する三角形 $\triangle \in T$ の探索を高速化することができます。根から探索を開始し、$p_r$ を内包する子の三角形に下っていき、子がない三角形にたどり着いたとき、それが欲しい $\triangle$ です。
これで、ドロネー三角形分割に対する期待計算量 $O(n \log n)$ を達成するアルゴリズムの完成です。やりましたね。詳細な解析は資料をじっくりと読んでください。
実装
基礎的な幾何処理を知っていれば点の内包判定や違反辺判定の実装ができて、BDD などの心得または Trie などへの親しみがあればドロネー木の実装ができるので、一般に実装できます。
私の GitHub で C++ による実装が公開されています。バグを見つけたらこっそり教えてください。
適用例
以下の図は GitHub にも載せている本実装の適用例です。なんかそれっぽいですね。キャッキャッ。
検証
重複なしでランダム生成した点集合におけるユークリッド最小全域木問題で検証しました。
- 座標が $[0,100]$ の範囲で小数点第3位まで
- 点の個数は 1000 で固定し 100 ケース作成
- 完全グラフ経由で得られる厳密解と比べた目的関数値の誤差を確認
結果、100 ケースすべてで誤差 0 であることを確認できました。しかし、これだけでは厳密なドロネー三角形分割であることの検証にはなっていません。少なくともユークリッド最小全域木問題にはある程度安心して適用できそうであることがわかっただけです。
計算時間の評価
$n = 100, 200, \ldots, 102400$ のランダムなインスタンスをそれぞれ 100 個ずつ用意して、それらを処理する平均時間を計測してみました。結果をプロットしたところ、期待計算量 $O(n \log n)$ を達成できていそうな曲線ができあがりました。
とりあえず、そこそこ速いです。もっと速くしたいけど、幾何アルゴリズムの限界だろうか。
苦労話
本アルゴリズムの実装において、実は「ある辺に対してそれを持つ隣り合う 2 つの三角形を記憶しておく」ことが地味に重要です。主に違反辺の判定で登場しますが、実装上は三角形の生成、三角形の削除でも必要になるので、速度のためには手を抜きたくありません。私が実装期間の大部分を割いた（悩みすぎて 5 日くらい虚無をしていた）のはここです。
実装に落とし込むと「key が辺で value は三角形のペア」的なものを管理すればよさそうです。これは $O(n^2)$ の空間計算量を割けば爆速ですが、期待計算量を完全に無駄にするので論外です。もちろん、挿入も検索も十分速いけど定数時間ではないデータ構造も（計算量にはおそらく影響しないですが）できれば使いたくないです。
となると、期待計算量が定数時間な適当なハッシュテーブルに辺をぶち込んで祈ることになりますが、Hack を防ぐことは最低限やるべきです。すなわち、期待計算量が達成されにくい特定の入力を容易には作れないようにしたいです。
そこである日、もう C++ の unordered_map で雑に済ませてええんちゃうかなどとツイッターで色々ぼやきました。すると、unordered_map への Hack の予防策1 が降ってきたので喜んで使わせていただきました。
予防策の内容は、初めてハッシュ関数が呼ばれたときの高精度タイマーの値といったほぼ予測不可能な値をハッシュ値に噛ませるというものです。完璧な対策ではありませんが、ものすごく簡単なので、今後この予防策は積極的に使っていきたいと思います。
おわりに
最後に、書いておきたいことを箇条書きして締めます。
- ドロネー三角形分割の本質は幾何アルゴリズムではない。
- 実際、初級幾何（点、円、内積、外積、ccw）しか使っていない。
- ただただ実装がめんどくさい。
- 記事やプログラムを書いているときに Delaunay の綴を無限に間違えた。
- Delaynay, Delauney, etc.
- 色んなパターンで間違えた。
- 分割統治により最悪計算量が $O(n \log n)$ になるアルゴリズムがあるらしい。
- こっちを実装したほうがよかったのでは？
- これを知ったのは記事公開の前日。。。
- 記事、濃すぎない？？？？？
- まじで実装から記事までプライベートの時間をかけすぎた。
- そろそろポケモンの厳選を再開します。
それでは、さようなら。