🌼
【ボロノイ分割】GISを作るデータ構造とアルゴリズム【計算幾何学】

2025/11/20に公開
 概要前回の記事では, GISの基本機能のひとつ, 地図の重ね合わせを説明した. 主要なテクニックとして, 平面走査法を用いた効率的な線分交差判定アルゴリズムを導入した.
今回は, GISの基本機能のひとつ, ボロノイ分割を説明する. 経済圏や商圏のエリア推定に利用するなど, 複数点間で均等となる境界線を引く領域分割法である.
今回も効率的なアルゴリズムを紹介するので, 平面走査法については過去の記事などを参考に,イメージをまず頭に入れておくことをお勧めする.
今回も「コンピュータ・ジオメトリ　第3版　計算幾何学：アルゴリズムと応用」や, レンセラー工科大学の講義資料を参考に説明する.
また, ボロノイ分割に関連する問題は,最適配置の数理という本にとてもわかりやすく整理されていいるのでこちらも参照していく.

 郵便ポスト配置問題まず, ボロノイ分割は何を解決したいものか？をイメージしよう. 代表的な問題が, 「郵便局問題」 である.[1] 次の例を考えよう.
手紙やハガキを書いたら近くの郵便局（ポスト）に出しに行くという経験をしたことはあるだろう. ある街で, 有限個の郵便局（ポスト）を設置するとしたら, どこに配置すると地域の人々にとって'便利'だろう...そんな郵便局（ポスト）の配置計画を考える.
ここで仮定をおく,
各郵便局（ポスト）のサービスは平等である. つまり, どこで手紙を出そうか郵送サービスには違いが出ない.
人は最も近い郵便局（ポスト）に手紙を出しに行くものだ.
2.の近さは「直線距離」で決定されている.
上記の仮定によれば, 人は自身との直線距離が最も小さい, 最近傍の郵便局（ポスト）に手紙を出す,という行動原理をもつ.
この前提のもと, 例えば, 全住民からの移動距離の合計が小さくなるように, 決められた数の郵便局（ポスト）の最適配置計画を考えるのが,この 「郵便局問題」 である.
いきなり最適化...といわれても難しいと思うので, 例えば上記の配置例での各郵便局（ポスト）が誘引できるエリアを作成してみる.
すると上記のような領域分割ができる. ここで, 郵便局（ポスト）P_5を真ん中にもつオレンジ色の領域をみる. これは, この中にいる住民は他の郵便局（ポスト）よりも, 郵便局（ポスト）P_5が最近傍となる領域を示している. これが郵便局（ポスト）P_5のボロノイ領域である.
各郵便局（ポスト）ごとにボロノイ領域を作成し, 各領域内に含まれる住民（の位置）をすべて特定する. その後, 住民ごとの郵便局（ポスト）までの移動距離を合計すれば, その配置計画の'便利さ'を定量的に求められる.
実際には,「郵便局問題」 を解くには,別途住民ごとの郵便局（ポスト）までの移動距離を最適化する計算式が必要であるが, このボロノイ領域を作成はその解に必要な処理となっている.
まとめると, ボロノイ分割とは, 複数の点,それぞれを中心とした領域を求めることを目的に, 他の点と均等距離を保った領域に分割する処理である.

 ボロノイ分割ボロノイ分割のイメージをもったところで, このボロノイ分割の性質を見ていこう.
もう少し一般化した図を以下におく.
前述した通り, ボロノイセルとは, この中の点（住民）は他の郵便局（ポスト）よりも, この領域の中心となる郵便局（ポスト）が最近傍となる領域を示している.(図中のオレンジ領域）
これを数式を使って定義する. つまり, 郵便局（ポスト）をP_iとし, そのボロノイ領域をV(P_i)とおく.いま, ある点pに対し, ボロノイセルV(P_i)の中では,

d(p,P_i) \leq d(p, P_j) (i \neq j)
が成り立っている. ここで, d(p,q)は2点間p,qのユークリッド距離である.

図中でいえば, ボロノイセルV(P_5)の中では,d(p,P_5) \leq d(p, P_1)であるし,  d(p,P_5) \leq d(p, P_4)である.
そのような性質をもつボロノイ領域を求めるのがボロノイ分割であるが, どう求めようか...
準備として, 図形は線分の集合で定義できることを思い出そう. ここで, ボロノイセルの境界部分を形成する線分をボロノイ辺,そして, ボロノイセルの頂点をボロノイ点と呼んでいこう. すると, 例えばボロノイ辺が何者かは想像できるだろう.

ボロノイ辺は, あるふたつの点の二等分線の一部である.
二等分線の一部というのに着目してほしい. 実際,ボロノイ辺が完全な直線になることはなく, 線分や半直線である.
では, 次に あるふたつの点の二等分線をどう分断していくかを考える.ここで新たな言葉として, 半平面を定義する. これは, あるふたつの点を, それらの二等分線で分断したうちの片方であり, Hで定義する. ここで, 2点A,Bの二等分線から半平面Hのうち, 点Aを含む方をH(A,B),点Bを含む方をH(B,A)と区別する.
さきの例をこの方法でH(P_5,P_1)を求めると下のピンク色の領域となる.
この半平面Hを使うと, ボロノイ領域は次のように定義できる.
ボロノイセルP_iによるボロノイセルV(P_i)は, 他の点との二等分線で生まれる半平面H(P_i,P_j) (i \neq j)の共通部分である
図の例だと, V(P_5) = H(P_5,P_1) \bigcap H(P_5,P_2) \bigcap H(P_5,P_3) \bigcap H(P_5,P_4)である.
なるほど, ボロノイセルP_5と他点の二等分線間の交点を求めて, それらを結ぶ線分のうち, ボロノイセルP_5を囲む最小領域を探せばよいと...となる.
ここで, 脳筋スタイルのBrute Force法を考える. つまり, 上記の考えをそのまま総当たりですべての点ペアで実施する. わかっていると思うが, 計算時間はO(n^2)のオーダーとなり非効率.
さてどう改善したものか...
!
 ボロノイ領域のその他の性質補足として, 他の面白い性質も簡単に述べておく.

まず, 各ボロノイ点は,特別な場合を除いて, ３つのボロノイ領域で共有されている. 特別な場合というのは,例えば, すべてのボロノイセルが格子状に並んでいる場合である. この場合は, 4つのボロノイ領域で共有されるが, それ以外のケースでは３つのボロノイ領域で共有されている.
このことから次のことが自動的にわかる.
あるボロノイ領域のボロノイ点を中心として,その2点を半径とした円を描くと, 3つのボロノイセルを通る. このとき,その円の内側には他のボロノイセルは必ず入らない
例として以下の図を見てほしい. ボロノイ領域V(P_5)のボロノイ点pを中心とした円は, ボロノイセルP_1,P_2,P_5を通るが, それ以外の点を円周上および内部にもたない.
この性質はドロネー三角分割という処理で重要になってくるが,今回は使用しないので, 参考程度に受け止めてほしい.

 平面走査法によるボロノイ図作成総当たり法の問題点はどこか, ずばり, 不要な線分間（直線間）の交点も求めていることだ.

ボロノイセルは, 近くにあるボロノイサイト同士で境界が構築される. つまり, 線分間（直線間）の交点を求めるにしても, 近くにあるボロノイサイト同士だけを考慮すればよい, とまず考える.
とすれば, 以前紹介した平面走査法が利用できそうだ.
詳細は省くが, 平面走査法では,走査線と呼ばれるラインをy軸の大きい方から小さい方へと移動させながら, 逐次的に入力データを処理していく.

その特徴として,
ある走査線の位置における状態構造Tを蓄える.
状態構造Tの更新タイミングは, イベント点と呼ばれる場所に走査線が到達した場合である
イベント点はイベントキューQで一元管理される

 線分交差判定との違い以前紹介した線分交差処理では, イベントキューQには, 線分の端点および交点が蓄えられ,状態構造Tには, 走査線と交わる線分間の隣接関係を管理していた.
では, ボロノイ図作成の場合はどうなるか🧐

例えば, 走査線と交わるボロノイセルに着目し, イベントキューQにはボロノイサイト, 状態構造Tには, そのボロノイサイト間の二等分線の隣接関係を管理すればいいのだろうか.
残念ながら, この方法は有効ではない. その理由として, 走査線よりも上のボロノイサイトの情報だけでは, イベント点PのボロノイセルVor(P)の形状が決定できないからだ. 逆にいえば, まだ走査していない走査線よりも下のボロノイサイトも知らないといけないが, 平面走査法では考慮できない.
つまり, ボロノイ図（ボロノイセル）を求めるために必要な頂点を探すために必要な情報が,走査線より上の情報だけでは足りないのだ(下図参考）
では,どうするか, というと, 走査線よりも下のボロノイサイトによる影響を受けない, 走査線よりも上のボロノイサイトに関するセルの部分だけを管理する
🧐🤔

 弧に着目する考え方を変えてみる. いま, あるボロノイサイトp_j周辺の領域のうち,走査線lよりも下にあるボロノイサイトの影響を受けないのはどこか.
ここで, 走査線lより上の閉半平面をl^{+}とかく. ある点q \in l^{+}が, 走査線lより下のボロノイサイトよりも, ボロノイサイトp_jに近くになる領域はどこか考える.
実は, この範囲はd(p_j, q) = d(q,l^{+})となる放物線で描けてしまう. （下図　左）

つまり, この放物線を境界として上側にあるある点q \in l^{+}は,走査線lより下のどのボロノイサイトよりも,  ボロノイサイトp_jに近くになる. つまり, このような範囲はボロノイセルVor(p_j)であり,そして, 走査線lより下のどのボロノイサイトに関わらず決定される.
この放物線を弧(Arc) と呼び, ボロノイサイトごとに作成される. 走査線に近い弧(Arc)の無武運を結んでできるラインを, ビーチラインとよぶ.
このビーチラインより上の部分は, 走査線lより下のどのボロノイサイトに関わらず決定されるため, 先述した管理対象になりうる.
いうまでもなく, このビーチラインの形状は走査線の移動により変化する.そして面白いことに, 各弧の交点は必ずボロノイ辺のうえを辿る.[2] (上図　右）
この交点の名前をブレークポイント（Break Point) とよぶ.このブレークポイントを追跡すれば, 走査線lより下の情報を気にせずとも, ボロノイ辺,つまりは,ボロノイセルが構築できそうだ.

 イベントキューと走査線状態構造
 Fortuneのアルゴリズムということで, 平面走査法によるボロノイ図作成では, このビーチライン,つまり, 各ボロノイサイトと走査線により決定される弧(arc)を状態として管理すればよい.

そして, ボロノイセルを決定していくには, 各弧(arc)での交点（ブレークポイント）に着目する.
これをアルゴリズムとして整理したものが, Fortuneのアルゴリズムである.

アルゴリズムの説明前に, 処理の雰囲気を以下の動画で確認するとよい. 走査線を移動させると,ビーチラインが変化,そして, ボロノイ辺/ボロノイセルが決定されていく過程がみれる.
https://www.youtube.com/watch?v=k2P9yWSMaXE
FortuneのアルゴリズムをPythonで忠実に再現したコードが有志により作成されている.

なので, 具体的な実装を確認したい場合は以下を参照してほしい.

https://github.com/Yatoom/foronoi/tree/master
以後, コードを解読するうえで重要な部分だけ詳細に説明する.

 イベントの管理方法**各ボロノイサイトと走査線により決定される弧(arc)を状態として管理すればよい.**と述べた.

このとき, 弧(arc)同士の隣接関係を状態構造としてツリーに記録する.[3]
しかし, 動画をみるとわかるが, 初期状態では弧はないし,処理の最後も弧は残らない. つまり, 走査線を移動させる過程で, 弧が新規に生まれたり,あるいは,消滅したりする.
弧同士の交点は, ボロノイ辺/ボロノイセルを決定する上で必要な情報であるため, 弧の発生・消滅のタイミングを知っておくことが重要である.それぞれ見ていこう.

 弧が生まれるタイミングまず弧が新たに発生するときだが,ずばり, 走査線があるボロノイサイトにかかったタイミングである. このとき, このボロノイサイトと走査線の関係性から, このボロノイサイト由来の新規の放物線が生まれる. つまり,間違いなくボロノイサイト群はイベントキューに蓄えるべきイベント点である.(これをSite Eventと呼ぶ）
同時に,ビーチラインも変化することに注意だ. 新規の弧の発生により,その前にできていたビーチラインに変化を与える. ビーチラインが変わるということは, 状態構造も変化するということで重要な情報である.
では, この変化をどのように表現するのか考える.
講義資料より
まず,そもそも新規の弧がビーチラインのどの箇所に影響を与えるのか, 厳密には,ビーチラインを構成するどの弧に影響を与えるのか求める. これは簡単で, 走査線上のボロノイサイトからビーチラインに向かって垂線を引いた時に交差する弧を求めたら良い.　（上図　真ん中）
そして, (動画をよく見ると気づくと思うが)新規の弧は以後この影響を与える弧を左右に分断するように広がっていく. （上図　右）すなわち. 新規の弧の発生により, 分断された左右の弧に由来する隣接関係, つまりは状態構造を更新しなければいけない.
また, 分断された左右の弧と新規の弧とのブレークポイントを結んだエッジを管理する. なぜならこれがボロノイ辺/ボロノイセルを決定するからである.

 弧が消滅するタイミング次に弧が消滅するタイミングを考える. 残念ながら, "弧が生まれるタイミング"と同じようにわかりやすいイベント点があるわけではない.
そこで, 弧が消滅するタイミングに何が起こるのか見てみる.
講義資料より
上図の例で説明する. いま, ボロノイサイトp_i, p_j, p_kに由来する弧\alpha,\alpha',\alpha''でビーチラインができている. 走査線がさらに移動すると, 弧\alpha'が縮小し,次第に消えてしまう.
実は,弧が消滅するその瞬間, その弧の左右にあった弧同士で新たなブレークポイント（点q）をもつ. ここでこの点に着目すると, ボロノイサイトp_i, p_j, p_kと点qの距離は等しく,そして, 点qと走査線との距離も等しい.
つまり, 点qを中心にボロノイサイトp_i, p_j, p_kを外周にもつ円を描け,そして, この円の下端と走査線は接する. 逆にいえば, ボロノイサイトp_i, p_j, p_kを外周にもつ円の下端点に走査線が到達したとき,真ん中にある弧が消滅する.
これで弧が消滅するタイミングもみえた. 実はそれぞれの弧が消滅するタイミングは, 弧が生まれるタイミングに合わせてさきに計算できる. なぜなら, ある弧が発生したタイミングで, 分断される弧はまず消えゆく運命にあるが, それがいつ消えるのかは,  分断された弧の両隣,それぞれ３つのボロノイサイトを通る円を描いたらわかるからである.
よって, 弧が生まれるタイミングに合わせて, 分断された弧が消えるタイミングをイベント点としてイベントキューに記録できる（これを,Circle Eventとよぶ）
このイベント点に走査線が到達したとき,当該の弧を消滅,つまり, 状態構造から排除すればよい.

同時に, このタイミングにできる交点qは,ボロノイ頂点になることにも注意がいる.

 処理イメージ最後に簡単に上記の処理イメージを共有しておく.


いまボロノイサイトPに走査線がかかったばかりの状況を考える.ボロノイサイトP由来の弧を\hat{P}とする. この弧\hat{P}は, 点Pから上に垂線を引くと弧\hat{B}と交差する.この交差により, 弧\hat{B}は左右に\hat{B}',\hat{B}''と分断される.
よって, 状態構造に新規の弧\hat{P}を加えつつ, 既存の弧の左右関係も更新する.


さらに走査線が移動すると, 弧\hat{P}が広がるとともに,弧\hat{B}を分断後にできた弧\hat{B}',\hat{B}''が次第に小さくなる. よって, この二つの消滅タイミング(Circle Event)をここで求めておく.
例えば, 弧\hat{B}'は, 弧\hat{A}と弧\hat{P}によって小さくなるため, それらを形成するボロノイサイトA,B,Pを通る円を求め,その下端をイベント点としてイベントキューに入れる.


このイベント点を走査線が通るとき,弧\hat{B}'は消滅するため, 状態構造から弧\hat{B}'を加えつつ, 既存の弧の左右関係も更新する.
基本的な流れは上記の通りである. もちろん, ある弧が消滅するタイミングであるCircle Eventは, 走査線の移動ともに変更される場合がある.[4]　そのため, Circle Eventも更新されうるということには注意が必要である.
そして, ブレークポイントを結んだ二重辺を記録しておくことも忘れてはいけない.

上図のように, 分断された弧と新規の弧とのブレークポイントは,各弧の元のボロノイサイトのボロノイセルの辺になっていく. そのため,このブレークポイントを端点とした二重辺も一緒に記録しておくと,最終的に求めたいボロノイセルの面も作成できる.
Fortuneのアルゴリズムの要点をまとめると,
イベントキューには,弧が生まれる"Site Event"と, 弧が消える"Circle Event"を格納
状態構造には, 弧同士の左右関係を記憶している
弧同士のブレークポイントをつないだものを二重辺として記録・更新していく

 まとめ今回は, 「ボロノイ分割」を説明した.

平面走査法＋弧の管理(+二重辺リスト)という組み合わせで効率的にボロノイ図を作成するFortuneのアルゴリズムを見てきた. このテーマでは,有志による動画や実装例もあるので理解はしやすいと思う.
ただ, 今回の記事が少しでも理解の助けになれば幸いである.

脚注
日本だと,郵便ポスト問題のほうが馴染みがあるかも ↩︎
孤は,それぞれのボロノイサイトと走査線との距離が等しい距離を結んだものである. 孤の交点では, それぞれの弧に対応するボロノイサイトと走査線との距離,つまり, 各ボロノイサイトと交点も常に等しくなるため ↩︎
線分交差の時と同じで, ある孤の左右だけみれば交点がわかるからだ. 左右ということは二分木である. ↩︎
新規の弧により,左右の弧がかわった場合など ↩︎