はじめに

processingでドロネー三角形分割(Delaunay Trianglation)を実装していきます。

Qiitaや他のサイトでは、processingやJava、C#(Unity)、JavaScript等での実装を挙げられています。
それぞれ分割方法の流れや、各種判定方法が違っていたので、自分でも実装してみました。

環境

Processing 3.5.3 (Java mode)
- openProcessingで動くprocessing(Javaモード)は1.6.6までなのが悲しい。

ドロネー三角形分割とは

各サイトによって色々書いてあったり書いてなかったりしますが、ざっくり言うと『平面上に点の集合Pがあるとき、各点をとあるルールに従ってバランスよく繋げ、三角形の集まりにする』のようなイメージです。

n個の点を含む点集合から三角形分割を作成するだけであれば、$\Omega(2.43n)$個の分割パターンがあるそうですが、ドロネー図は『三角形分割を構成する三角形の内角すべての中で、最小角が最大になる三角形分割』を指しています。

難しく書いてありますが、イメージとしては、長細い三角形(=内角が小さい)ものが含まれているとバランス悪く見える。内角を大きくしていると、どの三角形もそれなりにきれいな形に見える。ということです。

こちらに詳しく解説されています。

ちなみに平面上以外にも、3次元空間(3Dスキャナで得られた点群からメッシュを作成する等)や、それ以上の次元でも分割できます。

Wikipedia。 https://ja.wikipedia.org/wiki/%E3%83%89%E3%83%AD%E3%83%8D%E3%83%BC%E5%9B%B3

ドロネー図（Delaunay diagram）あるいはドロネー三角形分割は、距離空間内に離散的に分布した点の集合に対し得られる、それらをある方法に従い辺で結んだ図形である。
MathWorksのMATLABの解説。

Delaunay は、2 次元または 3 次元空間にある点集合のドローネ三角形分割を作成します。2 次元のドローネ三角形分割は、各三角形の外接円が内部に他の点を含まないようにします。この定義は、高次元にも自然に拡張できます。

wikipediaには"ある方法に従い"としか書いていませんが『各三角形の外接円が内部に他の点を含まないように、分割していく』ことが、ドロネー分割での肝となる部分です。

最終的に点群や三角形群を処理する方法はとりあえず置いておいて、まずは点一つ一つや三角形一つ一つに対して判定していく事を考えていきましょう。具体的にこの部分を図とコードで解説していきます。

余談になりますが、MATLAB関数(?)の公式ドキュメントは、数式やグラフ、図なども充実しているので、数学的、統計的な処理を学ぶ場合に参考になることが多い気がしますので、ぜひ参考にしてみてください。

三角形の外接円の内部に点があるかを判定する

シューティングゲーム等で当たり判定の実装などに用いられる方法で考えます。円を表すCircleクラスに、中心座標をcenter、半径をradiusの情報を持たせます。この円の半径radiusと点pからcenterの距離の大小比較をする事で、点pが、円の内部にあるかを判定できます。

processingにはx,y座標上の2点の距離を計算するdistメソッドが存在するのでこれを利用しましょう。
(Qiitaのハイライトが.pdeに対応していないので.javaにしていますが、内容は変わりません。)

Circleクラスサンプルコード

Circle.java

public class Circle
{
    private PVector center;
    private float radius;

    /* 解説のため、コンストラクタなどは省略 */

    // 円の半径radiusと中心centerから点pまでの距離を比較
    // 距離が半径未満の場合、円の内部にあると判定する。
    // 比較条件を<=とすると、円周上の点も含んでしまう為、注意。
    public boolean Contains(PVector p)
    {
        return this.center.dist(p) < this.radius;
    }
}

このままでは、円と点での判定になってしまいます。実際には上述のとおり、三角形の外接円と、点の判定が必要となるので、三角形Triangleクラスに、外接円を取得するメソッドを作成し、三角形(の外接円)と点で判定ができるようにしていきます。

外接円を求める計算式についてはこちら。

https://ja.wikipedia.org/wiki/%E5%A4%96%E6%8E%A5%E5%86%86

Triangleクラスサンプルコード

Triangle.java

public class Triangle
{
    // 各頂点はA,B,Cとする。
    PVector A;
    PVector B;
    PVector C;

    /* 中略 */

    public boolean Contains(PVector p)
    {
        return this.GetCircumscribedCircle().Contains(p);
    }

    // 三角形の外接円を取得するメソッド。
    private Circle GetCircumscribedCircle()
    {
        float x1 = this.A.x;
        float y1 = this.A.y;
        float x2 = this.B.x;
        float y2 = this.B.y;
        float x3 = this.C.x;
        float y3 = this.C.y;

        float c = 2.0 * ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1));
        float x = ((y3 - y1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)+(y1 - y2) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1))/c;
        float y = ((x1 - x3) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)+(x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1))/c;

        PVector center = new PVector(x, y);
        float r = PVector.dist(center, this.A);

        return new Circle(center, r);
    }
}

というわけで、三角形の外接円の内部に点が含まれるかどうかを判定する部分ができました。

では、早速このメソッドを用いて、プログラムの流れを考えてみましょう。

プログラム全体の流れを考える

愚直に、三角形分割のパターンを総当りして判定していく方法

n個の点からドロネー三角形分割を実施していくにあたり、とりあえず、簡単に思いつく方法として、総当たりを考えます。。
プログラムの流れは

n個の点すべてを並べる
n個の点を結んでできる三角形のパターンを1つ作り、各三角形が、外接円の内部に他の点を含んでいないか、調査する。
すべての三角形が条件を満たしていれば、完成。そうでなければ、2.に戻り、三角形のパターンを再度生成。再度チェックする。

三角形の点をどのように選ぶか、という点は難しいかもしれませんが、流れと判定部分についてはわりとシンプルだと思います。

ただ直感的にも分かる通り、点nの数量が増えていくごとに、三角形のパターン数もふえますので、計算量が指数関数的に増加していくこととなります。

※こちらの記事で、計算量について詳しく述べられています。

また、新規にあたらしい点を追加したい！となった場合も、再度総当りをし直すことにもなりそうです。数点の点群であれば問題なさそうですが、地図上の駅からドロネー図を作成するなど、現実的に使用する場合にはかなり厳しそうな計算量・実行時間になることが予想できるので、できれば総当り・貪欲法は避けたいです。。。

点を一つずつドロネー図に追加しながら、逐次ドロネー分割を行っていく方法

今回は、乱択による逐次追加法と呼ばれる方法に近い方法で考えていきましょう。

既存のドロネー図に、点pを追加
1. 点群Pから、ランダムに未追加の点pを1つ取得。
三角形の分割
1. 既存のドロネー図の三角形群から、内部にpが存在する三角形△ABCを探し、△ABP, △BCP, △CAP に分割
2. 新たに作成した△ABP, △BCP, △CAP について、外接円にドロネー図上の他の点を含んでいないかをチェックする。
  1. ここで、含んでいた場合、何らかの処理をして、点を含まないように分割を修正。
  2. 点を含んでいない場合、そのままとする。
点群Pから新たな点pを取得。1に戻り、これを繰り返していく。
完成

大体こんな感じでしょうか。ドロネー図を管理するクラスにTriangle型とPVector型の何らかのコレクションと、必要となりそうなメソッドをもたせてやれば、実現できそうですね。

ただ、実際は以下のような初期設定と、それに伴う整形処理が必要となるので、これらの機能も持たせておきましょう。

初期設定
1. ドロネー三角形の元にする点群Pの生成
2. 点群をすべて囲う事のできる三角形(SuperTriangle)を生成。
整形処理
1. 初期設定で作成していたSuperTriangleの頂点3つと、これにつながる辺をすべて削除する。

HashMap,HashTableなどを使い、適切にデータ構造を構築したアルゴリズムの場合、$n$個の頂点をもつドロネー図の計算は理想的な場合は$O(n \log_2 n)$となります。
（辺をキーにして三角形を持つ、親に分割前の三角形を持つ、等）

※ユークリッド最小全域木とFortuneのアルゴリズム、乱択、分割統治法などそれぞれ実現方法はあるようです。こちらを参照。

将来的には最速での実装を目指したいとは思いますが今回は簡易化の為、ハッシュマップなどは用いずに三角形を総当りしていくような形の実装をしていきます。このため、計算量は$O(n^2)$とかになると思います。

(点数が増えるとどんどんと計算に時間がかかっていくこととなりますが、今回はドロネー三角形分割の概要についてを主に解説したいので、ご勘弁を...。そのうち別の解説記事を書きたいと思います。)

ここまで整理してみて、はじめに書いた三角形の外接円の内部に点があるかを判定するメソッドだけでは足りなさそうな事がわかりました。

とりあえず必要なメソッドは以下です。

内部にpが存在する三角形を探す。
外接円に点を含んでいるものを見つけた場合、点を含まないように処理する。
SuperTriangleの作成
SuperTriangleの削除

以降は、これらのメソッドの実装について解説していきます。

三角形の内部に点があるか判定する

真円の内部に点があるかどうかの判定は、半径が一様に同じなのでシンプルでしたが、三角形の場合は同じようにはできません。

ここでは、内部に点があるかどうかをベクトルの外積によって判定していきます。
ベクトルAとベクトルBの外積は、AをBに向けて回転させたときの右ねじの親指の方向になる。点pが三角形の外側になるときは、全部の外積の方向が揃わないので、判定できる。という内容を応用しています。

外積ってなんじゃらほい。外積でなんで内外がわかるんや！という方はこちら。

Triangleクラスサンプルコード

Triangle.java

public class Triangle
{
    /* 中略 */
    //外積(cross)を用いて、三角形の内部にpがあるかを調べる。
    //内側にある場合、外積のベクトルはすべて同じ向きになる。
    public boolean IsInsideOfTriangle(PVector p)
    {
        PVector A = this.A;
        PVector B = this.B;
        PVector C = this.C;

        //AB x BP, BC x CP, CA x AP
        PVector ABxBP = Cross(A, B, B, p);
        PVector BCxCP = Cross(B, C, C, p);
        PVector CAxAP = Cross(C, A, A, p);

        return (ABxBP.z >=0 && BCxCP.z >=0 && CAxAP.z>=0) || (ABxBP.z <=0 && BCxCP.z <=0 && CAxAP.z<=0);
    }
}

PVector外積サンプルコード

Cross.java

public PVector Cross(PVector A, PVector B, , PVector C, PVector D)
{
    //return ABxCD
    PVector AB = PVector.sub(A, B);
    PVector CD = PVector.sub(C, D);
    return AB.cross(CD);
}

ここでの注意点は、ベクトルの外積を計算するための順番です。
時計回りか反時計回りに関わらず、回転方向を統一しなければなりません。

また、ベクトルの向きにも注意です。PVectorクラスのsubメソッドを使用して、頂点間の線分を表すベクトルを取得していますが、これの引数の順番によっても、ベクトルの向きが変わるので、統一しなければなりません。

このあたりの真面目な解説はまた別の機会などで...。

先程の円の内部に点があるかの判定はシューティングゲームなどでも使われていますが、こちらもゲームなどでよく使われている判定方法です。格ゲーや、円の組み合わせだけで表現できない複雑な形の物体の当たり判定などには、この方法を用いて、三角形を並べることで物体の当たり判定としています。

外接円に点を含んでいるものを見つけた場合、点を含まないようにする(Flip処理)

先に、三角形の外接円に点pを含むかどうかを判定する機能を実装したので、次は、これで見つけた点pを含む三角形に、点pを含まないようにする何らかの処理が必要です。
ここでは、この処理の事を辺フリップ(Flip)と呼ぶことにします。

点pが追加されたときに、新たにできた三角形のうち一つを△ABP、点pを外接円に含んでしまっている三角形を三角形△DBAとする。
辺ABが不正な辺なので、辺AB(=辺DB)を削除し、新たに辺PD(=辺DP)を作成。
これらの辺からできる三角形2つを三角形△ADP、△DBPとして、保存する。
新しくできた三角形2つに対しても、点pを外接円に含む三角形が周囲にないか確認。
これを繰り返し、点pを外接円に含む三角形が存在しなくなるまでFlipを繰り返す。

Flipを実装するサンプルコードはこのようになります。

Flip処理　サンプルコード

Flip.java

public class DelaunayTriangulation
{
    /* 中略 *

    // 点pによる分割でできた新しい三角形を格納するスタック dividedから、一つを取得。
    Triangle ABP = divided.pop();

    // 三角形ABCにおいて、点pの向かいの辺ABを取得。
    // この辺が不正な辺かどうかを考えていく。
    Edge AB = GetOppositeEdge(ABP, p);

    // ドロネー図の三角形群の中から、辺ABを共有する三角形△ADBを探して、取得。
    Triangle ADB = GetTriangleShareEdge(ABP,AB,baseTriangles);
    PVector D = GetVertexPoint(ADB, AB);

    // △ABPの外接円に点Dが含まれる場合は、辺ABをFlipする。
    // Flip後、新たにできた三角形は、再度チェックする必要があるため、dividedのスタックに積む。
    if(Contains(ABP,D))
    {
        //Flip処理を実行
        Deque<Triangle> FlipedTriangles = Flip(ADB, AB, p);
        for(Triangle t : FlipedTriangles) divided.push(t);
    }
    else
    {
        newTriangles.push(ABP);
        newTriangles.push(ADB);
    }

    Deque<Triangle> Flip(Triangle ADB, Edge AB, PVector p)
    {
        Deque<Triangle> FlipedTriangles = new LinkedList<Triangle>();

        PVector D = GetVertexPoint(ADB,AB);
        PVector A = AB.start;
        PVector B = AB.end;

        FlipedTriangles.push(new Triangle(A,D,p));
        FlipedTriangles.push(new Triangle(D,B,p));

        return FlipedTriangles;
    }
}

Flip処理を実装しなくとも、分割後の三角形をハッシュで管理し、衝突するものは不正な三角形であるという性質を利用する方法もあるのですが、今回は解説から省きたいと思います。

superTriangleの生成

大きめに作成します。具体的には下記の通り。

キャンバスを四角形として考え、この四角形の外接円を作成。
四角形の外接円を内接円とする巨大な三角形をsuperTriangleとする。

実装にあたっては、こちらを参考にしました。

http://tercel-sakuragaoka.blogspot.com/2011/06/processingdelaunay_3958.html

superTriangleの生成

superTriangle.java

public Triangle GenerateSuperTriangle()
{
    PVector center = new PVector(width/2, height/2);
    float radius = PVector.dist(center, new PVector(width, height));

    return  new Triangle(
        new PVector(center.x - radius*sqrt(3), center.y - radius),
        new PVector(center.x + radius*sqrt(3), center.y - radius),
        new PVector(center.x, center.y + 2*radius)
    );
}

初期化の段階で、ドロネー図の三角形リストにはこのsuperTriangleだけが存在している状態とし、1点追加したらこれを3つに分割していくこととなります。

SuperTriangleの頂点3つにつながる辺をすべて削除

最後の処理です。
これがないと、左上の方、右上の方、下の方に線がのびていってしまうことになるのできちんと綺麗にしましょう。

Finalize処理サンプル

finalize.java

処理としては、すべての三角形のリストの中から、superTriangleのいずれかの頂点を頂点に持つ三角形を除外していくといったものになります。
なので、頂点が同じ座標かどうかを確認するIsEqualメソッドを実装して処理に使っています。

完成

というわけで、ここまでの処理をなんやかんやして組み立てて、必要なクラスを設計して、出来上がったものがこちらになります。

(クラスの設計に関しても解説は書きたかったけど今回はここまでで手一杯だよ。。。)

https://github.com/sashisusesouyu/DelaunayTriangleByProcessing

上記の説明では、点群を順次追加していくように解説していましたが、今回のプログラムではわかりやすいように

左クリック : 点を、マウスの座標に追加(AddPoint)
右クリック : superTriangleを除去(Finalize処理)

としています。

また、分割がわかりやすいように適当な色を付けたりしています。
なお、デバッグ用で、各三角形の外接円を表示できるようにしています。(コメントアウトで消しています。)

おわり

いかがでしたｋ。。。

こういったアルゴリズムを学ぶには、学校で学んだような数学などの内容や、ゲームでよく使われているような処理なども学び直すことができるので、良い機会になると思います。

学校などでは、アルゴリズム = ソートというような学習になる場合もあるかと思いますし、行列の計算など、

今後の課題

同一座標に点が存在するとき、計算の処理が正確にできない。今回は、全部の点が違う場所に位置している前提で書いたが、実際にはバリデーター等を書いて、安全に処理する必要がありそう。
ボロノイ図をドロネー図から作成する。
各クラスの設計、などなど
- メソッドをどこに持たせるか
- ファクトリメソッドやコントローラーなどを書いて、引数の型をみて処理を変える、
- メソッドやクラスの命名方法
計算量が現状、O(n^2)なので、点の量が増えると高速に処理ができなくなってくる。各三角形に対して、分割の過程を残しておくことで、計算量をO(n log2,n)に抑える事ができる。

三角形分割と同時に木構造を作って再帰的に求めていく。

参考文献