ProcessingでDelaunay分割（実装篇）

前々回のあらすじ： 砕け散るエフェクトを作った。

……それはよいのだが、いかんせんやっつけ仕事だったため、砕け散る破片の形状生成がびっくりするほど適当だった。

もうちょっとうまく破片を作るために、今回は計算幾何のアルゴリズムの中でも比較的有名な Delaunay 分割に挑戦してみた。

まずはじめに、結果から示そう。

Delaunay 分割は、ランダムに与えられた点を結び、下のような無数の三角形を作る手法である。

さて。

Delaunay 三角分割法に関しては、Gary Bradski, Adrian Kaehler 著、松田 晃一 訳『詳解 OpenCV』にこんな解説がある。

Delaunay 三角分割法は、空間内の点を連結して三角形のグループにし、その三角形のすべての角に対する最小角度が最大になるようにするテクニックで、1934年に発明されました。　（中略）　与えられた任意の三角形の頂点に接する任意の円の内部には、その他の頂点が含まれないようにして三角分割が行われています。これは、外接円特性と呼ばれています。

正直、さっぱり意味がわからない。Delaunay 分割を全くご存知ない読者がこれを一度読んだだけで完全理解できたら大した読解力だと思う。少なくともぼくは数日間悩んだ。

---

ところで、この本の p.304 には、Delaunay 分割の『比較的単純なアルゴリズム』が概説されている。当該箇所を引用しよう。

今日、Delaunay 三角分割を計算するアルゴリズムはたくさんあります。そのいくつかはとても効率的ですが難解です。比較的単純なアルゴリズムの 1 つを要約すると以下のようになります。

1. 外部三角形を作り、その頂点の 1 つを開始点とする（これにより、必ず外側の点から開始されることになる）。
   
   
2. 内部の点を追加する。その後、すべての三角形の外接円を探索し、追加した点を含むような三角分割を削除する。
   
   
3. 今削除した三角分割の外接円の内部にある、新しい点も含めて、グラフを三角分割し直す。
   
   
4. 追加する点がなくなるまで、ステップ 2 に戻って繰り返す。

このアルゴリズムの複雑さのオーダーは、データ点の数に対して O(n²) です。（後略）

……。

なにこの不親切な解説。

特に 3. の『グラフを三角分割し直す』って、そもそもどうやって分割し直せばよいのだろう。

そして、アルゴリズムの解説書は数多く出版されているが、こういう計算幾何の領域を解説した本となると途端に少なくなるから困ったものである（あっても上述のように解説が果てしなく不親切だったりする）。

というわけで、同じように悩んでいる人がいるかも知れないなぁと思い、ぼくなりに作ったソースコードをノーカット・フルバージョンで公開してみる事にした。

なお、今回はいくつものクラスを作る事になるので、クラスごとにタブを作ってコードを分けている。

作ったのは、点を表す Point クラス、三角形を表す Triangle クラス、円を表す Circle クラス、そして Delaunay 三角分割のための DelaunayTriangles クラスである。

特に、Point クラスは、Processing の PVector を継承しており、equals メソッドによる同値判定ができるようになっている。

Point オブジェクトは Triangle に集約され、さらにそれらは Delaunay 分割の各領域を表すオブジェクトとして用いられる。

Circle は、外接円特性の判定の便宜上、用意したものである。

---

【使い方】

いきなりだが、DelaunayTriangles クラスの使い方から見ていこう。

DelaunayTriangles delaunay;

void setup() {
  size(400, 300);
  noLoop();
  
  // ランダムに頂点のリストを作成
  ArrayList ptList = new ArrayList();
  for(int i = 0; i < 500; i++) {
    ptList.add(new Point(random(width),random(height)));
  }
  
  // それを基に、DelaunayTrianglesオブジェクトを生成
  delaunay = new DelaunayTriangles(ptList);
}

void draw() {
  background(255);
  stroke(0);
  noFill();
  smooth(); 

  delaunay.draw();
}

たったこれだけで、Delaunay 分割ができてしまうのである。

ちなみに、ほんのちょっとパフォーマンスは落ちるが、10行目の

ptList.add(new Point(random(width),random(height)));

でリストに追加している Point オブジェクトを、おなじみの PVector に書き換えた

ptList.add(new PVector(random(width),random(height)));

であっても現実的な速度で動く。

以下、あんまり綺麗ではないが、実装例として各クラスのソースコードをご紹介する。

---

【たのしいソースコード】

Pointクラス

class Point extends PVector {
  // ======================================
  // コンストラクタ
  // ======================================
  public Point() {
    super();
  }
  public Point(float x, float y) {
    super(x, y);
  }
  public Point(float x, float y, float z) {
    super(x, y, z);
  }
  public Point(PVector v) {
    this.x = v.x;
    this.y = v.y;
    this.z = v.z;
  }
  
  // ======================================
  // 同値判定
  // ======================================
  public boolean equals(Object o) {
    boolean retVal;
    try {
      PVector p = (PVector)o;
      return (x == p.x && y == p.y && z == p.z);
    } catch (Exception ex) {
      return false;
    }
  }

  // ======================================
  // 描画
  // ======================================  
  public void draw() {
    point(x, y);
  }
}

---

Circle クラス

class Circle {
  // 中心座標と半径
  Point center;
  float radius;
  
  // ======================================
  // コンストラクタ
  // 中心座標と半径を与えて円をつくるよ
  // ======================================
  public Circle(PVector c, float r){
    this.center = new Point(c);
    this.radius = r;
  }
  
  // ======================================
  // 描画（デバッグ用）
  // ======================================  
  public void draw() {
    ellipse(center.x, center.y, radius * 2, radius * 2);
  }
}

---

Triangle クラス

class Triangle{
  public Point p1, p2, p3;  // 頂点
  
  // ======================================
  // コンストラクタ
  // 3頂点を与えて三角形をつくるよ
  // 頂点はPointで与えてもOK
  // ======================================
  public Triangle(PVector p1, PVector p2, PVector p3){
    this.p1 = p1 instanceof Point ? (Point)p1 : new Point(p1);
    this.p2 = p2 instanceof Point ? (Point)p2 : new Point(p2);
    this.p3 = p3 instanceof Point ? (Point)p3 : new Point(p3);
  }
  
  // ======================================
  // 同値判定
  // ======================================
  public boolean equals(Object obj) {
    try {
      Triangle t = (Triangle)obj;
      // ※ 同値判定に頂点を用いると、
      // 三角形の頂点の順番を網羅的に考慮する分条件判定が多くなる。
      return(p1.equals(t.p1) && p2.equals(t.p2) && p3.equals(t.p3) ||
             p1.equals(t.p2) && p2.equals(t.p3) && p3.equals(t.p1) ||
             p1.equals(t.p3) && p2.equals(t.p1) && p3.equals(t.p2) ||
              
             p1.equals(t.p3) && p2.equals(t.p2) && p3.equals(t.p1) ||
             p1.equals(t.p2) && p2.equals(t.p1) && p3.equals(t.p3) ||
             p1.equals(t.p1) && p2.equals(t.p3) && p3.equals(t.p2) );
    } catch (Exception ex) {
      return false;
    }
  }
    
  // ======================================
  // ハッシュ表で管理できるよう、hashCodeをオーバーライド
  // ======================================
  public int hashCode() {
    return 0;
  }
  
  // ======================================
  // 描画
  // ======================================  
  public void draw() {
    triangle(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y);
  }
  
  // ======================================
  // 他の三角形と共有点を持つか
  // ======================================  
  public boolean hasCommonPoints(Triangle t) {
    return (p1.equals(t.p1) || p1.equals(t.p2) || p1.equals(t.p3) ||
            p2.equals(t.p1) || p2.equals(t.p2) || p2.equals(t.p3) ||
            p3.equals(t.p1) || p3.equals(t.p2) || p3.equals(t.p3) );
  }
}

---

DelaunayTriangles クラス

class DelaunayTriangles {
  public HashSet triangleSet;  // 三角形リスト

  // ======================================
  // コンストラクタ
  // 与えられた点のリストを基にDelaunay分割を行う
  // ======================================
  public DelaunayTriangles(ArrayList pointList) {
    DelaunayTriangulation(pointList);
  }
  
  // ======================================
  // 点のリストを与えて、Delaunay三角分割を行う
  // ======================================
  public void DelaunayTriangulation(ArrayList pointList) {
    // 三角形リストを初期化
    triangleSet = new HashSet();
    
    // 巨大な外部三角形をリストに追加
    Triangle hugeTriangle = getHugeTriangle();
    triangleSet.add(hugeTriangle);

    try {
      // --------------------------------------
      // 点を逐次添加し、反復的に三角分割を行う
      // --------------------------------------
      for(Iterator pIter = pointList.iterator(); pIter.hasNext(); ) {
        Object element = pIter.next();
        Point p = element instanceof Point ? 
            (Point)element : new Point((PVector)element);
        
        // --------------------------------------
        // 追加候補の三角形を保持する一時ハッシュ
        // --------------------------------------
        // 追加候補の三角形のうち、「重複のないものだけ」を
        // 三角形リストに新規追加する
        //          → 重複管理のためのデータ構造
        // tmpTriangleSet
        //  - Key   : 三角形
        //  - Value : 重複していないかどうか
        //            - 重複していない : true
        //            - 重複している   : false
        // --------------------------------------
        HashMap tmpTriangleSet = new HashMap();

        // --------------------------------------
        // 現在の三角形リストから要素を一つずつ取り出して、
        // 与えられた点が各々の三角形の外接円の中に含まれるかどうか判定
        // --------------------------------------
        for(Iterator tIter=triangleSet.iterator(); tIter.hasNext();){
          // 三角形リストから三角形を取り出して…
          Triangle t = (Triangle)tIter.next();
              
          // その外接円を求める。
          Circle c = getCircumscribedCirclesOfTriangle(t);
              
          // --------------------------------------
          // 追加された点が外接円内部に存在する場合、
          // その外接円を持つ三角形をリストから除外し、
          // 新たに分割し直す
          // --------------------------------------
          if (Point.dist(c.center, p) <= c.radius) {
            // 新しい三角形を作り、一時ハッシュに入れる
            addElementToRedundanciesMap(tmpTriangleSet,
              new Triangle(p, t.p1, t.p2));
            addElementToRedundanciesMap(tmpTriangleSet,
              new Triangle(p, t.p2, t.p3));
            addElementToRedundanciesMap(tmpTriangleSet,
              new Triangle(p, t.p3, t.p1));
            
            // 旧い三角形をリストから削除
            tIter.remove();            
          }
        }
        
        // --------------------------------------
        // 一時ハッシュのうち、重複のないものを三角形リストに追加 
        // --------------------------------------
        for(Iterator tmpIter = tmpTriangleSet.entrySet().iterator();
            tmpIter.hasNext();) {

          Map.Entry entry = (Map.Entry)tmpIter.next();
          Object t = entry.getKey();
          
          boolean isUnique = 
              ((Boolean)entry.getValue()).booleanValue();

          if(isUnique) {
            triangleSet.add(t);
          }
        }
      }
      
      // 最後に、外部三角形の頂点を削除
      for(Iterator tIter = triangleSet.iterator(); tIter.hasNext();){
        // 三角形リストから三角形を取り出して
        Triangle t = (Triangle)tIter.next();
        // もし外部三角形の頂点を含む三角形があったら、それを削除
        if(hugeTriangle.hasCommonPoints(t)) {
          tIter.remove();
        }
      }
    } catch (Exception ex) {
      return;
    }
  }

  // ======================================
  // 描画
  // ======================================
  public void draw() {
    for(Iterator it = triangleSet.iterator(); it.hasNext(); ) {
      Triangle t = (Triangle)it.next();
      t.draw();
    }
  }
  
  
  // ======================================
  // ※ ここからprivateメソッド
  // ======================================
  
  
  // ======================================
  // 一時ハッシュを使って重複判定
  // hashMap
  //  - Key   : 三角形
  //  - Value : 重複していないかどうか
  //            - 重複していない : true
  //            - 重複している   : false
  // ======================================
  private void addElementToRedundanciesMap(HashMap hashMap, Object t)
  {
    if (hashMap.containsKey((Triangle)t)) {
      // 重複あり : Keyに対応する値にFalseをセット
      hashMap.put(t, new Boolean(false));
    } else {
      // 重複なし : 新規追加し、
      hashMap.put(t, new Boolean(true));
    }
  }
  
  // ======================================
  // 最初に必要な巨大三角形を求める
  // ======================================
  // 画面全体を包含する正三角形を求める
  private Triangle getHugeTriangle() {
    return getHugeTriangle(new PVector(0, 0), 
                           new PVector(width, height));    
  }
  // 任意の矩形を包含する正三角形を求める
  // 引数には矩形の左上座標および右下座標を与える
  private Triangle getHugeTriangle(PVector start, PVector end) {
    // start: 矩形の左上座標、
    // end  : 矩形の右下座標…になるように
    if(end.x < start.x) {
      float tmp = start.x;
      start.x = end.x;
      end.x = tmp;
    }
    if(end.y < start.y) {
      float tmp = start.y;
      start.y = end.y;
      end.y = tmp;
    }
    
    // 1) 与えられた矩形を包含する円を求める
    //      円の中心 c = 矩形の中心
    //      円の半径 r = |p - c| + ρ
    //    ただし、pは与えられた矩形の任意の頂点
    //    ρは任意の正数
    Point center = new Point( (end.x - start.x) / 2.0,
                              (end.y - start.y) / 2.0 );
    float radius = Point.dist(center, start) + 1.0;
    
    // 2) その円に外接する正三角形を求める
    //    重心は、円の中心に等しい
    //    一辺の長さは 2√3･r
    float x1 = center.x - sqrt(3) * radius;
    float y1 = center.y - radius;
    Point p1 = new Point(x1, y1);
    
    float x2 = center.x + sqrt(3) * radius;
    float y2 = center.y - radius;
    Point p2 = new Point(x2, y2);
    
    float x3 = center.x;
    float y3 = center.y + 2 * radius;
    Point p3 = new Point(x3, y3);

    return new Triangle(p1, p2, p3);    
  }
  
  // ======================================
  // 三角形を与えてその外接円を求める
  // ======================================
  private Circle getCircumscribedCirclesOfTriangle(Triangle t) {
    // 三角形の各頂点座標を (x1, y1), (x2, y2), (x3, y3) とし、
    // その外接円の中心座標を (x, y) とすると、
    //     (x - x1) * (x - x1) + (y - y1) * (y - y1)
    //   = (x - x2) * (x - x2) + (y - y2) * (y - y2)
    //   = (x - x3) * (x - x3) + (y - y3) * (y - y3)
    // より、以下の式が成り立つ
    //
    // x = { (y3 - y1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)
    //     + (y1 - y2) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)} / c
    //
    // y = { (x1 - x3) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1)
    //     + (x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1)} / c
    //
    // ただし、
    //   c = 2 * {(x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)}
    
    float x1 = t.p1.x;
    float y1 = t.p1.y;
    float x2 = t.p2.x;
    float y2 = t.p2.y;
    float x3 = t.p3.x;
    float y3 = t.p3.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;
    Point center = new Point(x, y);
    
    // 外接円の半径 r は、半径から三角形の任意の頂点までの距離に等しい
    float r = Point.dist(center, t.p1);
    
    return new Circle(center, r);
  }
}

---

Processing と計算幾何ってけっこう親和性の高い分野だと思うのだけど、日本国内ではなかなか浸透しないなぁ、ぼくも興味あるんだけどなぁ……と少々残念だったりする。

『Processingで学ぶ計算幾何』とか出たら、きっとぼくだったらものすごく欲しがると思うのだけど……。

ちなみに、この実装方針とか、アルゴリズムの解釈とかは、長くなるので今日は書かない。

ただ、ソースコードには一部わかりにくすぎる箇所も存在するので、気が向いたらそこらへんの事も書き留めておきたいと思う（見出しに『実践篇』と書いたのは、その予告的な意味も含まれている……かも）。

---

【おまけ: 3次元 Delaunay 分割】

ちなみに、応用編として、3次元 Delaunay 分割というものも存在する。

こちらは、3次元の点群に対して Delaunay 分割を施すというものである。

例として、以下のような3次元点が与えられたとする。

これを、Delaunay 分割すると、以下のようになる。

さらに、得られた三角形に対して面を張ると、こんな感じになる。

これの実際に動くサンプルは、ぼくのサイトに載せているので、もしよかったらご覧ください。まる。

---

【おまけ（その2）: 簡易版実装】

実は、『どの頂点も、Delaunay 分割されたどの三角形の外接円の内部にも含まれない』という条件を満たすだけならば、もっと（人間にとって）簡単なアルゴリズムが存在する（ソースコードは展開してご覧ください）。

int N = 300;     // 点の数
PVector[] pVec;  // 頂点配列

void setup() {
  size(400, 300);
  pVec = new PVector[N];
  for (int i = 0; i < N; i++) {
    pVec[i] = new PVector(random(width), random(height));
  }  
  smooth();
  noLoop();
}

void draw() {
  background(255);
  
  stroke(0);
  for(int i = 0; i < N-2; i++) {
    PVector v1 = pVec[i];
    for(int j = i+1; j < N-1; j++) {
      PVector v2 = pVec[j];
      for(int k = j+1; k < N; k++) {
        PVector v3 = pVec[k];
        
        float tmp = 2.0*((v2.x-v1.x)*(v3.y-v1.y)-(v2.y-v1.y)*(v3.x-v1.x));
        PVector center = new PVector(
          ((v3.y-v1.y)*(v2.x*v2.x-v1.x*v1.x+v2.y*v2.y-v1.y*v1.y)+
           (v1.y-v2.y)*(v3.x*v3.x-v1.x*v1.x+v3.y*v3.y-v1.y*v1.y))/tmp,
          ((v1.x-v3.x)*(v2.x*v2.x-v1.x*v1.x+v2.y*v2.y-v1.y*v1.y)+
           (v2.x-v1.x)*(v3.x*v3.x-v1.x*v1.x+v3.y*v3.y-v1.y*v1.y))/tmp           
        );        
        float r = PVector.dist(center, v1) - 0.01;
        
        boolean flg = false;
        for (int l = 0; l < N; l++) {
          if (PVector.dist(center, pVec[l]) < r) {
            flg = true;
            break;
          }
        }
        if (!flg) {
          line(v1.x, v1.y, v2.x, v2.y);
          line(v2.x, v2.y, v3.x, v3.y);
          line(v3.x, v3.y, v1.x, v1.y);
        }
      }
    }
  }
}

このソースコードは単独で期待した振る舞いを見せるのだが、計算効率が恐ろしく悪いため（まさかのO(n⁴)）、決しておすすめはできない。

---

【8月3日追記】

『まつにおまかせ』の I☆MAT さんが、ActionScript バージョンを作成されたようです。

記事中で当ブログを参照なさっています。ありがとうございました！