|
ひっくり返して線分の集合にバラして双対となる多角形の集合を求めて外周を決めるだけ。
ライブラリ除いて50行くらい。サンプルは通った。
// 最近ライブラリを書きかえてあちこちにバグを埋め込んだのでコピペしないほうが良い #include <algorithm> #include <complex> #include <iomanip> #include <iostream> using namespace std; typedef long double Real; typedef complex<Real> Point; struct Line{Point a, b;}; typedef vector<Point> Polygon; const Real EPS = 1e-8; const Real PI = acos(-1); // 比較関数 inline int sgn(const Real& a, const Real& b = 0) { if (a > b + EPS) return 1; if (a < b - EPS) return -1; return 0; } inline bool near(const Point& a, const Point& b) { return !sgn(abs(a - b)); } namespace std { inline bool operator<(const Point& a, const Point& b) { if (sgn(a.real(), b.real())) return a.real() < b.real(); return a.imag() < b.imag(); } inline bool operator<(const Line& a, const Line& b) { if (!near(a.a, b.a)) return a.a < b.a; return a.b < b.b; } } inline Real dot(const Point& a, const Point& b) { return real(conj(a) * b); } inline Real det(const Point& a, const Point& b) { return imag(conj(a) * b); } // 線分のベクトル inline Point vec(const Line& a) { return a.b - a.a; } // 線分abに対する点cの位置 enum CCW{FRONT = 1, RIGHT = 2, BACK = 4, LEFT = 8, ON = 16}; inline int ccw(const Point& a, const Point& b, const Point& c) { if (near(a, c) || near(b, c)) return ON; int s = sgn(det(b - a, c - a)); if (s) return s > 0 ? LEFT : RIGHT; return (a < b) == (b < c) ? FRONT : (c < a) == (a < b) ? BACK : ON; } // 有向角度 inline Real arg(const Point& base, const Point& a, const Point& b) { return arg((b - base) / (a - base)); } // 射影 inline Point proj(const Point& a, const Point& b) { return a * dot(a, b) / norm(a); } inline Point perp(const Line& l, const Point& p) { return l.a + proj(vec(l), p - l.a); } inline Point refl(const Line& l, const Point& p) { return perp(l, p) * (Real)2 - p; } // 交差判定 inline bool iLL(const Line& a, const Line& b, bool strict = false) { if (sgn(det(vec(a), vec(b)))) return true; return !strict && !sgn(det(vec(a), b.a - a.a)); } inline bool iLS(const Line& a, const Line& b, bool strict = false) { if (strict) return sgn(det(vec(a), b.a - a.a)) * sgn(det(vec(a), b.b - a.a)) < 0; return sgn(det(vec(a), b.a - a.a)) * sgn(det(vec(a), b.b - a.a)) <= 0; } inline bool iSS(const Line& a, const Line& b, bool strict = false) { int cwa = ccw(a.a, a.b, b.a) | ccw(a.a, a.b, b.b); int cwb = ccw(b.a, b.b, a.a) | ccw(b.a, b.b, a.b); if ((cwa & cwb) == (LEFT | RIGHT)) return true; return !strict && ((cwa | cwb) & ON); } // 交点 inline Point pLL(const Line& a, const Line& b) {return a.a + vec(a) * (det(vec(b), b.a - a.a) / det(vec(b), vec(a)));} // 多角形の面積 inline Real aPol(const Polygon& vp) { Real res = 0; for (int i = 0; i < (int)vp.size(); ++i) res += det(vp[i], vp[(i + 1) % vp.size()]); return res / 2; } // 凸多角形カット inline Polygon convexCut(const Polygon& p, const Line& l) { vector<Point> q; for (int i = 0; i < (int)p.size(); ++i) { if (ccw(l.a, l.b, p[i]) != RIGHT) q.emplace_back(p[i]); Line s = {p[i], p[(i + 1) % p.size()]}; if (iLS(l, s, true)) q.emplace_back(pLL(l, s)); } return q; } // 線分をマージする inline void merge(vector<Line>& vs) { for (Line& s : vs) if (s.b < s.a) swap(s.a, s.b); sort(vs.begin(), vs.end()); for (int i = 0; i < (int)vs.size(); ++i) { for (int j = 0; j < i; ++j) { if (iSS(vs[i], vs[j]) && !iLL(vs[i], vs[j], true)) { vs[j].b = max(vs[i].b, vs[j].b); vs.erase(vs.begin() + i--); break; } } } } // 線分アレンジメント 隣の点への辺のみを持つ inline void sArr(vector<Line>& vs, vector<Point>& vp, vector<vector<int>>& g) { merge(vs); for (const Line& s : vs) { vp.emplace_back(s.a); vp.emplace_back(s.b); } for (int i = 0; i < (int)vs.size(); ++i) { for (int j = 0; j < i; ++j) { if (iSS(vs[i], vs[j])) vp.emplace_back(pLL(vs[i], vs[j])); } } sort(vp.begin(), vp.end()); vp.erase(unique(vp.begin(), vp.end(), near), vp.end()); g = vector<vector<int>>(vp.size()); for (const Line& s : vs) { vector<pair<Real, int>> v; for (int j = 0; j < (int)vp.size(); ++j) { if (ccw(s.a, s.b, vp[j]) == ON) v.emplace_back(make_pair(norm(vp[j] - s.a), j)); } sort(v.begin(), v.end()); for (int j = 0; j < (int)v.size() - 1; ++j) { g[v[j + 1].second].emplace_back(v[j].second); g[v[j].second].emplace_back(v[j + 1].second); } } } // 線分の集合の双対である多角形の集合を返す inline vector<Polygon> dual(vector<Line>& s) { vector<Point> vp; vector<vector<int>> edge; sArr(s, vp, edge); vector<vector<bool>> used(edge.size()); for (int i = 0; i < (int)edge.size(); ++i) used[i] = vector<bool>(edge[i].size(), false); vector<Polygon> vpol; for (int i = 0; i < (int)edge.size(); ++i) for (int j = 0; j < (int)edge[i].size(); ++j) if (!used[i][j]) { Polygon pol; pol.emplace_back(vp[i]); used[i][j] = true; int pre = i, now = edge[i][j]; while (now != i) { pol.emplace_back(vp[now]); Real a = -4; int next = -1; for (int k = 0; k < (int)edge[now].size(); ++k) { Real aa = arg(vp[now], (Real)2 * vp[now] - vp[pre], vp[edge[now][k]]); if (!sgn(aa, PI)) aa = -PI; if (a < aa) { a = aa; next = k; } } used[now][next] = true; pre = now; now = edge[pre][next]; } vpol.emplace_back(pol); } return vpol; } // ここまでライブラリ Real len(const Polygon& pol) { Real res = 0; for (int i = 0; i < (int)pol.size(); ++i) { res += abs(pol[i] - pol[(i + 1) % pol.size()]); } return res; } Real size(const Polygon& pol) { int res = pol.size(); for (int i = 0; i < (int)pol.size(); ++i) { if (ccw(pol[i], pol[(i + 2) % pol.size()], pol[(i + 1) % pol.size()]) == ON) --res; } return res; } int main() { while (true) { int n; cin >> n; if (n == 0) break; Polygon pol(n); for (auto& p : pol) { Real x, y; cin >> x >> y; p = Point(x, y); } Polygon res; for (int i = 0; i < n; ++i) { for (int j = i + 1; j < n; ++j) { Point mid = (pol[i] + pol[j]) / (Real)2; Point mid2 = mid + (pol[j] - pol[i]) * Point(0, 1); Polygon pol1 = convexCut(pol, (Line){mid, mid2}); Polygon pol2 = convexCut(pol, (Line){mid2, mid}); for (auto& p : pol2) p = refl((Line){mid2, mid}, p); vector<Line> vl; for (int i = 0; i < (int)pol1.size(); ++i) vl.emplace_back((Line){pol1[i], pol1[(i + 1) % pol1.size()]}); for (int i = 0; i < (int)pol2.size(); ++i) vl.emplace_back((Line){pol2[i], pol2[(i + 1) % pol2.size()]}); auto d = dual(vl); for (auto pp : d) { if (sgn(aPol(pp)) >= 0) continue; if (size(res) < size(pp) || (size(res) == size(pp) && len(res) < len(pp))) res = pp; } } } cout << fixed << setprecision(6) << len(res) << endl; } }
