AOJ 2385 - Shelter
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2385
問題
頂点数Mの多角形領域の中に,N個のシェルターがある。
多角形領域の中に等確率で位置すると考えたとき、最も近いシェルターまでの距離の期待値を求めよ。
制約
M,N<=100
-1000<=多角形の頂点座標<=1000
解法
あるシェルターが最も近くなる領域は、そのシェルターを母点とするボロノイ領域である。
そこで、各ボロノイ領域で期待値を計算して(重み付き平均として)足し合わせる。
ここでさらに、各ボロノイ領域を、母点とボロノイ領域の頂点で作られる三角形に分割する。
三角形に分割した後、下図のように、母点が原点、残りの2点のx座標が一致するように平行移動、回転する。
この三角形領域上で原点からの距離を積分する。
積分を計算すると、
のようになる。(追記:上の図のlとkの定義逆)
後は実装するだけ。
あと解説に”簡単な積分”と書いてあったけど自分はどう積分するかで小一時間悩んだorz。
解説を見るとグリーンの定理とか書いてあるけどどう使うんだろう…
ソース
#include <iostream> #include <vector> #include <algorithm> #include <complex> #include <cstdio> #include <cassert> using namespace std; #define REP(i,n) for(int i=0;i<(int)n;++i) #define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i) #define ALL(c) (c).begin(), (c).end()) const double EPS = 1e-8; const double PI = acos(-1); typedef complex<double> P; namespace std { bool operator < (const P& a, const P& b) { // if (abs(a-b)<EPS) return 0; return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b); } } double cross(const P& a, const P& b) { return imag(conj(a)*b); } double dot(const P& a, const P& b) { return real(conj(a)*b); } struct L : public vector<P> { L(const P &a, const P &b) { push_back(a); push_back(b); } L() {} }; typedef vector<P> G; #define curr(P, i) P[i] #define next(P, i) P[(i+1)%P.size()] int ccw(P a, P b, P c) { b -= a; c -= a; if (cross(b, c) > 0) return +1; // counter clockwise if (cross(b, c) < 0) return -1; // clockwise if (dot(b, c) < 0) return +2; // c--a--b on line if (norm(b) < norm(c)) return -2; // a--b--c on line return 0; } P crosspoint(const L &l, const L &m) { double A = cross(l[1] - l[0], m[1] - m[0]); double B = cross(l[1] - l[0], l[1] - m[0]); if (abs(A) < EPS && abs(B) < EPS) return m[0]; // same line if (abs(A) < EPS) assert(false); // !!!PRECONDITION NOT SATISFIED!!! return m[0] + B / A * (m[1] - m[0]); } double area(const G& g) { double A = 0; for (int i = 0; i < g.size(); ++i) { A += cross(g[i], next(g, i)); } return abs(A/2); } G convex_cut(const G& g, const L& l) { G Q; REP(i, g.size()) { P A = curr(g, i), B = next(g, i); 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(const P &a, const 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, const vector<P> &v, int s) { REP(i, v.size()) if (i!=s) g = convex_cut(g, bisector(v[s], v[i])); return g; } P rotate(P p, double ang) { return p * P(cos(ang), sin(ang)); } double calc(P a, P b, P c) { // aを原点とし,bとcのx座標が同一になるように変換する b -= a; c -= a; double ang = PI/2 - arg(c-b); b = rotate(b,ang); c = rotate(c,ang); double p = b.real(); double q = b.imag(); double r = c.imag(); double k = q/p; double l = r/p; double A = (l+l*l*l/3) - (k+k*k*k/3); return p*p*p*p * A / 4; } int main() { int m,n; cin>>m>>n; G g(m); REP(i,m) cin>>g[i].real()>>g[i].imag(); vector<P> v(n); REP(i,n) cin>>v[i].real()>>v[i].imag(); double ans = 0; REP(i,n) { G voronoi = voronoi_cell(g, v, i); REP(j,voronoi.size()) { ans += calc(v[i],curr(voronoi,j),next(voronoi,j)); } } ans /= area(g); printf("%.10f\n", ans); }