ボロノイ図、ボロノイ領域

概要

ボロノイ図関連について、定義はボロノイ図 - Wikipedia参照

高速にボロノイ図を求めるアルゴリズム（フォーチュンのアルゴリズム？）は存在しているのだが、ここでは1つのボロノイ領域をO(n(n+m))くらい?で求める方法を考える。（n:点の数、m:凸領域の頂点の数）

アルゴリズム

全体の領域をRとする。
ボロノイ領域を求めたい点と、他の点に対して、垂直二等分線でRを切断し、注目点を含む方の領域を新たにRとする。
これを全ての点に対して繰り返し行えば、ボロノイ領域を得ることができる。

実装

データ構造、ccw、crosspointは全てSpaghetti Sourceのもの。下のconvex_cutはSpaghettiSourceを少し修正したもの。
convex_cutは凸多角形を直線で切断し左側の部分を残す。
bisectorは、点a,bを結ぶ線分の垂直二等分線を求める。この垂直二等分線の向きはa->bを半時計回りに９０度回転した向き。これで凸多角形を切断すれば、点aを含む方を残すことが出来る。

G convex_cut(const G& g, const L& l) {
  G Q;
  REP(i, g.size()) {
    P A = g[i], B = g[(i+1)%g.size()];
    if (ccw(l[0], l[1], A) != -1) Q.push_back(A);
    if (ccw(l[0], l[1], A)*ccw(l[0], l[1], B) < 0)
      Q.push_back(crosspoint(L(A, B), l));
  }
  return Q;
}
L bisector(P a, P b) {
  P A = (a+b)*P(0.5,0);
  return L(A, A+(b-a)*P(0, PI/2));
}
G voronoi_cell(G g, vector<P> v, int s) {
  REP(i, v.size())
    if (i!=s)
      g = convex_cut(g, bisector(v[s], v[i]));
  return g;
}

もっとスマートな方法があったら教えてください