断り書き
ドロネー図やボロノイ図を扱った論文や記事は多いですが、今回のこの記事は下の論文
に沿って基本的に書いています。また、記事中の図の多くは上記論文を元にした(論文の図からそのまま引用していいような)内容を大量に使わせて頂いています。
上記論文は３０年以上前の論文ということで現在では内容が最新でないかもしれないですが、ボロノイ図やドロネー図の論文としては参考として取り上げられていることが多いようです。
ただし申し訳ないですが、自分は主に実装目的で読み進めたので、全文理解できているとは言い難いですし、また完璧に再現した実装でもないです。
読める方は上記論文を直接読んだ方が理解が正確だと思います。
また、下の投稿記事も参考にしています。
内容が非常に分かりやすいです。特にquad-edgeの実装はこの投稿記事を参考にしています。
今回実装したコードは以下にあります。終わりの方は力尽きてたのでリファクタリングのし甲斐があります
以下のページにサンプルもおいてます。よければ遊んで行ってみてください。
ボロノイ図とドロネー図とは何か
wikiの方が詳しいのでささっと
ボロノイ図
与えられた母点に対して最寄りとなる領域で分割された図をボロノイ図といいます。
https://ja.wikipedia.org/wiki/%E3%83%9C%E3%83%AD%E3%83%8E%E3%82%A4%E5%9B%B3
ドロネー図
与えられた母点に対して、それを三角形で分割した際に、その任意の三角形の外接円が他の母点を含まないもの、となります。
これはテクスチャを作る際に使用されたりします。
https://ja.wikipedia.org/wiki/%E3%83%89%E3%83%AD%E3%83%8D%E3%83%BC%E5%9B%B3
ボロノイ図とドロネー図の関係
このドロネー図とボロノイ図は双対の関係にあり、ドロネー図にとっての面がボロノイ図の点であり、ボロノイ図の面がドロネー図の点になります。
ボロノイ図の作図に関してはFortuneのアルゴリズムや母店を逐次追加する方法がありますが、今回はquad-edgeと呼ばれるデータ構造を用いることで、ドロネー図の作図をボロノイ図を書くための中間ステップとします。
なので、タイトルにはドロネー＝ボロノイ図としていますが、内容としてはドロネー図が中心となります。ボロノイ図の話はドロネー図からの構成方法について語る最後の段階で出てきます。
データ構造(quad-edge)
今回はドロネー=ボロノイ図のデータ構造としてquad-edgeを採用しています。
quad-edgeというのは一個の辺に対して起点、目的点。及び辺を回転してできる辺への参照がある構造です。[1]
辺に対しては起点、目的地、右、左があり、さらに他の辺への遷移があります。
論文には$Onext$や$Lnext$等の複数の辺の遷移を辺に対する演算として記述されています。
図にするとこんなかんじです。[2]
いろいろ種類があります。名称のquadは４つの参照先があるという意味ではあるのですが、辺のネットワーク上の遷移を考える上では参照先として絶対に必要なのは$Onext$と$rot$だけで、他はその組み合わせで全部表現できます。
例えば後々の実装上で必要なものを取り上げると、
\begin{align}
sym &= rot\ rot \\
tor &= rot\ rot \ rot \\
Lnext &= tor\ Onext\ rot\\
Oprev &= rot\ Onext\ rot \\
Rprev &= sym\ Onext \\
\end{align}
という形です。
$Onext$は$next$と略し、この$next$と$rot$の演算を元に、他の演算をプログラム上で構成していきます。
$next$という処理は時計回りに針を回す、一方で$rot$という処理は$90$度回りながら面の世界、すなわちボロノイ図の世界に行くという感覚です。定義や実装から見るようにquad-edgeは参照系です。座標に関してはあまり意識しすぎない方がいいかもしれません。大事なのは位置関係とネットワークです。このquad-edge上の演算のネットワークを使えば、ドロネー＝ボロノイ系を辿っていきやすい、というのがquad-edgeの主張です。
今回は以下の実装を行いました。実装に関しては断り書きで入れたブログの記事を大いに参考にさせて貰っています。
pub type QuadEdge<T> = Arc<RwLock<QuadEdgeBase<T>>>;
type QuadEdgeOption<T> = Option<QuadEdge<T>>;
#[derive(Clone, Debug)]
pub struct QuadEdgeBase<T> {
origin: Option<PointRef<T>>,
rot: QuadEdgeOption<T>,
next: QuadEdgeOption<T>,
}
impl<T> QuadEdgeBase<T> {
pub fn _rot(&self) -> QuadEdgeOption<T> {
match &self.rot {
Some(quad_edge_rc) => Some(Arc::clone(quad_edge_rc)),
None => None,
}
}
pub fn _next(&self) -> QuadEdgeOption<T> {
match &self.next {
Some(quad_edge_rc) => Some(Arc::clone(quad_edge_rc)),
None => None,
}
}
fn rot(&self) -> QuadEdge<T> {
self.read().unwrap()._rot().unwrap()
}
fn next(&self) -> QuadEdge<T> {
self.read().unwrap()._next().unwrap()
}
fn sym(&self) -> QuadEdge<T> {
self.rot().rot()
}
fn tor(&self) -> QuadEdge<T> {
self.sym().rot()
}
fn lnext(&self) -> QuadEdge<T> {
self.tor().next().rot()
}
fn prev(&self) -> QuadEdge<T> {
self.rot().next().rot()
}
fn rprev(&self) -> QuadEdge<T> {
self.sym().next()
}
}
RwLockとArcを使用しているのは半分は読み込み管理のため、半分は描画に使用するbevy側で要求されるからです。$next$や$rot$の操作は重ねて頻発するため、このまま使うとunwrap地獄になるので、そうならないように_rot()や_next()をかませています。
張ったコードは少し最終的な実装と違いますが、今回のプログラムで使う基本操作は入っています。
さて、quad-edge上では、ネットワーク作成のため辺の接続、辺の消去の操作を行いたいです。これらを実装するため、次の２つの基幹的な操作を定義します。
makeEdge
辺の作成です。実装と一緒にquad-edgeのイメージ画像を貼ります。
pub fn make_quad_edge<T>(origin: PointRef<T>, dest: PointRef<T>) -> QuadEdge<T> {
let start_end_ref = QuadEdgeBase::ref_new();
let end_start_ref = QuadEdgeBase::ref_new();
let left_right_ref = QuadEdgeBase::ref_new();
let right_left_ref = QuadEdgeBase::ref_new();
let mut start_end_mut = start_end_ref.write().unwrap();
let mut end_start_mut = end_start_ref.write().unwrap();
let mut left_right_mut = left_right_ref.write().unwrap();
let mut right_left_mut = right_left_ref.write().unwrap();
start_end_mut.origin = Some(origin);
end_start_mut.origin = Some(dest);
start_end_mut.rot = Some(right_left_ref.clone());
end_start_mut.rot = Some(left_right_ref.clone());
left_right_mut.rot = Some(start_end_ref.clone());
right_left_mut.rot = Some(end_start_ref.clone());
start_end_mut.next = Some(start_end_ref.clone());
end_start_mut.next = Some(end_start_ref.clone());
left_right_mut.next = Some(right_left_ref.clone());
right_left_mut.next = Some(left_right_ref.clone());
return start_end_ref.clone();
}
実装にある通り、右左と左右の辺は同一接続元ですが、方向は区別されます。
また、作成した時点では他のquad-edgeとの参照関係もできていないことに注意してください。つまり辺$ab$と辺$ac$のquad-edgeを作ったとしても、この関数だけでは$ab$と$ac$の参照関係ができてないので、$ab$の$next$が$ac$に対応しません。
参照関係の作成は次の$splice$という処理を使用します。
splice
$splice$という処理は直感的ではないです。基幹処理と捉えてなんとなく理解する方が良いかもしれません。
辺$a$及び$b$に対し$splice[a,b]$は
(1) 辺$a,b$が異なる起点の場合は同一起点にする
(2) 辺$a,b$が同じ起点の場合はそれを分割する(1の逆)
という動作をします。
(3) もし2つが同じリングで、向きが逆なら、Spliceはそのリングの順序を逆転させる(直訳)
という内容の３番目の効果もあるのですが、特殊なケースなので基本的に今回実装する処理には1と2があればなんとかなったことと、それ以上に自分の方で理解が十分におよばなかったので括弧書きで付け加えておきます。詳しくは原文の論文96ページにあるので、もし把握できた方がいたらコメント頂けると助かります。
具体的な$splice$の内容はquad-edgeの$next$と$next\ rot$の参照先を入れ替える処理です。実装を見た方が理解しやすいと思います。
pub fn splice<T>(a: QuadEdge<T>, b: QuadEdge<T>) {
swap_nexts(a.next().rot(), b.next().rot());
swap_nexts(a.clone(), b.clone());
}
fn swap_nexts<T>(a: QuadEdge<T>, b: QuadEdge<T>) {
let anext = a.clone().next();
let bnext = b.clone().next();
let mut a = a.write().unwrap();
//NOTE: This checks whether a and b are same points.
// If a = b, b.try_write() will return Err because it can't get write lock.
match b.try_write() {
Ok(mut b_w) => {
a.next = Some(bnext);
b_w.next = Some(anext);
}
Err(_) => {
return;
}
}
}
$splice$という同じ操作内容で切断と接合ができるのは不思議ですが、内容は入れ替えですので、入れ替えを２回繰り返せば確かに元に戻ります。
$makeEdge$と$splice$を合わせて図に書くと下図のようになります。
これを使用してさらに高次元の操作を定義します。
connect
２つのquad-edgeの間にさらにquad-edgeを追加する処理です。
実装は以下です。
pub fn connect<T>(a: QuadEdge<T>, b: QuadEdge<T>) -> QuadEdge<T>
where
Point<T>: VectorCalc,
QuadEdge<T>: QuadEdgeCalc<T>,
{
let a_dest = a.get_dest();
let b_orig = b.get_orig();
let ab = make_quad_edge(a_dest, b_orig);
splice(ab.clone(), a.clone().lnext().clone());
splice(ab.sym(), b);
ab
}
deleteEdge
quad-edgeの削除です。実装は以下です。
pub fn delete_edge<T>(e: QuadEdge<T>) {
splice(e.clone(), e.prev());
splice(e.sym(), e.sym().prev());
}
ちなみに$rot$中のボロノイ図の描画で必要となる座標は必要に応じて作成されます。ボロノイ図のグラフはドロネー図から逆再生できるので、ドロネー図の構造が作成された後で十分だからです。具体的にはドロネー図の各三角形の外心を求めればいいです。ボロノイ図の座標の生成に関しては今回の記事の最後の方に改めて記述します。
分割統治法によるドロネー図の構造決定
ドロネー図の書き方は点を逐次追加する方法と領域分割する方法、Fortuneのアルゴリズム等があります。
このうち、分割統治の計算量のオーダー的には期待的には$O(nlogn)$で逐次追加法と変わりないです。点を一個ずつ追加する方法(逐次添加法)の場合は、一度生成された後からさらに点が増える場合に強いのですが、一方で事前に点が所与であるという前提であれば、最悪計算量が$O(nlogn)$になる分割統治の方が安定しています。
ちなみに今回紹介する分割統治方法を改良した平均$O(nloglog(n))$の方法[3]があるようですが今回は踏み込まずにベーシックな方法になります。
分割統治のアルゴリズム
atcoderみたいなプロコンに参加している方は名前からピンとくるかもしれないですが、与えられた母点の集合$S$を$x$軸を基準に左右に分けて、下から順番に統合することの繰り返しを考えます。
そのために、母点はx軸で昇順にソートされていることを前提にします。x軸の値が同一ならばy軸の昇順です。
このときに左の母点の集合を$L$、右を$R$として、$L-R$の橋渡しをする基底となる辺を$basel$とします。この$basel$から上に向かって、ドロネー三角分割のもう一点の候補を$L$と$R$の端の辺の中から一つずつ求めることを繰り返します。$basel$はドロネー分割の三角形の一辺なので、そこから繋ぐべき母点は$L$か$R$に属す一点です。なので、その母点に向かって$basel$の両端からさらに２辺を引いて一つのドロネー分割された三角形が形成されます。
では繋ぐべき母点をどうやって見つければよいでしょうか？イメージは$basel$から上に向かって円を膨らませていく形になります。膨らませた円が最初に触れた母点が次に接続するべき点になります。
この風船を膨らます動作に対応するアルゴリズムを考えましょう。
まずLとRの合計が2点の間はその間のquad-edgeを引き、3点の場合は素直にquad-edgeの三角形を作ればいいです。なので合計が4点以上の場合を考えます。このとき$L$や$R$中の三角の分割の候補については、$basel$の起点や終点から、辺を順繰りに接合側から順番に辿っていけば見つかります。
この候補となる辺を$lcond$や$rcond$と書きます。左なら$lcond$。右ならば$rcond$です。$lcond$の起点は$basel$の起点、$rcond$の起点は$basel$の目的地です。
ドロネー三角分割の定義から、もし正しい分割がされているならば、候補地点($lcond$や$rcond$の目的地)と$basel$の起点と終点の外接円には、他の母点を含みません。実はこの時に、$basel$の起点と目的地、および現在の$lcond$($rcond$)の目的地の３点で作られる外接円がその次の$lcond$($rcond$)の目的地に含まれるかどうかの判定[4]を繰り返していき、最終的に次の$lcond(rcond)$の目的地が現在の$lcond(rcond)$と$basel$の外接円に含まれない点で$lcond$($rcond$)が決定できます。[5]
この作業を$L$と$R$に行うと、条件を満たす候補がそれぞれに最大１点ずつ求められます。最後は$L$と$R$の点のうち一番$basel$に近い点を次の三角分割の対象として設定します。
また、候補地点に到達するまでに辿った辺($lcond$と$rcand$)については、結合時にドロネー三角の辺として不適になるので削除します。[6]
そして、基底となる$basel$を新しく設定して同じ操作を繰り返します。
もし求めた$lcond$も$rcond$も$basel$より下に来て不適ならば[7]、$basel$が最上部に来たと判定してループが終わり、$L$と$R$の統合を完了します。
アルゴリズムの実装
実装に関しては論文中の41ページに疑似コードがあり、それをrustにして後々必要な処理を追加して実装しています。疑似コードについては論文中にあるので、自分の実装のうち分割統治に関係する所を抜き出したものを投下します。
渡されている点は事前に$x$座標から$y$座標の順番にソートされていることに注意してください。
pub fn divide_conquer(mut points: Vec<PointRef<Vec2>>) -> (QuadEdge<Vec2>, QuadEdge<Vec2>)
{
match points.len() {
1 => {
panic!("points should have more than 2 elements.")
}
//点２つなので直線
2 => {
let (p1, p2) = (points[0].clone(), points[1].clone());
let a = make_quad_edge(p1.clone(), p2.clone());
let a_sym = a.sym();
(a, a_sym)
}
//点が３つしかないので三角形を作成
3 => {
let (p1, p2, p3) = (points[0].clone(), points[1].clone(), points[2].clone());
let a = make_quad_edge(p1.clone(), p2.clone());
let b = make_quad_edge(p2.clone(), p3.clone());
splice(a.sym(), b.clone());
if ccw(p1.as_ref(), p2.as_ref(), p3.as_ref()) {
let _ = connect(b.clone(), a.clone());
(a, b.sym())
} else if ccw(p1.as_ref(), p3.as_ref(), p2.as_ref()) {
let c = connect(b.clone(), a.clone());
(c.sym(), c)
} else {
(a, b.sym())
}
}
_ => {
let latter = points.split_off(points.len() / 2);
let (ldo, mut ldi) = divide_conquer(points);
let (mut rdi, rdo) = divide_conquer(latter);
//初期のldiとrdiの候補を取得[8]
loop {
if ldi.is_left_point(rdi.get_orig().clone()) {
let ldi_next = ldi.lnext();
ldi = ldi_next;
continue;
} else if rdi.is_right_point(ldi.get_orig().clone()) {
let rdi_prev = rdi.rprev();
rdi = rdi_prev;
continue;
}
break;
}
let mut basel = connect(rdi.clone().sym(), ldi.clone());
//結合後の次のldiとrdiの候補を調整
let ldo = if ldi.is_same_origin(ldo.clone()) {
basel.sym()
} else {
ldo
};
let rdo = if rdi.is_same_origin(rdo.clone()) {
basel.clone()
} else {
rdo
};
//統合処理
loop {
//L側のdelaunay分割の候補(lcand)を特定,通過したquad-edgeは削除
let mut lcand = basel.sym().next();
if basel.is_valid(lcand.clone()) {
while in_circle(
(
basel.get_dest().as_ref(),
basel.get_orig().as_ref(),
lcand.get_dest().as_ref(),
),
lcand.next().get_dest().as_ref(),
) {
let t = lcand.next();
delete_edge(lcand.clone());
lcand = t;
}
}
//R側のdelaunay分割の候補(rcand)を特定,通過したquad-edgeは削除
let mut rcand = basel.prev();
if basel.is_valid(rcand.clone()) {
while in_circle(
(
basel.get_dest().as_ref(),
basel.get_orig().as_ref(),
rcand.get_dest().as_ref(),
),
rcand.prev().get_dest().as_ref(),
) {
let t = rcand.prev();
delete_edge(rcand.clone());
rcand = t;
}
}
//baselが最上部なら動作停止
if !basel.is_valid(lcand.clone()) && !basel.is_valid(rcand.clone()) {
break;
}
//lcand,rcandのうちより近い外接円に含まれる方を選択し,次のbaselを設定
basel = if !basel.is_valid(lcand.clone())
|| in_circle(
(
lcand.get_dest().as_ref(),
lcand.get_orig().as_ref(),
rcand.get_orig().as_ref(),
),
rcand.get_dest().as_ref(),
) {
let resp = connect(rcand.clone(), basel.sym());
resp
} else {
let resp = connect(basel.sym(), lcand.sym());
resp
};
}
(ldo, rdo)
}
}
}
これを元に書きくだしてプログラムを動かして、後述する座標の取り出しのプログラムを得て生成された図をみてみます。
できました。
ちなみに分割統治にかかる時間の大半は統合の時間です。なのでいかに分割方法を工夫するかが高速化の鍵です。
今回はベーシックに$x$座標に応じて左右に分け、$y$軸の下から行う実装にしていますが、他の分割方法に関しては例えばピザみたいに上下の次は左右と分割していく方法があります。
ボロノイ図の逆算
最後に、ドロネー図からボロノイ図を逆算します。
ボロノイ図用の点の追加
ボロノイ図用の点に関しては点の座標が用意されていないので、埋める処理を行います。基本的にはドロネー三角分割の外心を入れていけばいいです。
注意しなくてはいけないのは、$splice$の処理でボロノイ図の点は接合されてないという点です(下図の左)。$p$をボロノイ図上のquad-edgeの起点とすると、$p.next$の起点$q$はまだ$q\neq p$です。なのでこれを接合する処理をするか、もしくは接合しなくても代入時に一致させる処理を入れるかです。ここで代入漏れがあると、続く処理でNoneが出力されてエラーになります。今回は分割統治の三角形の作成時に接合する処理を追加で入れました。
ただし、端の場合には注意が必要です。この場合は三角形の外心がなく、これは接合してはいけません。$q\neq p$で正しいのです。接合すると下図の右のようになってしまいます。
端点に関しては分割統治の完了後にしかわからないので、探査を行い、接合せずに描写のために遠方の点をquad-edgeの端点とし、この点にのみ値を代入します。
ボロノイ図のネットワークの抽出
こちらも分割統治から得られたquad-edgeの一辺から探査します。
またボロノイ図とドロネー図の区別はrotの使用回数の偶奇により決まります。なので、nextとsymのみを使用してネットワークを抽出すれば、入力したquad-edgeがボロノイ図対応かドロネー図対応かで、自動的にそれに対応したネットワークを抽出できます。
例えば、
let daunay = extract_network(l);
let voronoi = extract_network(l.rot());
という形です。ただしextract_network(l:QuadEdge)をネットワークを抽出する関数としています。
ちなみにボロノイ図に点を追加する際にもl.rot()を入口にしています。
結果
描画ができました。完成したプログラムのコード以下にあります。
wasm化したものが自分のおもちゃ箱においてます。よければ遊んで行ってください。
性能
計算に必要な時間を計算しました。
下の出力用調整は、分割統治＋ドロネー図の座標抽出＋ボロノイ図の逆算の座標出力までの時間です。
| 母点数 | 分割統治の完了時間 | 出力用調整完了時間 |
|---|---|---|
| 1000 | 0.149s | 0.178s |
| 2000 | 0.317s | 0.375s |
| 3000 | 0.488s | 0.572s |
| 5000 | 0.871s | 1.013s |
| 10000 | 1.899s | 2.191s |
なので時間が母点数にほぼ比例するのを見ると計算量はきちんと$O(NlogN)$に収まっていそうですね。
おまけ
母点が常に動くケース
母点が動く場合は、辺の変動を計算するよりは全部を再度生成しても問題ないのではないかと考えているので、より良いアルゴリズムがあるかもしれないですが、それでも各母点毎の計算が必要なことを考えると$O(n)$は割らないはずなので、今回は都度都度の領域分割に頼ることにします。
一方で描画処理の際の辺のオブジェクトに関しては再利用しない場合は、辺の削除と生成を毎フレームごとに繰り返すことになり、愚直に更新を使わずに実行すると、今度はテーブルへの変更回数自体がオーバーヘッドになりそうです。
辺の数が一定だったらそのまま上書きすればいいので話は早いのですが、母点の位置関係により変化します。これは点が４つの例でも明らかです。
今回の使用上では見栄え重視で点があまり多い事例で実施するつもりもなく、十分高速だったので、妥協して、基本は値の更新で、過不足ある場合はオブジェクトの生成や削除をすることにしました。
辺数の変化は適切に捕らえられていることと、辺の対応は維持されていないものの値の変更自体は問題ないはずなので、テーブル削除・追加の回数でのボトルネックの危惧は回避されているはずです。
bavyというかECSの自体の使用方法にも依存しますが、よりよい方法があるかもしれません。
終わり
というわけで記事としてはこんな感じです。
読んでくださってありがとうございます。
補足
[1]
ここに[2]で後述するの$Flip$の裏表の情報を含めて本当は８パターンありえますが、ドロネー＝ボロノイ図の作図という目的において、この$Flip$の情報は必須ではないので本文中では省きます。
また、Ｃで以下の実装が公開されています。
４つの参照があるのが分かると思います。この実装は論文共作者のStolfi教授監修らしく、本来想定されているQuadEdgeの実装です。
[2]
また論文中では他に$Flip$という動作があります。ここでの$Flip$は進行方向に向かって180度上下回転するという動作です。進行方向に対して表向きか裏向きかの概念と言った方がいいかもしれないです。$Flip$は概念上はあると便利です。$Flip$があると$rot$や$next$の方向を反転させることができます。
論文中には、他にも双対性を示す$Dual$という演算があります。双対性を示す$Dual$という演算は２回で元の辺に戻る必要があります。この$Dual$は、ドロネー図からボロノイ図、ボロノイ図からドロネー図の対応でもあります。つまり、点と面の入れ替えです。後述するquad-edge的には上下と左右ですね。
点から面への移動に対応する演算は$rot$です。しかし$rot$だけだと元の辺に戻るには$rot$が4回必要になります。一方で$Dual$は2回で戻る必要があります。なので$Dual≠rot$です。操作を工夫する必要があります。もちろん、$rot\ rot$を$Dual$とするわけにはいきません。これは$sym$であり、上下間の反転なわけですから。
ここで$Flip$があれば、$Flip\ rot$の演算が$Dual$に対応できます。$rot$が$90$度回転だと仮定すれば、$Flip\ rot$を2回行うということは、裏の状態で$90$度回転、つまり$-90$度回転し、また裏返ると表の状態になり$90$度回転することです。結果、表向きで$0$度回転した状態で、最初の状態になります。
$Flip\ rot$が$Dual$に対応することにより、$Dual$が定義できる。そのために$Flip$が必要なのだと自分は$Flip$の必要性を解釈しています。他にも論文中では$Flip$は$next$の逆操作を作る際に使われています。
が、実は他の記事や論文中では、$Flip$をその意味合いではあまり使いません。大抵はドロネー分割中に向かい合う三角形の辺を正しく入れ替える動作を$Flip$と言っています。
ちなみに余談の余談ですが、同論文中では他の記事では$Flip$になっている辺の入れ替えは$swap$と呼ばれています。
[3]
ちなみに$O(NlogNlogN)$の方法は1996年時点ですと大体のケースで最良みたいです。
[4]
二次元平面上の点$(a,b,c)$の外接円に対して点dが含まれるかの判定は次の公式を利用します。
{\displaystyle {\begin{aligned}&{\begin{vmatrix}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{vmatrix}}>0\end{aligned}}}
証明
1.三次元上の点$(a,b,c,d)$に対して四面体の体積の公式は以下の公式で計算できる
{\displaystyle {\begin{aligned}&\frac{1}{6}\ {\begin{vmatrix}A_{x}&A_{y}&A_{z}&1\\B_{x}&B_{y}&B_{z}&1\\C_{x}&C_{y}&C_{z}&1\\D_{x}&D_{y}&D_{z}&1\end{vmatrix}}\end{aligned}}}
2.円上の点と円内の点を放物面に移した四面体の体積の考察から外接円の判定条件を導出する
の２つのパートに分けます。
1.四面体の体積の公式の証明
前提として、三次元上の点$(a,b,c)$及び原点を中心にした四面体の体積の公式は以下になります。
{\displaystyle {\begin{aligned}&\frac{1}{6}\ {\begin{vmatrix}A_{x}&A_{y}&A_{z}\\B_{x}&B_{y}&B_{z}\\C_{x}&C_{y}&C_{z}\end{vmatrix}}\end{aligned}}}
この公式に帰着します。
行列の基本変形から元の式の四行目の値を各行から引く変形をしても行列式の値は変わりません。
よって、
{\displaystyle {\begin{aligned}&\frac{1}{6}\ {\begin{vmatrix}A_{x}-D_{x}&A_{y}-D_{y}&A_{z}-D_{z}&0\\B_{x}-D_{x}&B_{y}-D_{y}&B_{z}-D_{z}&0\\C_{x}-D_{x}&C_{y}-D_{y}&C_{z}-D_{z}&0\\0&0&0&1\end{vmatrix}}\end{aligned}}}\\
={\displaystyle {\begin{aligned}&\frac{1}{6}\ {\begin{vmatrix}A_{x}-D_{x}&A_{y}-D_{y}&A_{z}-D_{z}\\B_{x}-D_{x}&B_{y}-D_{y}&B_{z}-D_{z}\\C_{x}-D_{x}&C_{y}-D_{y}&C_{z}-D_{z}\\\end{vmatrix}}\end{aligned}}}
これは三次元上の点$(a-d,b-d,c-d)$及び原点の四面体の体積の公式に一致します。
余談ですがこの行列式の判定は０以上かどうかの判定になり結構シビアです。この判定を失敗するとアルゴリズムが成立しません。使用していたライブラリが四次元以上の場合はLU分解してから行列式を求めていたせいか、今回用いるためには、精度がかなり悪かったので、今回の実装では三次元に直し、行が同一の場合は０を早期リターンして使ってます。
さて、1はこれだけなのですが、四面体の体積の公式のこのバリアントは検索してもなかなか出てこないのと、次の証明の鍵なので書いておきます。
2.円上の点と円内の点を放物面に移した四面体の体積の考察から外接円の判定条件を導出する
円の公式
(x-p)^2+(x-q)^2=r^2
を展開して整理すると
-2px -2qy+(x^2+y^2)+(p^2+q^2-r^2)=0
になります。これを変形すると
\displaylines{
-2px -2qy+(x^2+y^2)+(p^2+q^2-r^2)=0 \\
x^2+y^2=2px+2qy-(p^2+q^2-r^2)
}
です。この上の式は左辺は放物面の公式、右辺は平面の公式に対応します。つまり、二次元平面での同一円周上の点を$(x,y,x^2+y^2)$で表される放物面上に移した場合は、必ずその放物面上の点も同一平面上にあるということです。
立体絡むとさすがに汚い図だと通らないのでPrimitives for the Manipulation of General Subdivisions and the Computation of Voronoi Diagrams p108より引用
ここで例えば、４点$(a,b,c,d)$が同一円周上にある場合、それを放物面上に移した４点$(a',b',c',d')$で構成される四面体の体積は$0$になります。４点は同一平面上にあるからです。
一方で、もし４点のうち$d$が円内にある場合は、平面と放物面の位置関係より$d'$は平面部分よりも下に位置します。
つまり４点$(a',b',c',d')$で構成される四面体は今度は値を持ちます。
このとき四面体の面積の正負を考えると、$a,b,c$が反時計回りという前提より正になります。
よって1で導出した公式を用いて四面体の体積の正負を表現すると、
{\displaystyle {\begin{aligned}&{\begin{vmatrix}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{vmatrix}}>0\end{aligned}}}
になることが示せました。
[5]
便宜上,左側($lcond$)のケースのみ考えます。右の場合は左の場合を反対にすればよいです。
まず$basel$を延長した直線を$l$とします。この時に、ドロネー分割の辺を反時計回りに辿ります。その$N$番目の$lcond$の目的地の候補地点を$Ni$とします。このとき、ある点$A$より右にある任意の点$X$に対し、$AXN_{i}$の$l$以上の部分を$Γ_{i}$とすると、$i$がある$k$以下では$Γ_{i}⊇Γ_{i+1}$となり、それ以上では$Γ_{i}⊆Γ_{i+1}$になることが示せます。
つまり、端から順番に辿っていけばやがて次の候補地点を含まない外接円が見つかって、それが求める点になる、ということが言えます。
証明は以下のようにできます。
三角形$AN_{i}N_{i+1}$の外接円と$basel$を延長した直線の交点を$X_{i}$とします。$X_{i+1}$は、$X_{i}$よりも左側にあり、$AN_{i}N_{i+1}$の外接円に含まれています。$X_{i}$は$i$が増えるにつれ単調に点$A$に近づいていきます。
$X$が$X_{i}$より左にある限り、$X$は$A_{x}N_{i}N_{i+1}$の外接円に含まれます。言い換えると、$AXN_{i}$の外接円に$N_{i+1}$が含まれているということです。なので$Γ_{i}⊇Γ_{i+1}$です。
もし$X_{i}$が$X$よりも左になった場合(つまり$i$がある$k$以上)、逆に考えればいいので、今度は$Γ_{i}⊆Γ_{i+1}$になります。
上で証明した内容をAをbaselの起点(目的地)、Xを目的地(起点)とすれば良いです。$Γ_{i}$の抱合関係はドロネー分割の目的地として不適であることを示しているので、抱合関係が成立しなくなる点が次のドロネー分割の対象です。
[6]
これはそれまでに辿った辺($lcond$や$rcond$)が候補地点と$basel$を結んだ直線と交差することや、$basel$と$lcond$の目的地で構成される三角形の外接円が次の$lcond$の候補地点を含むことを考えると良いです。
[7]
また、辺に対する左右の判定は、次の公式が成立するならば二次元平面上の３点$(a,b,c)$半時計回りであることを利用します。
{\displaystyle {\begin{aligned}&{\begin{vmatrix}A_{x}&A_{y}&1\\B_{x}&B_{y}&1\\C_{x}&C_{y}&1\end{vmatrix}}>0\end{aligned}}}
中身は[4]の1番目の公式の二次元版です。証明もほぼ同様ですので省略です。
これを利用して３点が半時計回りかどうかを判定し、それを通して辺と点の関係の左右を区別することで、求められた$lcand$の目的地がドロネー三角分割の対象として有効か無効かの判定を行います。例えば$a$に$lcand$の起点、$b$に$lcand$の目的地、$c$に判定の対象となる点を適応し、もし式が成立するならば、$(a,b,c)$は反時計回りなので、判定対象の点は$lcand$より左に位置することが分かります。
[8]
ここで取得したい辺は結合の開始となる辺、つまり向かい合わせになったドロネー図を結ぶ最下部の辺です。
言い換えると、すべての左右のドロネー図の点が同一方向にあるような辺です。
まず、右側の点を固定したケースを考えます。
そして左側のドロネー図上の点が直線に対して同一方向にあるようなドロネー図外周上の点を求めます。
そして、左側のドロネー図上で最も右にある点から直下外周上の点に向かう直線(つまりプログラム中でループされている$rdo$)から、ドロネー図外周を時計周りに辺を回していきます。
このとき、もし外周上の辺の向きが右側の点に対して逆になれば、向きが変わった２辺の共通点と右側の点との直線が、条件を満たすような左側の点です。ドロネー図は凸だからです。
右側でも同様のことを考えます。右側に対しての起点は$ldo$で半時計回りです。
一回だけの繰り返しだと下まで降りてないかもしれないので、左側と右側両方に対して変化がなくなるまでループすれば、求める辺が導出できます。
divide_conquer関数の帰り値中の$ldo$と$rdo$はこの部分ために保持されています。