2011/06/03

ProcessingでDelaunay分割(解説篇)

今日は、前回実装した Delaunay 分割のアルゴリズムをわかりやすく解説したいと思う。

まずは、前知識として 『詳解 OpenCV』 の記述をふたたび引用しよう。
  1. 外部三角形を作り、その頂点の 1 つを開始点とする(これにより、必ず外側の点から開始されることになる)。

  2. 内部の点を追加する。その後、すべての三角形の外接円を探索し、追加した点を含むような三角分割を削除する。

  3. 今削除した三角分割の外接円の内部にある、新しい点も含めて、グラフを三角分割し直す。

  4. 追加する点がなくなるまで、ステップ 2 に戻って繰り返す。
これをいかに解釈し、ソースコードに落とし込んでいくか。

計算幾何に興味がある方はもちろん、普段ネット上のソースをコピペして『動きさえすればそれでよい』と思っている方にも、この記事がプログラミングの楽しみを知るきっかけになれば幸いである(偉そうなこと言ってごめんなさいごめんなさい)。



【必要なデータ構造】

以下に示すように、Delaunay 分割では、点の集合を三角形の集合に変えるため、まず前提としてそれらを表すクラスが必要になる。


それらのために、ぼくは Point クラスと Triangle クラスを用意した。

Point クラス

Triangle クラス

いずれも、さして特別な性質を持たない。

しいて言えば、Triangle クラスに「他の三角形と共有点を持つか」の判定を加えた点がチャームポイント。

先にタネを明かすと、hasCommonPoints メソッドは、処理の最後で外部三角形を除去する際の判定に使うものである。

さて。

点も三角形も、何らかのコレクションに格納して管理される事になる。

少なくとも Delaunay 分割を行っている間は、与えられた点集合には何の変更も与えない(つまり要素数が増減する事もない)が、三角形の集合は処理が完了するまで絶えず変動する事になる。

そこで、静的な点集合は、配列と同じ感覚で扱える ListArrayList)で、一方動的な三角形集合は、追加・削除・検索(参照)が高速な HashSet で管理する事にし、それぞれ pointListtriangleSet と名付けた。

これで下ごしらえは完了である。いよいよ処理の内容に踏み込んでいく。



【逐次添加法】

冒頭のアルゴリズムは、『もともと分割されている Delaunay 三角形領域の内部に、新しい点を一つずつ加えていき、その都度分割をやり直す』というものであり、この事から逐次添加法とも呼ばれる。

例を示そう。

以下のように、既に領域が分割されているところに、黄色い点を加えるとする。

このとき、既に分割が完了している各三角形の外接円を視覚化すると、以下のようになる。

これから追加しようとする点は、上図より明らかに複数の三角形の外接円に包含されているため、外接円特性を満たすように、分割し直す必要がある。

以下の図は、分割し直した結果である。ご覧の通り、どの三角形の外接円も、新たに追加した点を包含する事はなくなっている。

これを、追加する点がなくなるまで繰り返せば、Delaunay 分割が完了する。



【外部三角形の作り方】

逐次添加法で Delaunay 分割を行う場合、はじめに与えられた点集合をすべて包含する架空の『外部三角形』を作らなければならない(そうしないと、『既に Delaunay 分割されている領域の内部に点を追加する』事ができないからである)。
外部三角形は大きすぎて見えないので、各頂点から外部三角形の各頂点に線を引いてみた。
では、このような三角形は、どのようにして求めればよいだろうか。

今回であれば、この問題を『スクリーン(表示領域)全体を包含する正三角形を求める』ものと見なして考えても差し支えなかろう(ちなみにこの方法は、一般化可能な理論である)。

スクリーンは以下のような矩形をなしている。そこでまずは、スクリーンとなる矩形領域を求める事から始める。
Processing の場合、この矩形の左上の座標値は (0, 0) で、右下の座標値は (width, height) となる。

ゆえに、中心座標 c は (width / 2, height / 2) である。

次に、その矩形を包含する円を求める。
中心は c に一致し、半径 r はスクリーン矩形の対角線の長さの半分以上であればよい。

r = (width2 + height2)1/2/2 + λ
(ただし、λ ≧ 0

で。

最後に、その円に外接する正三角形を求める。これが、外部三角形である。
ここで、外部三角形の各点の座標を求める。

外部三角形の各角に対して二等分線を引くと、下図のように 6 つの直角三角形に分割できる。

この直角三角形における辺長比は明らかに 1 : 2 : √3 であり(※忘れてしまった良い子のみんなは中学2年生の数学の教科書を押し入れから引っ張り出してこよう)、さらに短辺の長さは円の半径 r に等しい。

従って、外部三角形の頂点 A, B, C の座標値は、円の中心座標 c = (width / 2, height / 2) = (cx, cy)用いて以下のように求められる。

A = (cx - 2r, cy - √3r)
B = (cx + 2r, cy - √3r)
C = (cx, cy + 2r)

この方法で求めた外部三角形は、スクリーン矩形を包含するには少し大きすぎるものだが、本質ではないためこれ以上の最適化は行わない。

DelaunayTriangles.getHugeTriangle メソッド




【三角形の外接円】

さて、新たに点を追加する際、既に分割されているどの Delaunay 三角形のセットに対しても、外接円特性が満たされなければならない。

つまり、もしも新たな点に対して、それを内包するような外接円をもつ Delaunay 三角形が存在する場合、それをセットから 削除する必要があるのだ(削除後、後述する再分割を施す)。

ここでは、新たに追加する点 p が、任意の三角形の外接円の内側に存在するかどうかの判定を行う議論を行う。

まず円を表現する Circle クラスを作る。これは、中心座標と半径を保持するものである。

Circle クラス

外接円の中心座標 c と半径 r が判れば、点 pc の距離 |p - c| が r より大きいか小さいかを判定する事で、p が外接円の内部にあるかどうかを調べる事ができる。

もし r より小さい場合、その外接円を持つ三角形はセットから除外する必要がある事を意味する。

DelaunayTriangles.DelaunayTriangulation メソッド(一部)

次に、3点 (x1, y1), (x2, y2), (x3, y3) で定義される任意の三角形の外接円の求め方について述べる。

ここで、外接円の中心を (x, y) とすると、円の定義より以下が成り立つ。

(x1 - x)2 + (y1 - y)2 = (x2 - x)2 + (y2 - y)2 = (x3 - x)2 + (y3 - y)2

この方程式を xy について解くと、

x = { (y3 - y1)(x22 - x12 + y22 - y12) + (y1 - y2)(x32 - x12 + y32 - y12) } / c
y = { (x1 - x3)(x22 - x12 + y22 - y12) + (x2 - x1)(x32 - x12 + y32 - y12) } / c

となり、中心座標が得られる。ただし

c = 2 { (x2 - x1)(y3 - y1) - (y2 - y1)(x3 - x1) }

である。

また、その半径 r は、内接している三角形の任意の頂点から円の中心までの距離に他ならない。

以上より、三角形を与えて外接円を得るメソッドは以下のように書ける。

DelaunayTriangles.getCircumscribedCirclesOfTriangle メソッド



【再分割】

新たな点の追加に対して外接円特性を満たさない三角形は、三角形セットから削除し、再分割を施す。

削除に関しては、先日紹介したイテレータを使えばよい。

では、再分割はどうすればよいだろう。これがけっこう難しい。

ひとつひとつ考えてみよう。

追加しようとする点を内包するような外接円がひとつだけの場合は簡単である。

『追加する点と、それまでの三角形の各頂点のうち任意の 2 点を結んでできるすべての三角形』を、新たな分割結果として採用すればよい。

たとえば、新たに追加しようとしている点 p が、点 <p1, p2, p3> から成る三角形 T の外接円の内部に入っている場合、セットから T を消去し、代わりに <p, p1, p2>, <p, p2, p3>, <p, p2, p3> をセットに追加する。

もちろん、これで Delaunay 分割の条件である外接円特性は満たされている。

では、以下のように、新たな点を内包する外接円が複数ある場合はどうだろうか。

この場合は、先ほどのように簡単にはいかない。

『追加した点と、それを内包する三角形をの任意の 2 頂点からできる三角形』を新たな領域として見なすと、以下のように、明らかに不適当なものも含まれる。

実際には、赤色で示した三角形は除外すべきであり、これならばどの点も Delaunay 三角形の外接円の内部に入る事はない。

どうやって不適切な三角形を追加候補から除外すればよいだろう。

実は、三角形のセットに対して、前述の通り反復的に外接円特性を判定し、新たな追加候補の三角形を(不適当である可能性のものも含めて)新たな Delaunay 領域としてセットに追加していった場合、Delaunay 領域として不適切な三角形のみが重複する。

従って、新たに分割し直した追加候補の三角形は、とりあえず一時コレクションに格納しておき、重複判定の後に元の三角形セットへ追加すればよい。

……とはいうものの、必ずしもこれが最もスマートな方法とは言えず、ぼく自身、悩んだ末の苦肉の策である。

重複管理ができるデータ構造としては、HashMap を用いるのが定石だ。管理対象となるオブジェクト(ここでは三角形)を key に、重複しているかどうかの判定を value にする。

今回は、valueBoolean のオブジェクトにし、重複していなければ TRUE、重複していたら FALSE を持たせるようにした。

この重複管理用のメソッドは、以下のように書ける。

DelaunayTriangles.addElementToRedundanciesMap メソッド


このメソッドは、以下のように使われている。

DelaunayTriangles.DelaunayTriangulation メソッド(一部)

重複判定ののち、三角形セットに正当な領域を追加する処理は、以下のように書ける。

DelaunayTriangles.DelaunayTriangulation メソッド(一部)


【後始末】

この一連の処理を、すべての点に対して反復的に行った後、一番はじめに追加した外部三角形のいずれかの頂点を持つ Delaunay 三角形を削除すれば完了だ。

DelaunayTriangles.DelaunayTriangulation メソッド(一部)

以上が、Delaunay 分割のアルゴリズムと実装方針である。

ソースコードをすべて見たい人は昨日の記事を読んでね!

というか、ここまで読んでくれてありがとうありがとう! おつかれさま (バンバン

4 件のコメント:

  1. 初めまして、一応フロントエンドを担当していますが、数学が苦手でこの手の解説をずっと探してました。。わかりやすい解説ありがとうございます。

    返信削除
  2. 外部三角形の求め方で躓いていたのですが次に進めそうです。ありがとうございました。

    返信削除
  3. 貴重な記事をありがとうございます!
    非常に分かりやすいです。

    返信削除
  4. やっと大体コードの書き方等まで理解できた〜

    返信削除

ひとことどうぞφ(・ω・,,)