【ドロネー三角形分割】GISを作るデータ構造とアルゴリズム【計算幾何学】

概要

前回の記事では, GISの基本機能のひとつ, ボロノイ分割を説明した.

今回は, ドロネー三角形分割について説明する. 地理空間情報では等高線の３Dモデルの作成, コンピューターグラフィックスでは有限点のメッシュ生成, などが利用例である.

ドロネー三角分割は, 有限個の点をもとに三角形の面を作成する, 三角形の領域に分割する処理である. (通常の）三角分割に加え, 綺麗な三角分割をつくる, そんな目的の分割方法である.

今回も「コンピュータ・ジオメトリ　第3版　計算幾何学：アルゴリズムと応用」や, レンセラー工科大学の講義資料を参考に説明する.

ドロネー三角分割

点集合$P$にはいくつかの分割方法が存在する.どの分割方法でも良いかと言えばそうではない. では, 望ましい分割とはなにか？まず目的物をはっきりしよう.

角度最適

ここで$T$を$P$の三角形分割として, そこに$m$個の三角形が含まれるとする. $m$個の三角形には計$3m$個の内角が存在し,それぞれを$\alpha_i (i=1,..,.3m)$とする. さらに,$3m$個の内角を大小で昇順にソートした列を取得し, $A(T):= (\alpha_1,\alpha_2,...,\alpha_{3m})$のように書く. これを$T$の角度ベクトルと呼ぶ.

次に, $P$の別の三角形分割として$T'$を考える. 同様に,$3m$個の内角を$\alpha_i' (i=1,..,.3m)$, そして$A(T'):= (\alpha_1',\alpha_2',...,\alpha_{3m}')$を求める.

ここで, 分割方法に序列を導入したい. 具体的には, $A(T)$が辞書式順序で $A(T')$より後にある, ことを$A(T) > A(T')$という関係式で表す.

辞書式順序とはつまり, 次式を満たすようなインデックス$i$(ただし,$1 \leq i \leq 3m$)が存在するとき, $T$の角度ベクトルは$T'$の角度ベクトルより大きい, とする.

$$すべての j < i に対して \alpha_j = \alpha_j' かつ \alpha_i = \alpha_i'$$

そして, $P$の考えられる三角形分割$T'$に対し, $A(T) \geq A(T')$が成り立つとき, その三角形分割は角度最適(angle-optimal) とよぶ.

このような状態の三角形分割を探索したい, というのがドロネー三角分割の狙いである.

正当な三角分割

つぎに,どのような条件を満たすときに三角形分割が角度最適になるのか見ていく. 情報として,ターレスの定理を導入する.

ターレスの定理

$C$を円として, $l$を$C$と点$a,b$で交差する直線とする. $p$と$q$が$C$上にあり, $r$は$C$の内側に,$s$は$C$の外部にあると仮定すると, 以下の式が成り立つ. (高校数学でも習ったと思うが）

$$\angle arb> \angle apb = \angle aqb > \angle asb$$

次に, 以下の例を考える. 四つの点$p_i, p_j,p_k,p_l$を三角形分割した$A(T)$ができている. この四つの点の三角形分割にはもう一つパターンがある. それは, 辺$\bar{p_ip_j}$を, 辺$\bar{p_kp_l}$へと繋ぎ直した三角形分割$A(T')$である.

同じ頂点群に対し, 頂点の結び方を入れ替える方法を辺フリップ(Edge Flip) とよぶ.
つまり, 四つの点$p_i, p_j,p_k,p_l$を三角形分割すると$A(T)$と$A(T')$が考えうるが,さてどちらが求める三角形分割であるか.

角度最適の式を言い換えると, 次のような条件を満たすとき, その三角形分割は不当であると呼ぶ.

$$\min_{1 \leq i \leq 6} \alpha_i < \min_{1 \leq i \leq 6} \alpha_i'$$

式を口語的に読むと, 辺フリップの結果, 三角形分割後の内角のうち最小のものの大きさが局所的に大きなる場合, その元の辺は不当である,ということだ.

この条件に従えば, いくつかの候補がある三角形分割について,良いものを抽出できる.

さらに良いこと, 不当な辺を検出するために, わざわざ内角の大きさを計算で求める必要はない.
便利な方法をとして,

1. 三角形分割を構成する各三角形の3つの頂点を結ぶ外接円を求める.
2. 1のとき, この外接円に3点以外の点が含まれている場合, その三角形には不当な辺を含んでいると判定する.

実際にこの方法で判定できるのかさきほどの図形でみてみる. (上図）　すると, 不当な辺を含む$A(T)$では, それぞれの外接円が４番目の頂点を内部に含んでいるが, 不当な辺を含まない$A(T')$では, それぞれの外接円の内部に別の点が含まれていない.

なぜ外接円の内外判定で不当な辺を検出できるのか

ここで上述したターレスの定理が活用できる. 以下の例を見て欲しい.

三角形分割$A(T)$の三頂点$p_i, p_j, p_k$を結ぶ赤い外接円$C_1$を扱う. ここから辺フリップにより, 辺$\bar{p_ip_j}$を辺$\bar{p_kp_l}$へと繋ぎ直すと内角はどうなるだろうか.

新しい三角形分割での外接円は一旦考えず, 赤い外接円$C_1$を見て欲しい.
点$p_l$はこの赤い外接円$C_1$の内部にある. そのため内角の大小関係は次のようになっている.

$$\alpha_2' (= \angle p_i p_l p_k) > \alpha_5 (= \angle p_i p_j p_k)$$

これを他の内角にも適用すると以下のような関係になる（図中の番号を考慮）

$$\alpha_2' > \alpha_5, \alpha_5' > \alpha_1, \alpha_3' > \alpha_3, \alpha_4' > \alpha_4$$

全体を通して, 元の三角形分割$A(T)$の内角よりも, 辺フリップ後の内角の最小値が大きくなっていることがわかる.

つまり, 外接円の内外判定で検出されたある辺を辺フリップすれば, 内角の最小は大きくなり, より望ましい三角形分割に近づけるのだ. よって, 外接円の内外判定で検出されたある辺は不当であるとみなせる.

アルゴリズム

基本戦術

さて, 上述した辺フリップを駆使すれば, ある初期状態の三角形分割をスタートに, 徐々に角度最適な三角形分割$A(T)$を作っていけそうだ.

簡単な話, とりあえず何かしらの三角形分割を構成後, すべての辺に対し,正当か不当な辺であるかを判定し, 辺フリップを繰り返せばよい. 結果, これ以上辺フリップの余地がなくなればそのとき, 角度最適な三角形分割$A(T)$が得られている.

この方法は実のところ適切に稼働する. そして, 必ず収束するのだ. 詳細な証明は書籍に任せるが, 辺フリップにより角度ベクトルは単調に増大していくが, 上界があるために, どこかで必ず停止せざるを得ないからだ（無限大には発散できない）

この方法はとてもシンプルに実装できる.
が, 問題としては計算効率が非常に悪いということだ.

乱択逐次構成法

実際のところ, ドロネー三角形分割を求める方法はいくつか考えられる. ここでは, 乱択逐次構成法(randomized incremental algorithm)　を用いた実装方法を扱う.

この乱択逐次構成法は, 入力となる点集合$P={p_1,p_2,...,p_n}$をシャッフルし, ランダムに選択した点を既存の三角形分割に加え, 新たな三角形分割を構成するといった逐次的にドロネー三角形分割を管理する, というものだ.

既存の三角形に新たな点$p_r$を加える, その三角形の３頂点と点$p_r$を結んだ辺を加える. あるいは, 既存の三角形のどこかの辺$e$上に新たな点$p_r$が乗った場合, 辺$e$を共有している２つの三角形の対角頂点と$p_r $を結ぶ辺を加える.　図示例は以下のとおりだ.

点$p_r$が既存三角形の内部あるいは辺上に加わると, 周辺の点と点$p_r$を結ぶことで,この時点での三角形分割ができる. しかし, この三角形分割はドロネー三角形分割か,というとそうではない. 既存の三角形分割はドロネー三角形分割だったとしても, 点を追加後はどこかの辺が不当な辺に変わっている可能性がある. そのため,辺フリップを駆使して再度ドロネー三角形分割を構成しないといけない.

例えば先ほどの例だと, 少なくとも追加した点を内包する三角形および四角形の外周辺は不当なものである可能性がある. そのため, 外周円判定を行い, 必要に応じて辺フリップを行う必要がある.

さらに, 辺フリップを行うことで, これらの三角形および四角形の外にあった三角形にまで不当なものにしてしまう可能性がある. そこで今度は不当な辺を含んでいるかもしれない外側の三角形もチェックする... といった再起処理を, 全体で不当な辺がなくなるまで行うのだ.^[1]

以下にこのアルゴリズムの疑似コードを載せておく.また, 処理のイメージを掴みたい読者は後の章を読むと理解の助けになると思う.

---

アルゴリズム: $DelaunayTriangulation(P)$
入力：　平面上の$n$個の点の集合$P$
出力: $P$のドロネー三角形分割

1. $p_0$を辞書式順序で最も高い位置にある$P$の点とする. すなわち, 最大の$y$座標をもつ点の中で最も右にある点である.
2. $p_{-1}$と$p_{-2}$を十分に遠くにある$R^2$の２点で, $P$全体が三角形$p_0p_{-1}p_{-2}$に含まれるように作成する.
3. ただ１つの三角形$p_0p_{-1}p_{-2}$からなる三角形分割として$T$を初期化する.
4. $P \ {p_0}$のランダム順列$p_1,p_2,...,p_n$を求める.
5. for $r \leftarrow$ to $n$
6. ---- do ($p_r$を$T$に挿入する):
7. -------- $p_r$を含む三角形$p_ip_jp_k \in T$を求める.
8. -------- if $p_r$を含む三角形$p_ip_jp_k$の内部にある.
9. ------------- then $p_r$から$p_ip_jp_k$の三頂点への辺を加える. これによって,$p_ip_jp_k$は３つの三角形にさらに分割される.
10. ------------ $LegalizeEdge(p_r, \bar{p_ip_j}, T)$
11. ------------ $LegalizeEdge(p_r, \bar{p_jp_k}, T)$
12. ------------ $LegalizeEdge(p_r, \bar{p_kp_i}, T)$
13. ------- else $p_r$が$p_ip_jp_k$の辺, たとえば$\bar{p_ip_j}$上にある.
14. ------------ $p_r$から$p_k$への辺を加え, さらに$\bar{p_ip_j}$に接するもう１つの三角形の第3の頂点$p_l$と$p_r$を辺で結ぶことによって, $\bar{p_ip_j}$に接するこれら２つの三角形を４つの三角形に分割する.
15. ------------ $LegalizeEdge(p_r, \bar{p_ip_l}, T)$
16. ------------ $LegalizeEdge(p_r, \bar{p_lp_j}, T)$
17. ------------ $LegalizeEdge(p_r, \bar{p_jp_k}, T)$
18. ------------ $LegalizeEdge(p_r, \bar{p_kp_i}, T)$
19. 2点$p_{-1},p_{-2}$と,これらに接続する辺をすべて$T$から取り除く.
20. Return $T$

---

---

関数: $LegalizeEdge(p_r, \bar{p_ip_j}, T)$

1. (挿入しようとしている点が$p_r$で, $\bar{p_ip_j}$はフリップする必要があるかもしれない$T$の辺である.)
2. if $\bar{p_ip_j}$が正当ではない.
3. ---- then $p_ip_jp_k$を$\bar{p_ip_j}$で$p_rp_ip_j$と隣接する三角形とする.
4. ---- ($\bar{p_ip_j}$のフリップ）$\bar{p_ip_j}$を$\bar{p_rp_k}$で置き換える.
5. ---- $LegalizeEdge(p_r,\bar{p_ip_k},T)$
6. ---- $LegalizeEdge(p_r,\bar{p_kp_j},T)$

---

補足点として,

• $DelaunayTriangulation(P)$内の2-4行目では,点$P$の外界となる巨大な三角形を先に作っている. その理由として, 新規の点を既存三角形に加えて分割する再帰的な処理を動かすための初期状態をつくるためである. もちろん, この巨大な三角形は点$P$のドロネー三角形分割では本来ないため, 19行目で排除している.
• $LegalizeEdge(p_r, \bar{p_ip_j}, T)$の2行目の$\bar{p_ip_j}$が正当性判定は, 上述した外接円判定を行う. 外接円判定についてはとても簡単な行列式があるので利用するとよいだろう.

外接円内外判定

以下のような3点$a,b,c$を通る外接円が, 他の点$d$を含むかどうかは以下の行列式で定まる. ただし, 3点$a,b,c$は,この外接円上で反時計回りに並んでいるものとする.

$$InCircle(a,b,c,d)=det \begin{pmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\ b_x & b_y & b_x^2 + b_y^2 & 1\\ c_x & c_y & c_x^2 + c_y^2 & 1\\ d_x & d_y & d_x^2 + d_y^2 & 1 \end{pmatrix} > 0$$

イメージ

乱択逐次構成法によるドロネー三角形分割の処理の流れを確認する. 以下の事例はFelkel先生の講義資料から引用している.

スタートとして,新規の点$p$を加えることを考える.
いま,既存の三角形分割のうち,三角形$abc$の内部に入ったとする.(①） すると, 三角形$abc$の各頂点と点$p$を結んだ辺を作る. これらは正当な辺である. (②）
次に, 三角形$abc$外周の辺のチェックをする（例えば,辺$ab$). 不当な辺となった場合, 辺フリップにより正当な辺へ切り替える（③, 例では, 辺$ab$->辺$pd$).
この辺フリップにより, さらに隣接する三角形の辺に影響がないかチェックする(④）

以降, 全体で不当な可能性のある辺（図では赤い辺）を全てチェックしていく. 最後に, 不当な可能性のある辺がなくなれば, 新規の点$p$を追加後のドロネー三角形分割が構成されている.(⑧）

まとめ

今回は, 「ドロネー三角形分割」 を説明した. 乱択逐次構成法と辺フリップによる三角形分割の修正処理を構成したアルゴリズムを紹介した.

実は, ドロネー三角形分割は以前と紹介したボロノイ分割と密接な関係があり, 言ってしまえばボロノイ分割が求まればドロネー三角形分割も同時に構成できる.
このように, ドロネー三角形分割はほかにも作り方がある. 例えば, Bowyer-Watson アルゴリズムと呼ばれる処理もある. (コード付き）

https://www.gorillasun.de/blog/bowyer-watson-algorithm-for-delaunay-triangulation/

次回はいよいよ 「地理情報の検索」の仕組みに焦点を当てていく.

脚注

1. この再帰処理は必ず収束する.なぜなら辺フリップにより角度ベクトルは大きくなっていくが, それには上界があるためどこかで止まらざるを得ないからだ