- Fortune’s Algorithm はボロノイ図を O(n log n) 時間で作成できるが、実装が難しく、大規模な図を繰り返し生成するのでなければ O(n²) 実装やライブラリのほうが現実的
- ボロノイ図は複数の site を基準に平面を最も近い領域へ分割し、境界は 2 つの site と等距離にある点で構成される
- アルゴリズムは左から右へ動く sweep line と、放物線弧の前線である beachline を維持し、site event と circle event だけを処理する
- site event では新しい弧を挿入して既存の弧を分割し、incomplete edge を作る。circle event では中央の弧を取り除きながら Voronoi Vertex と half edge を完成させる
- 実用的な実装ではイベントキュー、beachline、incomplete edge のマップ、DCEL をまとめて扱う必要があり、無効な circle event の除去や残った edge の整理が複雑さを大きく増す
実装難度と適用範囲
- Fortune’s Algorithm はボロノイ図を O(n log n) 時間で生成するアルゴリズム
- 実際に使う目的なら、すぐに実装するより先に必要な規模を見極めるほうがよい
- 大きな図を毎秒何枚も生成する必要がないなら O(n²) 実装のほうが簡単な選択肢になりうる
- より現実的な代替案は既存の ライブラリ を使うこと
- アルゴリズム自体の動作結果は視覚的におもしろいが、実装過程は難しく、挫折しやすい
ボロノイ図の基本概念
- ボロノイ図は平面を複数の領域に分割する方法で、手続き型の地図生成でよく使われる
- 入力として選んだ点は site または seed と呼ばれる
- 各 site に対応する cell は、平面上の点のうちその site に最も近い点の集合
- セルの境界は 2 つの site と等距離にある点で構成される
- セルの角が交わる Voronoi Vertex は、3 つの site と等距離にある点
sweep line、beachline、event
- Fortune’s Algorithm は左から右へ移動する垂直線である sweep line を使う
- sweep line が site に達すると、その site を焦点とする放物線弧ができ、sweep line が離れるほど弧は大きくなる
- 異なる site に対応する 2 つの弧が交わる点は 2 つの site と等距離なので、セルの境界になる
- 2 つの境界が出会うと、図の頂点が生成される
- 活性な弧の前線は beachline と呼ばれる
- 実装では sweep line をピクセル単位で動かすのではなく、計算可能な特定の地点である event だけを処理する
- site event: あらかじめ分かっている site 座標で定義され、処理時に beachline に新しい弧が追加される
- circle event: beachline 上の 3 つの弧で定義され、処理時に 1 つの弧が取り除かれ、Voronoi Vertex と half edge が生成される
放物線で境界を見つける
- このアルゴリズムでは放物線を一般的な
y = ax^2 + bx + c の形ではなく、locus definition で扱う
- 放物線は 1 つの focus point と 1 本の directrix で定義される
- focus point は site になる
- directrix は sweep line になる
- 同じ sweep line を directrix に持つ 2 つの放物線の交点は、2 つの site と等距離にある
- したがって 2 つの放物線の交点を見つければ、2 つの site 間の equiedge を見つけられる
- 放物線の x 座標を計算する疑似コードと、sweep line の位置が変わると 2 つの放物線の交点が境界に沿って移動する例が使われる
beachline の表現と site event の処理
- beachline の各弧は、その site の座標だけで表現できる
- sweep line はすべての弧に共通して適用される
- 実装では弧を個別オブジェクトではなく 2D 座標 として扱う
- beachline は単純な点の並びとして表現できる
- 例:
[arc1, arc2], [arc1, arc2, arc3]
- 同じ site の弧が beachline に複数回現れることがある
- 例:
[arc1, arc3, arc1, arc2]
- site event が発生したら、新しい site から左に線を引いたときに最初にぶつかる beachline の弧を見つけ、その弧を新しい弧が 分割 する
- 新しい site
L が既存の beachline [.., i, j, k, ..] の j を分割すると、構造は [.., i, j, L, j, k, ..] になる
- site は x 座標順にキューへ入り、処理のたびに beachline とイベント候補が更新される
circle event と circumcircle
- beachline の 3 つの弧
[.., i, j, k, ..] で 2 本の境界が交わる状況が起こると、中央の弧 j が消える
- このとき 3 つの site を通る circumcircle が存在し、その中心は 3 つの site と等距離にある
- circumcircle の中心は Voronoi Vertex になる
- circle event は円の右端点である circle point を基準にイベントキューへ入る
- 新しい site が circle point に到達する前に円の内部で見つかると、既存の circle event は無効になる
- 新しい site が先に中央の弧を分割するため、3 つの弧の組み合わせがもう保たれない
- 元の triple
i, j, k は消え、i, j, L や L, j, k のような新しい triple を調べる必要がある
incomplete edge と half edge
- incomplete edge は、一方の端点は固定されているが、もう一方の端点が 2 つの放物線弧の交点で定義される線
- site event で新しい弧が挿入されると 2 本の incomplete edge が生成される
- 固定点は、新しい弧が既存の beachline と出会った座標
- 新しい弧
j が既存の弧 i を分割すると、[i, j], [j, i] の交点に対応する edge が生じる
- circle event で 2 本の incomplete edge が衝突すると、その衝突点が Voronoi Vertex になる
- 既存の incomplete edge はこの点で half edge として完成し、新たに隣接した 2 つの弧の間には新しい incomplete edge が生成される
反時計回りの円だけが circle event になる
- beachline に
[i, j, k, j, i] があるとき、ijk と kji はどちらも円を作れるが、両方が有効な circle event になるわけではない
- 中央の弧が消えるのは、境界が実際に収束する側だけ
- プログラムでは 3 点の orientation を determinant で判定する
- determinant が負なら反時計回りで、circle event になる
- determinant が正なら時計回りで、circle event にはならない
- determinant が 0 なら 3 点は一直線上にあり、円は存在しない
アルゴリズム全体の流れ
- 入力 site を x 座標でソートして site event としてキューへ入れる
- キューが空になるまで次の event を取り出して処理する
- site event の処理:
- これから先の circle event のうち、新しい site が円の内部に入るものを削除する
- 新しい site が分割する beachline の弧を見つける
- 新しい弧を挿入して既存の弧を分割する
- 2 本の incomplete edge を追加する
- 新しくできた triple が circle event を作れるか調べる
- circle event の処理:
- circumcircle の中心を Voronoi Vertex として追加する
- 中央の弧を beachline から取り除く
- 取り除かれた弧のために無効になる今後の circle event を削除する
- 新たに隣接した弧の triple を調べて circle event を追加する
- キューが空になったら、残っている incomplete edge を図の境界まで延長し、境界と交わる点に Voronoi Vertex を生成する
Odin 実装のデータ構造
- サンプル実装は C の代替言語である Odin で書かれている
- 全コードは RedPenguin101/voronoi リポジトリにある
- 基本型:
V2: [2]int 形式の 2D 点PointPair: 2 つの V2 の組Event: {site: bool, a, b, c: V2} 構造体
Event の意味は型によって変わる
- site event では
a が site 座標で、b, c は使わない
- circle event では
a, b, c が event を作った beachline 上の 3 つの弧
Fortune 構造体は次の状態を管理する
beachline: V2 配列queue: Event 配列incomplete_edges: PointPair -> V2 マップvd: ボロノイ図を保存する DCEL
実装で省略または単純化されている点
- beachline はベクタで表現しているが、効率を上げるなら binary tree のほうが適している
- event queue も概念上は priority queue だが、サンプル実装では配列にソート挿入する方式で扱っている
- circle event の無効化は将来のイベントを走査して調べる方式で、より高速な方法が必要だという TODO がある
clean_beachline_edges は beachline の両端で不要な弧を切り落とす手続き- 実装には、同じ x 座標を持つ site、circle point と site が同じ場合、参照点の衝突といった例外処理が含まれる
- キューが空になった後に残る incomplete edge、twin を持たない half edge、vertex を整理する最後の段階は、単純な数学的処理としてのみ扱われている
DCEL でボロノイ図を保存する
- ボロノイ図は通常 Doubly Connected Edge List(DCEL) で保存される
- DCEL は vertex と edge からなる cell-complex を操作しやすく表現するデータ構造
- edge 中心の表現だが、vertex と face の情報も併せて保持する
- 通常の edge には向きがないが、DCEL では各 edge を双方向の 2 つの half edge として保存する
- ボロノイ図で DCEL に保存される vertex は site ではなく Voronoi Vertex
- edge
E の終点は E.twin.origin で取得でき、右側の face は E.twin.left で取得できる
1件のコメント
Hacker News のコメント
以前、ClojureScriptで Fortune アルゴリズムの進行をアニメーション表示する実装を作ったことがあります: https://voronoi.ajwerner.net/#/app-diagrams
本当に美しいアルゴリズムです
ただ、そのプロジェクト以降は Fortune アルゴリズムが少し嫌いになりました。浮動小数点の数値安定性が良くないからです
点が一直線上にある場合や、浮動小数点の基準でほぼ一直線に近い場合には壊れることがあります
記憶が正しければ、この点では delaunator のほうが優れています: https://github.com/mapbox/delaunator
参考ページには「old」実装へのリンクが見えますが、現在のアニメーション版もオープンソースとして公開する可能性があるのか気になります
数年前に、こういう 3D 可視化を作りました: https://x.com/KangarooPhysics/status/1253336959755251716
uBlock Originで有名な Raymond Hill による JavaScript 実装があります: https://github.com/gorhill/Javascript-Voronoi
ここでは少し手を加えて動くようにしてみました: https://animations.adgent.com/voronoi.html
動画を入力として受け取り、Voronoi 方式で表示するアルゴリズムに通せるのか気になります
そのくらいになると厳密には Voronoi 図ではないかもしれませんが、かなり格好よく見えそうです
D3.jsには新しい実装があります: https://github.com/d3/d3-delaunay
そのページの下のほうに、使われているスイープアルゴリズムの説明と、JavaScript 以外の言語による実装一覧があります
既存の d3-voronoi は廃止予定ですが、こちらで見られます: https://github.com/d3/d3-voronoi
辺に関心がなく、各地点を別々の色で塗るだけでよいなら、シード点から始めるフラッドフィルの変種を使えます
すでにそのピクセルに塗られている色より、その色の距離のほうが短い場合にだけピクセルをスタックに入れればよいです
頂点の上側から 2D 正射影でレンダリングすると、z バッファが最も近い頂点のピクセルを保持します
シェーダーで行う方法もあるでしょうが、古典的な 3D 円錐デモは理解して実装するのが非常に簡単です
D3 が Fortune アルゴリズムから https://mapbox.github.io/delaunator/ に移行したのは興味深いです
理由は「Delaunay 三角形分割や Voronoi 図を作る際に d3-voronoi より5〜10倍高速で、数値的により堅牢であり、Canvas レンダリングを内蔵し、Delaunay グラフの走査や複数の改善点を提供する」ためです
今タイルを計算しているコードは、痛々しいほど単純です
新しい議論: https://github.com/KaliedaRik/Scrawl-canvas/discussions/120
この記事のせいで、Steve が最近どこにいるのか調べることになりました
何十年も前に知り合いだった人です
関連して読む価値のある記事: https://news.ycombinator.com/item?id=37998923 - Fortune アルゴリズムで O(n log n) で Voronoi 図と Delaunay 三角形分割を作る（2020）
以前の記事と議論には、ほかのアルゴリズムの短い要約もあります
個人的には、今でも Jump Flooding Algorithm がいちばん気に入っています: https://en.wikipedia.org/wiki/Jump_flooding_algorithm