まずは、前知識として 『詳解 OpenCV』 の記述をふたたび引用しよう。
これをいかに解釈し、ソースコードに落とし込んでいくか。
- 外部三角形を作り、その頂点の 1 つを開始点とする（これにより、必ず外側の点から開始されることになる）。
- 内部の点を追加する。その後、すべての三角形の外接円を探索し、追加した点を含むような三角分割を削除する。
- 今削除した三角分割の外接円の内部にある、新しい点も含めて、グラフを三角分割し直す。
- 追加する点がなくなるまで、ステップ 2 に戻って繰り返す。
計算幾何に興味がある方はもちろん、普段ネット上のソースをコピペして『動きさえすればそれでよい』と思っている方にも、この記事がプログラミングの楽しみを知るきっかけになれば幸いである（偉そうなこと言ってごめんなさいごめんなさい）。
【必要なデータ構造】
以下に示すように、Delaunay 分割では、点の集合を三角形の集合に変えるため、まず前提としてそれらを表すクラスが必要になる。
↓
それらのために、ぼくは Point クラスと Triangle クラスを用意した。
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);
}
}
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) );
}
}
いずれも、さして特別な性質を持たない。
しいて言えば、Triangle クラスに「他の三角形と共有点を持つか」の判定を加えた点がチャームポイント。
先にタネを明かすと、hasCommonPoints メソッドは、処理の最後で外部三角形を除去する際の判定に使うものである。
さて。
点も三角形も、何らかのコレクションに格納して管理される事になる。
少なくとも Delaunay 分割を行っている間は、与えられた点集合には何の変更も与えない（つまり要素数が増減する事もない）が、三角形の集合は処理が完了するまで絶えず変動する事になる。
そこで、静的な点集合は、配列と同じ感覚で扱える List （ArrayList）で、一方動的な三角形集合は、追加・削除・検索（参照）が高速な HashSet で管理する事にし、それぞれ pointList、triangleSet と名付けた。
これで下ごしらえは完了である。いよいよ処理の内容に踏み込んでいく。
【逐次添加法】
冒頭のアルゴリズムは、『もともと分割されている Delaunay 三角形領域の内部に、新しい点を一つずつ加えていき、その都度分割をやり直す』というものであり、この事から逐次添加法とも呼ばれる。
例を示そう。
以下のように、既に領域が分割されているところに、黄色い点を加えるとする。
このとき、既に分割が完了している各三角形の外接円を視覚化すると、以下のようになる。
これから追加しようとする点は、上図より明らかに複数の三角形の外接円に包含されているため、外接円特性を満たすように、分割し直す必要がある。
以下の図は、分割し直した結果である。ご覧の通り、どの三角形の外接円も、新たに追加した点を包含する事はなくなっている。
これを、追加する点がなくなるまで繰り返せば、Delaunay 分割が完了する。
【外部三角形の作り方】
逐次添加法で Delaunay 分割を行う場合、はじめに与えられた点集合をすべて包含する架空の『外部三角形』を作らなければならない（そうしないと、『既に Delaunay 分割されている領域の内部に点を追加する』事ができないからである）。
 |
| 外部三角形は大きすぎて見えないので、各頂点から外部三角形の各頂点に線を引いてみた。 |
今回であれば、この問題を『スクリーン（表示領域）全体を包含する正三角形を求める』ものと見なして考えても差し支えなかろう（ちなみにこの方法は、一般化可能な理論である）。
スクリーンは以下のような矩形をなしている。そこでまずは、スクリーンとなる矩形領域を求める事から始める。
ゆえに、中心座標 c は (width / 2, height / 2) である。
次に、その矩形を包含する円を求める。
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 メソッド
// ======================================
// 最初に必要な巨大三角形を求める
// ======================================
// 画面全体を包含する正三角形を求める
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);
}
【三角形の外接円】
さて、新たに点を追加する際、既に分割されているどの Delaunay 三角形のセットに対しても、外接円特性が満たされなければならない。
つまり、もしも新たな点に対して、それを内包するような外接円をもつ Delaunay 三角形が存在する場合、それをセットから 削除する必要があるのだ（削除後、後述する再分割を施す）。
ここでは、新たに追加する点 p が、任意の三角形の外接円の内側に存在するかどうかの判定を行う議論を行う。
まず円を表現する Circle クラスを作る。これは、中心座標と半径を保持するものである。
Circle クラス
class Circle {
// 中心座標と半径
Point center;
float radius;
// ======================================
// コンストラクタ
// 中心座標と半径を与えて円をつくるよ
// ======================================
public Circle(PVector c, float r){
this.center = c instanceof Point ? (Point)c : new Point(c);
this.radius = r;
}
// ======================================
// 描画（デバッグ用）
// ======================================
public void draw() {
ellipse(center.x, center.y, radius * 2, radius * 2);
}
}
外接円の中心座標 c と半径 r が判れば、点 p と c の距離 |p - c| が r より大きいか小さいかを判定する事で、p が外接円の内部にあるかどうかを調べる事ができる。
もし r より小さい場合、その外接円を持つ三角形はセットから除外する必要がある事を意味する。
DelaunayTriangles.DelaunayTriangulation メソッド（一部）
// --------------------------------------
// 追加された点が外接円内部に存在する場合、
// その外接円を持つ三角形をリストから除外し、
// 新たに分割し直す
// --------------------------------------
if (Point.dist(c.center, p) <= c.radius) {
次に、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
この方程式を x と y について解くと、
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 メソッド
// ======================================
// 三角形を与えてその外接円を求める
// ======================================
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);
}
【再分割】
新たな点の追加に対して外接円特性を満たさない三角形は、三角形セットから削除し、再分割を施す。
削除に関しては、先日紹介したイテレータを使えばよい。
では、再分割はどうすればよいだろう。これがけっこう難しい。
ひとつひとつ考えてみよう。
追加しようとする点を内包するような外接円がひとつだけの場合は簡単である。
『追加する点と、それまでの三角形の各頂点のうち任意の 2 点を結んでできるすべての三角形』を、新たな分割結果として採用すればよい。
たとえば、新たに追加しようとしている点 p が、点 <p1, p2, p3> から成る三角形 T の外接円の内部に入っている場合、セットから T を消去し、代わりに <p, p1, p2>, <p, p2, p3>, <p, p2, p3> をセットに追加する。
もちろん、これで Delaunay 分割の条件である外接円特性は満たされている。
では、以下のように、新たな点を内包する外接円が複数ある場合はどうだろうか。
この場合は、先ほどのように簡単にはいかない。
『追加した点と、それを内包する三角形をの任意の 2 頂点からできる三角形』を新たな領域として見なすと、以下のように、明らかに不適当なものも含まれる。
実際には、赤色で示した三角形は除外すべきであり、これならばどの点も Delaunay 三角形の外接円の内部に入る事はない。
どうやって不適切な三角形を追加候補から除外すればよいだろう。
実は、三角形のセットに対して、前述の通り反復的に外接円特性を判定し、新たな追加候補の三角形を（不適当である可能性のものも含めて）新たな Delaunay 領域としてセットに追加していった場合、Delaunay 領域として不適切な三角形のみが重複する。
従って、新たに分割し直した追加候補の三角形は、とりあえず一時コレクションに格納しておき、重複判定の後に元の三角形セットへ追加すればよい。
……とはいうものの、必ずしもこれが最もスマートな方法とは言えず、ぼく自身、悩んだ末の苦肉の策である。
重複管理ができるデータ構造としては、HashMap を用いるのが定石だ。管理対象となるオブジェクト（ここでは三角形）を key に、重複しているかどうかの判定を value にする。
今回は、value を Boolean のオブジェクトにし、重複していなければ TRUE、重複していたら FALSE を持たせるようにした。
この重複管理用のメソッドは、以下のように書ける。
DelaunayTriangles.addElementToRedundanciesMap メソッド
// ======================================
// 一時ハッシュを使って重複判定
// hashMap
// - Key : 三角形
// - Value : 重複していないかどうか
// - 重複していない : true
// - 重複している : false
// ======================================
private void addElementToRedundanciesMap
(HashMap<Triangle, Boolean> hashMap, Triangle t) {
if (hashMap.containsKey(t)) {
// 重複あり : Keyに対応する値にFalseをセット
hashMap.put(t, new Boolean(false));
} else {
// 重複なし : 新規追加し、
hashMap.put(t, new Boolean(true));
}
}
このメソッドは、以下のように使われている。
DelaunayTriangles.DelaunayTriangulation メソッド（一部）
// --------------------------------------
// 追加された点が外接円内部に存在する場合、
// その外接円を持つ三角形をリストから除外し、
// 新たに分割し直す
// --------------------------------------
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();
}
重複判定ののち、三角形セットに正当な領域を追加する処理は、以下のように書ける。
DelaunayTriangles.DelaunayTriangulation メソッド（一部）
// --------------------------------------
// 一時ハッシュのうち、重複のないものを三角形リストに追加
// --------------------------------------
for(Iterator tmpIter = tmpTriangleSet.entrySet().iterator();
tmpIter.hasNext();) {
Map.Entry entry = (Map.Entry)tmpIter.next();
Triangle t = (Triangle)entry.getKey();
boolean isUnique =
((Boolean)entry.getValue()).booleanValue();
if(isUnique) {
triangleSet.add(t);
}
}
【後始末】
この一連の処理を、すべての点に対して反復的に行った後、一番はじめに追加した外部三角形のいずれかの頂点を持つ Delaunay 三角形を削除すれば完了だ。
DelaunayTriangles.DelaunayTriangulation メソッド（一部）
// 最後に、外部三角形の頂点を削除
for(Iterator tIter = triangleSet.iterator(); tIter.hasNext();){
// 三角形リストから三角形を取り出して
Triangle t = (Triangle)tIter.next();
// もし外部三角形の頂点を含む三角形があったら、それを削除
if(hugeTriangle.hasCommonPoints(t)) {
tIter.remove();
}
}
以上が、Delaunay 分割のアルゴリズムと実装方針である。
ソースコードをすべて見たい人は昨日の記事を読んでね！
というか、ここまで読んでくれてありがとうありがとう！ おつかれさま (ﾊﾞﾝﾊﾞﾝ
初めまして、一応フロントエンドを担当していますが、数学が苦手でこの手の解説をずっと探してました。。わかりやすい解説ありがとうございます。
返信削除外部三角形の求め方で躓いていたのですが次に進めそうです。ありがとうございました。
返信削除貴重な記事をありがとうございます！
返信削除非常に分かりやすいです。
やっと大体コードの書き方等まで理解できた〜
返信削除