#include "geometry/polygon.hpp"
#pragma once #include "geometry-base.hpp" using Polygon = vector<Point>; // 多角形の内部に点があるか? // OUT : 0, ON : 1, IN : 2 int contains_polygon(const Polygon &Q, const Point &p) { bool in = false; for (int i = 0; i < (int)Q.size(); i++) { Point a = Q[i] - p, b = Q[(i + 1) % Q.size()] - p; if (imag(a) > imag(b)) swap(a, b); if (sign(imag(a)) <= 0 && 0 < sign(imag(b)) && sign(cross(a, b)) < 0) in = !in; if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1; } return in ? 2 : 0; } // 多角形の面積 Real area(const Polygon &p) { Real A = 0; for (int i = 0; i < (int)p.size(); ++i) { A += cross(p[i], p[(i + 1) % p.size()]); } return A * 0.5; } // 頂点集合から凸包を生成 // boundary : 周上の点も列挙する場合 true template <bool boundary = false> Polygon convex_hull(vector<Point> ps) { sort(begin(ps), end(ps)); ps.erase(unique(begin(ps), end(ps)), end(ps)); int n = ps.size(), k = 0; if (n <= 2) return ps; vector<Point> ch(2 * n); // 反時計周り Real th = boundary ? -EPS : +EPS; for (int i = 0; i < n; ch[k++] = ps[i++]) { while (k >= 2 && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k; } for (int i = n - 2, t = k + 1; i >= 0; ch[k++] = ps[i--]) { while (k >= t && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k; } ch.resize(k - 1); return ch; } // 凸包の内部に点があるか? // OUT : 0, ON : 1, IN : 2 int contains_convex(const Polygon &C, const Point &p) { int N = C.size(); auto b1 = cross(C[1] - C[0], p - C[0]); auto b2 = cross(C[N - 1] - C[0], p - C[0]); if (b1 < -EPS or b2 > EPS) return 0; int L = 1, R = N - 1; while (L + 1 < R) { int M = (L + R) / 2; (cross(p - C[0], C[M] - C[0]) >= 0 ? R : L) = M; } auto v = cross(C[L] - p, C[R] - p); if (equals(v, 0)) { return 1; } else if (v > 0) { return equals(b1, 0) or equals(b2, 0) ? 1 : 2; } else { return 0; } } // 凸包が与えられるので最遠点対を返す // 返り値:頂点番号のペア pair<int, int> convex_polygon_diameter(const Polygon &p) { int N = (int)p.size(); int is = 0, js = 0; for (int i = 1; i < N; i++) { if (imag(p[i]) > imag(p[is])) is = i; if (imag(p[i]) < imag(p[js])) js = i; } Real maxdis = norm(p[is] - p[js]); int maxi, maxj, i, j; i = maxi = is; j = maxj = js; do { if (cross(p[(i + 1) % N] - p[i], p[(j + 1) % N] - p[j]) >= 0) { j = (j + 1) % N; } else { i = (i + 1) % N; } if (norm(p[i] - p[j]) > maxdis) { maxdis = norm(p[i] - p[j]); maxi = i; maxj = j; } } while (i != is || j != js); return minmax(maxi, maxj); }
#line 2 "geometry/polygon.hpp" #line 2 "geometry/geometry-base.hpp" #include <algorithm> #include <cassert> #include <cmath> #include <complex> #include <iostream> #include <vector> using namespace std; using Real = long double; constexpr Real EPS = 1e-10; constexpr Real pi = 3.141592653589793238462643383279L; bool equals(Real a, Real b) { return fabs(b - a) < EPS; } int sign(Real a) { return equals(a, 0) ? 0 : a > 0 ? 1 : -1; } template <typename R> struct PointBase { using P = PointBase; R x, y; PointBase() : x(0), y(0) {} PointBase(R _x, R _y) : x(_x), y(_y) {} template <typename T, typename U> PointBase(const pair<T, U>& p) : x(p.first), y(p.second) {} P operator+(const P& r) const { return {x + r.x, y + r.y}; } P operator-(const P& r) const { return {x - r.x, y - r.y}; } P operator*(R r) const { return {x * r, y * r}; } P operator/(R r) const { return {x / r, y / r}; } P& operator+=(const P& r) { return (*this) = (*this) + r; } P& operator-=(const P& r) { return (*this) = (*this) - r; } P& operator*=(R r) { return (*this) = (*this) * r; } P& operator/=(R r) { return (*this) = (*this) / r; } bool operator<(const P& r) const { return x != r.x ? x < r.x : y < r.y; } bool operator==(const P& r) const { return x == r.x and y == r.y; } bool operator!=(const P& r) const { return !((*this) == r); } P rotate(R rad) const { return {x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad)}; } R real() const { return x; } R imag() const { return y; } friend R real(const P& p) { return p.x; } friend R imag(const P& p) { return p.y; } friend R dot(const P& l, const P& r) { return l.x * r.x + l.y * r.y; } friend R cross(const P& l, const P& r) { return l.x * r.y - l.y * r.x; } friend R abs(const P& p) { return sqrt(p.x * p.x + p.y * p.y); } friend R norm(const P& p) { return p.x * p.x + p.y * p.y; } friend R arg(const P& p) { return atan2(p.y, p.x); } friend istream& operator>>(istream& is, P& p) { R a, b; is >> a >> b; p = P{a, b}; return is; } friend ostream& operator<<(ostream& os, const P& p) { return os << p.x << " " << p.y; } }; using Point = PointBase<Real>; using Points = vector<Point>; // ccw, 点の進行方向 int ccw(const Point& a, const Point& b, const Point& c) { Point x = b - a, y = c - a; if (cross(x, y) > EPS) return +1; // 反時計回り if (cross(x, y) < -EPS) return -1; // 時計回り if (min(norm(x), norm(y)) < EPS * EPS) return 0; // c=a または c=b if (dot(x, y) < EPS) return +2; // c-a-b の順で一直線 if (norm(x) < norm(y)) return -2; // a-b-c の順で一直線 return 0; // a-c-b の順で一直線 } #line 4 "geometry/polygon.hpp" using Polygon = vector<Point>; // 多角形の内部に点があるか? // OUT : 0, ON : 1, IN : 2 int contains_polygon(const Polygon &Q, const Point &p) { bool in = false; for (int i = 0; i < (int)Q.size(); i++) { Point a = Q[i] - p, b = Q[(i + 1) % Q.size()] - p; if (imag(a) > imag(b)) swap(a, b); if (sign(imag(a)) <= 0 && 0 < sign(imag(b)) && sign(cross(a, b)) < 0) in = !in; if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1; } return in ? 2 : 0; } // 多角形の面積 Real area(const Polygon &p) { Real A = 0; for (int i = 0; i < (int)p.size(); ++i) { A += cross(p[i], p[(i + 1) % p.size()]); } return A * 0.5; } // 頂点集合から凸包を生成 // boundary : 周上の点も列挙する場合 true template <bool boundary = false> Polygon convex_hull(vector<Point> ps) { sort(begin(ps), end(ps)); ps.erase(unique(begin(ps), end(ps)), end(ps)); int n = ps.size(), k = 0; if (n <= 2) return ps; vector<Point> ch(2 * n); // 反時計周り Real th = boundary ? -EPS : +EPS; for (int i = 0; i < n; ch[k++] = ps[i++]) { while (k >= 2 && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k; } for (int i = n - 2, t = k + 1; i >= 0; ch[k++] = ps[i--]) { while (k >= t && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k; } ch.resize(k - 1); return ch; } // 凸包の内部に点があるか? // OUT : 0, ON : 1, IN : 2 int contains_convex(const Polygon &C, const Point &p) { int N = C.size(); auto b1 = cross(C[1] - C[0], p - C[0]); auto b2 = cross(C[N - 1] - C[0], p - C[0]); if (b1 < -EPS or b2 > EPS) return 0; int L = 1, R = N - 1; while (L + 1 < R) { int M = (L + R) / 2; (cross(p - C[0], C[M] - C[0]) >= 0 ? R : L) = M; } auto v = cross(C[L] - p, C[R] - p); if (equals(v, 0)) { return 1; } else if (v > 0) { return equals(b1, 0) or equals(b2, 0) ? 1 : 2; } else { return 0; } } // 凸包が与えられるので最遠点対を返す // 返り値:頂点番号のペア pair<int, int> convex_polygon_diameter(const Polygon &p) { int N = (int)p.size(); int is = 0, js = 0; for (int i = 1; i < N; i++) { if (imag(p[i]) > imag(p[is])) is = i; if (imag(p[i]) < imag(p[js])) js = i; } Real maxdis = norm(p[is] - p[js]); int maxi, maxj, i, j; i = maxi = is; j = maxj = js; do { if (cross(p[(i + 1) % N] - p[i], p[(j + 1) % N] - p[j]) >= 0) { j = (j + 1) % N; } else { i = (i + 1) % N; } if (norm(p[i] - p[j]) > maxdis) { maxdis = norm(p[i] - p[j]); maxi = i; maxj = j; } } while (i != is || j != js); return minmax(maxi, maxj); }