ボロノイ図を描く！ ~Fortune's Algorithm~

こんにちは！ Ysmr-Ry です．
Twitter の知り合いにボロノイ図の論文を見せてもらって興味が湧いたのでボロノイ図の描き方を勉強してみました！
ボロノイ図とはなんなのか，そしてその描き方を見てみましょう．
「ビーチに波が押し寄せる」ようにボロノイ図ができてゆきます．
長くなってしまいましたが，なぜ放物線なのかがわかれば十分な気がします．

ボロノイ (Voronoi) 図とは

ボロノイ図とは，与えられたいくつかの母点のうちどれに一番近いかで領域を分けた図のことです．
実際にボロノイ図を見てみましょう．

これがボロノイ図です．黒い点が母点 (ボロノイ図を描くために与えられた点) で色分けされてるのがそれぞれのボロノイ領域です．
よく見てると気づくとは思いますが，隣接するボロノイ領域の境界線は対応する母点を結んだ線分の垂直二等分線です．

球面上で作ったり，いろいろ拡張はあるようです．
次のリンク先は地球上の空港を母点として球面上にボロノイ図を描いています．

World Airports Voronoi

Fortune's Algorithm

ボロノイ図の描き方は当然いろいろあるとは思いますが，今回紹介するのは Fortune's Algorithm という方法です．
Fortune's Algorithm は平面走査法の一種で，円錐曲線としての放物線に着想を得ています．
簡単言うと，ボロノイ図を描きたい平面上を横線が上から下に動いていってスキャンし，それに放物線を重ね合わせた "波" が続いて，波の継ぎ目がボロノイ領域の境界をなぞっていきます．
実際見た方が早いですね．こんな感じです．(点の配置によってはバグるのでそういう時はうまくいくまでページを更新してください)

assqr.github.io

水色の線は放物線の組み合わせでできている "波" で，それよりも上の領域が "海" と考えると良いです．
黒い横線が平面上を走査して行って，緑色のボロノイ図が描かれていくのがわかります．
赤い丸はなんなのか，右の木 (こういう形のデータ構造を木と言います) は何を表しているのか，はこの後説明します．

前提知識 1 (円錐曲線としての放物線)

何と言ってもこのアルゴリズムの面白いところは放物線が生きてくるということです．
放物線，2 次関数自体は中学生でもやると思いますが，ここで使う放物線の性質は現在の課程では数 III の 2 次曲線のところで出てきます．
それは，「ある定直線 とその上にない定点 があるとき， と から等距離にある点の軌跡は放物線である」です．
定直線 を準線 (directrix)，定点 を焦点 (focus) と言います．
アニメーションで見るとこんな感じです．

紫の定点と定直線がそれぞれ焦点と準線，赤いのがそれによってできる放物線，オレンジ色の円は焦点を通り定直線に常に接しているような円です．
ということはこの円の中心は焦点と準線から同じ距離 (半径) にあるので，放物線の上を動くというわけです．
これをより一般化して焦点と準線からの距離の比が となるような点の軌跡を離心率 の二次曲線といい，放物線・楕円・双曲線になります．

前提知識 2 (Binary Search Tree)

先ほどの Voronoi Visualizer の右側に表示されていた木を二分探索木 (Binary Search Tree; BST) といいます．
木というのは先ほど見ていただいた通りの形をしたデータ構造なのですが，中でも二分木はそれぞれの頂点 (丸のこと) から伸びる枝が高々二股以下に分かれているものを言います．
それぞれの頂点には例えば数字が入っています．二分探索木は，「ある頂点に注目した時にその中に入っている値に対して，そこから左に生えている枝の先の部分木に入っている値は全てそれよりも小さく，右は全て大きいという性質を保つもの」です．
図で見ていきましょう．次のような木は二分探索木と言えます．

一番上の 5 の入った頂点 (根と言います) に注目します．すると，左の部分木は下の四角で囲った部分です．

見てみるとこの中のどの値も 5 より小さくなっていることがわかります．(1, 2, 3 < 5)
また右の部分技についてはどの値も 5 より大きくなっていることがわかります．(5 < 6, 8, 11)
今回は根に注目しましたがこれの性質はどの頂点についても成り立っています．(例えば 2 の頂点なら 1 < 2 で 2 < 3 です)
こういう風に数値のデータを持っていて何がうれしいかというと，探索木ですから，この中からある値を探す時に効率がいいのです．
例えばこの二分探索木から 6 を探したいとしましょう．
根から順にたどっていきます．根に入っている値は 5 です．5 < 6 なので，先ほど見た性質から，6 は少なくとも右の部分木にあることがわかるので，右の枝を辿ります．
これを繰り返せば，最悪の場合でも木の深さ分だけ比較をすれば見つかることがわかります．

木の深さ に対してその深さにある頂点はだいたい 個あるので計算量は (平衡状態であれば) となります．
Fortune's Algorithm では二分探索木を放物線の波の最前線を管理するのに使います．

Fortune's Algorithm (なぜ放物線を考えるのか)

では本題の Fortune's Algorithm について考えていきましょう．
Fortune's Algorithm は平面走査法なので，走査線 (Sweep Line) が平面上を上から下に動きます．線が掃いていくから Sweep です．
ただこのアルゴリズムの場合は特殊で，もう一つ放物線の組み合わせの "波" である汀線 (Beach Line) を考えます．
今回はボロノイ図を構成するためのアルゴリズムなので，走査線よりも上の掃き終わった領域においてはボロノイ図の一部が組みあがっていて欲しいのです．
しかしながら，走査線よりも上ならば下にある母点たちの影響を全く受けないとは言えません．
そこで，それよりも上ならば下の母点たちから一切影響を受けないというもう一つの境界線，汀線を考えるのです．
先ほど言ったように汀線の上側の領域を "海" と呼称すると，ボロノイ図の出来上がった領域である "海" がじわじわと平面全体を覆っていくということです．
ではなぜその汀線が放物線となるのか．それは前提知識 1 で見た放物線の性質によるものです．
汀線は走査線が掃き終わった，上にある全ての母点をそれぞれ焦点とし，現在の走査線を準線とした放物線よりも上の領域のうち，もっとも下にある境界線をつぎはぎしてできた最前線として決めます．
例えば，先ほどの Visualizer で放物線を全て表示すると

となるのですが，実際に考えるのは一番下にある，最前線の放物線の組み合わせです．
そうするとそれらの放物線は「走査線と対応する母点から等距離にある点の軌跡」なので，放物線の上の領域は「走査線よりも母点の方が近い点の集まり」になります．走査線よりも下の点は当然走査線よりも遠いですから，"海" は走査線よりも下の点から影響を受けないことがわかります．
また，隣接する放物線の交点を Breakpoint と呼ぶのですが，その点は準線と二つの母点からの距離が等しいのでボロノイ図の辺の上に乗ります．

Fortune's Algorithm (2 つのイベント)

平面走査法はアニメーションとしては連続的に走査線が動いていくのですが，実際のプログラムでは特別なイベントが起こる時点をとびとびにたどっていきます．
今回の場合そのイベントは 2 つあります．

Circle Event

図のように 3 つの母点が同一円周上にあり，その円の下側に走査線が接する時点を言います．(したがって上の図は Circle Event が起こる少し前)
このまま走査線が下に行くと，赤い円の中心に集まってきている 3 つの放物線の一部のうち，真ん中が完全に汀線からなくなります．
また，円の中心は Circle Event のまさにその時点では 3 つの放物線に同時に乗っていると言うことができ，それに対応する 3 つの母点 (と走査線) から等距離にあるため，ボロノイ図における点 (ボロノイ点) になります．
したがってこのイベントが起きた時は真ん中の放物線を汀線の放物線のリストの中から消去し，また，ボロノイ点で終わっている辺の終わりにボロノイ点をセットし，ボロノイ点から始まる辺の始点をセットします．

Site Event (NewPoint Event)

Site Event は母点の一つに走査線が触れた時点で起きます．上の図は走査線が 3 の母点を掃いた直後です．
母点 3 を焦点とする放物線が生じるため，それを汀線に組み込まねばなりません．

Fortune's Algorithm (実装の詳細)

先ほどのイベントは上から順番に並べて順に処理します．この時 y 座標が最小のものを効率よく取り出すデータ構造が必要で，Priority Queue を使うみたいです．最初は母点から Site Event だけを作って並べて置いて，それを処理ながら Circle Event を発生させてイベントリストに追加します．
このアルゴリズムで使うもう一つのデータ構造は先ほども紹介したように二分探索木です．
頂点の持つ値として何を入れているかというと，Breakpoint の x 座標です．
Visualizer で (1,2) と頂点に表示されていたら，母点 1 と母点 2 を焦点とする放物線が左からこの順で交わっている時の交点の x 座標が入っています．
Site Event のところの図がわかりやすいと思います．(再掲)

この場合は，(0,2) (2,3) (3,1) (1,0) の順で左から Breakpoint が並んでいるのがわかります．(このうち画面内に見えているのは (2,3) (3,1) だけですが…)
ちなみにこの順序は二分探索木を深さ優先探索 (Depth-First Search) すると得られます．

Circle Event においては真ん中の放物線を消し (つまり真ん中と左，真ん中と右の作る Breakpoint を消す)，新たな Breakpoint として円の中心を追加します．
Site Event では新たに放物線ができるので，既存の汀線とどこでぶつかるかを調べて，新たに Breakpoint を追加しなければなりません．
関係する汀線の部分を探索するのに効率が良いのです．

また実装の際には母点の座標と準線の位置さえあれば放物線も Breakpoint も計算できるのでそれさえ持っておけば OK です．
辺は始点だろうが終点だろうがどのみち Circle Event の時の円の中心になるのでその時決定すれば良いです．

詳しくは Haskell で実装してる人がいるのでそちらでどうぞ．

github.com

また，他の言語での実装もあります．(英語だけど)

Fortune’s algorithm and implementation

おわり

日本語の資料少ないので苦労しました…

http://www.cs.sfu.ca/~binay/813.2011/Fortune.pdf
http://cs.nyu.edu/yap/bks/egc/09/7Vor.pdf

ちなみに僕の Visualizer のソースコードは以下にあります．Haskell の Haste というライブラリを使って JavaScript のコードを生成しています．
したがって描画は HTML5 の Canvas です．
僕が確認しただけでも以下のようなバグがあるので，Pull-Request をいただけると助かります…

・最初の母点 0, 1 だけの放物線の場合，片方が表示されずかくかくする (治った？)
・点の配置による原因不明のフリーズ
・放物線の継ぎ目がはみだしてる (CircleEvent が発生していない) ← CircleEvent がおきその先は画面によって切り取られているときにおかしくなる

これアルゴリズムの性質上 Circle Event が起きないと辺が終わらないので画面に切り取られているような辺は自分で計算しなきゃいけないんですよね…
僕が改造したんですが，ボロノイ辺がありえないようなものになる場合があるみたいで…

だいぶ下調べはしたはずなんですけど，間が空いちゃって今は Haskell のソースちらちら見ながら書いてるだけなので，間違い等あったらご指摘お願いします…