多角形
(geometry/polygon.hpp)

View this file on GitHub
Last update: 2026-08-11 23:01:09+09:00
Include: #include "geometry/polygon.hpp"

二次元平面上の多角形に関する包含判定，凸包，直径，切断を扱う．

頂点を反時計回りに並べた二次元多角形を Polygon で表す．

contains_polygon(P, p)：単純多角形 P と点 p の位置関係を返す．外部，周上，内部ならそれぞれ $0,1,2$．P は凸でなくてもよい．$O( P )$ 時間．

area(P)：多角形 P の符号付き面積を返す．$O( P )$ 時間．

is_convex(P)：反時計回りの多角形 P が凸か判定する．周上に連続する三点が一直線上にある場合も許す．$O( P )$ 時間．

convex_hull<boundary>(ps)：点集合 ps の凸包を反時計回りに返す．boundary = true なら周上の点をすべて含める．$O( ps \log ps )$ 時間．

contains_convex(C, p)：反時計回りの凸多角形 C と点 p の位置関係を返す．外部，周上，内部ならそれぞれ $0,1,2$．$O(\log C )$ 時間．

convex_polygon_diameter(C)：反時計回りの凸多角形 C の最遠点対を頂点番号の組で返す．$O( C )$ 時間．

argmax_dot(C, p)：反時計回りの狭義凸多角形 C に対し，dot(C[i], p) を最大化する最小の頂点番号 i を返す．p は零ベクトルでないものとする．$O(\log C )$ 時間．

convex_polygon_cut_info(C, a, b)：反時計回りの狭義凸多角形 C と有向直線 a $\to$ b の左閉半平面との共通部分を ConvexPolygonCutResult で返す．$O(\log C )$ 時間．

ConvexPolygonCutResult::to_polygon(C)：切断結果を反時計回りの Polygon として構築する．出力頂点数を $K$ として $O(K)$ 時間．
convex_polygon_cut(C, a, b)：convex_polygon_cut_info(C, a, b).to_polygon(C) を返す．$O(\log C +K)$ 時間．

ConvexPolygonCutResult の [first, last) は，切断後に残る元の頂点 C[i % C.size()] の巡回区間を表す．front, back はそれぞれ区間の前後に挿入する交点であり，交点が元の頂点に一致する場合は nullopt となる．

資料

geometry/polygon.hpp - Nyaan’s Library

Depends on

二次元幾何の基本要素 (geometry/geometry-base.hpp)

Required by

Verified with

Code

#pragma once

#include "geometry/geometry-base.hpp"

// 反時計周り
using Polygon = vector<Point>;

// 多角形の内部に点があるか
// OUT : 0, ON : 1, IN : 2
// 凸でなくてもよい, 自己交差はもたない
// crossing number algorithm : x 軸正方向に伸ばした半直線が奇数回交われば内部
int contains_polygon(const Polygon& P, const Point& p) {
  bool in = false;
  for (int i = 0; i < (int)P.size(); i++) {
    Point a = P[i] - p, b = P[(i + 1) % P.size()] - p;
    if (a.y > b.y) swap(a, b);
    if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1;
    if (sign(a.y) <= 0 && 0 < sign(b.y) && sign(cross(a, b)) < 0)
      in = !in;
  }
  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;
}

bool is_convex(const Polygon& P) {
  int n = P.size();
  if (n < 3) return true;
  for (int i = 0; i < n; ++i) {
    if (ccw(P[i], P[(i + 1) % n], P[(i + 2) % n]) == -1) return false;
  }
  return true;
}

// 頂点集合から凸包を生成
// 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;
  if constexpr (boundary) {
    bool collinear = true;
    for (int i = 2; i < n; ++i) {
      if (!equals(cross(ps[1] - ps[0], ps[i] - ps[0]), 0)) {
        collinear = false;
        break;
      }
    }
    if (collinear) return ps;
  }
  vector<Point> ch(2 * n);
  // 反時計周り
  const 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();
  if (N == 0) return 0;
  if (N == 1) return abs(C[0] - p) < EPS ? 1 : 0;
  if (N == 2) return ccw(C[0], C[1], p) == 0 ? 1 : 0;
  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& C) {
  int N = (int)C.size();
  assert(N > 0);
  if (N == 1) return {0, 0};
  if (N == 2) return {0, 1};
  int is = 0, js = 0;
  for (int i = 1; i < N; i++) {
    if (C[i].y > C[is].y) is = i;
    if (C[i].y < C[js].y) js = i;
  }
  Real maxdis = norm(C[is] - C[js]);

  int maxi, maxj, i, j;
  i = maxi = is;
  j = maxj = js;
  do {
    if (cross(C[(i + 1) % N] - C[i], C[(j + 1) % N] - C[j]) >= 0) {
      j = (j + 1) % N;
    } else {
      i = (i + 1) % N;
    }
    if (norm(C[i] - C[j]) > maxdis) {
      maxdis = norm(C[i] - C[j]);
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return minmax(maxi, maxj);
}

// min argmax_i dot(C[i], p)
int argmax_dot(const Polygon& C, const Point& p) {
  int n = C.size();
  assert(n > 0);
  assert(p.x != 0 || p.y != 0);
  if (n <= 2) {
    if (n == 1 || dot(C[0], p) >= dot(C[1], p)) return 0;
    return 1;
  }
  auto arg_half = [](const Point& a) {
    return a.y < 0 || (a.y == 0 && a.x < 0);
  };
  auto arg_less = [&](const Point& a, const Point& b) {
    int ha = arg_half(a), hb = arg_half(b);
    return ha != hb ? ha < hb : cross(a, b) > 0;
  };
  auto edge = [&](int i) { return C[(i + 1) % n] - C[i]; };
  int l = 0, r = n - 1;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge(r), edge(m)))
      l = m + 1;
    else
      r = m;
  }
  int s = l;
  Point q = p.rotate90();
  l = 0, r = n;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge((s + m) % n), q))
      l = m + 1;
    else
      r = m;
  }
  int i = (s + (l == n ? 0 : l)) % n;
  Point e = edge(i);
  if (!arg_less(e, q) && !arg_less(q, e)) {
    int low = l;
    l = low, r = n;
    while (l < r) {
      int m = (l + r) / 2;
      if (!arg_less(q, edge((s + m) % n)))
        l = m + 1;
      else
        r = m;
    }
    int count = l - low;
    return i + count >= n ? 0 : i;
  }
  return i;
}

struct ConvexPolygonCutResult {
  int first = 0, last = 0;
  optional<Point> front, back;

  bool empty() const { return first == last && !front && !back; }
  Polygon to_polygon(const Polygon& C) const {
    Polygon ret;
    auto push = [&](const Point& p) { if (ret.empty() || abs(ret.back() - p) >= EPS) ret.push_back(p); };
    if (front) push(*front);
    for (int i = first; i < last; ++i) push(C[i % C.size()]);
    if (back) push(*back);
    if (ret.size() >= 2 && abs(ret.front() - ret.back()) < EPS) ret.pop_back();
    return ret;
  }
};

ConvexPolygonCutResult convex_polygon_cut_info(const Polygon& C, const Point& a, const Point& b) {
  int n = C.size();
  assert(n > 0 && a != b);
  Point d = b - a, normal = d.rotate90();
  auto side = [&](int i) { return sign(cross(d, C[i % n] - a)); };
  int mx = argmax_dot(C, normal), mn = argmax_dot(C, -normal);
  if (side(mx) < 0) return {};
  if (side(mn) >= 0) return {0, n, nullopt, nullopt};

  int l = mn, r = mx;
  if (r < l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) >= 0 ? r : l) = m;
  }
  int first = r;
  l = mx, r = mn;
  if (r <= l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) < 0 ? r : l) = m;
  }
  int last = r;
  if (last < first) last += n;
  int count = last - first;
  first %= n;
  ConvexPolygonCutResult ret{first, first + count, nullopt, nullopt};
  auto cross_point = [&](const Point& p, const Point& q) {
    return p + (q - p) * (cross(d, a - p) / cross(d, q - p));
  };
  if (side(first) > 0) ret.front = cross_point(C[(first + n - 1) % n], C[first]);
  int back = (ret.last - 1) % n;
  if (side(back) > 0) ret.back = cross_point(C[back], C[(back + 1) % n]);
  return ret;
}

Polygon convex_polygon_cut(const Polygon& C, const Point& a, const Point& b) {
  return convex_polygon_cut_info(C, a, b).to_polygon(C);
}

/**
 * @brief 多角形
 * @docs docs/geometry/polygon.md
 */

#line 2 "geometry/polygon.hpp"

#line 2 "geometry/geometry-base.hpp"

#include <bits/stdc++.h>

using Real = long double;
constexpr Real EPS = 1e-10;
constexpr Real PI = 3.141592653589793238462643383279L;
bool equals(Real x, Real y) { return fabs(x - y) < EPS; }
int sign(Real a) { return equals(a, 0) ? 0 : (a > 0 ? 1 : -1); }

template <class 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 P{x + r.x, y + r.y}; }
  P operator-(const P& r) const { return P{x - r.x, y - r.y}; }
  P operator-() const { return P{-x, -y}; }
  P operator*(R r) const { return P{x * r, y * r}; }
  P operator/(R r) const { return P{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)};
  }
  P rotate90() const { return {-y, x}; }

  R real() const { return x; }
  R imag() const { return y; }

  friend P operator*(R r, const P& p) { return p * r; }
  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>;

// relative position of c from a->b
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;        // counter-clockwise
  if (cross(x, y) < -EPS) return -1;       // clockwise
  if (dot(x, y) < -EPS) return +2;         // collinear in the order c-a-b
  if (norm(x) + EPS < norm(y)) return -2;  // collinear in the order a-b-c
  return 0;                                // collinear in the order a-c-b
}

/**
 * @brief 二次元幾何の基本要素
 * @docs docs/geometry/geometry-base.md
 */
#line 4 "geometry/polygon.hpp"

// 反時計周り
using Polygon = vector<Point>;

// 多角形の内部に点があるか
// OUT : 0, ON : 1, IN : 2
// 凸でなくてもよい, 自己交差はもたない
// crossing number algorithm : x 軸正方向に伸ばした半直線が奇数回交われば内部
int contains_polygon(const Polygon& P, const Point& p) {
  bool in = false;
  for (int i = 0; i < (int)P.size(); i++) {
    Point a = P[i] - p, b = P[(i + 1) % P.size()] - p;
    if (a.y > b.y) swap(a, b);
    if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1;
    if (sign(a.y) <= 0 && 0 < sign(b.y) && sign(cross(a, b)) < 0)
      in = !in;
  }
  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;
}

bool is_convex(const Polygon& P) {
  int n = P.size();
  if (n < 3) return true;
  for (int i = 0; i < n; ++i) {
    if (ccw(P[i], P[(i + 1) % n], P[(i + 2) % n]) == -1) return false;
  }
  return true;
}

// 頂点集合から凸包を生成
// 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;
  if constexpr (boundary) {
    bool collinear = true;
    for (int i = 2; i < n; ++i) {
      if (!equals(cross(ps[1] - ps[0], ps[i] - ps[0]), 0)) {
        collinear = false;
        break;
      }
    }
    if (collinear) return ps;
  }
  vector<Point> ch(2 * n);
  // 反時計周り
  const 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();
  if (N == 0) return 0;
  if (N == 1) return abs(C[0] - p) < EPS ? 1 : 0;
  if (N == 2) return ccw(C[0], C[1], p) == 0 ? 1 : 0;
  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& C) {
  int N = (int)C.size();
  assert(N > 0);
  if (N == 1) return {0, 0};
  if (N == 2) return {0, 1};
  int is = 0, js = 0;
  for (int i = 1; i < N; i++) {
    if (C[i].y > C[is].y) is = i;
    if (C[i].y < C[js].y) js = i;
  }
  Real maxdis = norm(C[is] - C[js]);

  int maxi, maxj, i, j;
  i = maxi = is;
  j = maxj = js;
  do {
    if (cross(C[(i + 1) % N] - C[i], C[(j + 1) % N] - C[j]) >= 0) {
      j = (j + 1) % N;
    } else {
      i = (i + 1) % N;
    }
    if (norm(C[i] - C[j]) > maxdis) {
      maxdis = norm(C[i] - C[j]);
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return minmax(maxi, maxj);
}

// min argmax_i dot(C[i], p)
int argmax_dot(const Polygon& C, const Point& p) {
  int n = C.size();
  assert(n > 0);
  assert(p.x != 0 || p.y != 0);
  if (n <= 2) {
    if (n == 1 || dot(C[0], p) >= dot(C[1], p)) return 0;
    return 1;
  }
  auto arg_half = [](const Point& a) {
    return a.y < 0 || (a.y == 0 && a.x < 0);
  };
  auto arg_less = [&](const Point& a, const Point& b) {
    int ha = arg_half(a), hb = arg_half(b);
    return ha != hb ? ha < hb : cross(a, b) > 0;
  };
  auto edge = [&](int i) { return C[(i + 1) % n] - C[i]; };
  int l = 0, r = n - 1;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge(r), edge(m)))
      l = m + 1;
    else
      r = m;
  }
  int s = l;
  Point q = p.rotate90();
  l = 0, r = n;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge((s + m) % n), q))
      l = m + 1;
    else
      r = m;
  }
  int i = (s + (l == n ? 0 : l)) % n;
  Point e = edge(i);
  if (!arg_less(e, q) && !arg_less(q, e)) {
    int low = l;
    l = low, r = n;
    while (l < r) {
      int m = (l + r) / 2;
      if (!arg_less(q, edge((s + m) % n)))
        l = m + 1;
      else
        r = m;
    }
    int count = l - low;
    return i + count >= n ? 0 : i;
  }
  return i;
}

struct ConvexPolygonCutResult {
  int first = 0, last = 0;
  optional<Point> front, back;

  bool empty() const { return first == last && !front && !back; }
  Polygon to_polygon(const Polygon& C) const {
    Polygon ret;
    auto push = [&](const Point& p) { if (ret.empty() || abs(ret.back() - p) >= EPS) ret.push_back(p); };
    if (front) push(*front);
    for (int i = first; i < last; ++i) push(C[i % C.size()]);
    if (back) push(*back);
    if (ret.size() >= 2 && abs(ret.front() - ret.back()) < EPS) ret.pop_back();
    return ret;
  }
};

ConvexPolygonCutResult convex_polygon_cut_info(const Polygon& C, const Point& a, const Point& b) {
  int n = C.size();
  assert(n > 0 && a != b);
  Point d = b - a, normal = d.rotate90();
  auto side = [&](int i) { return sign(cross(d, C[i % n] - a)); };
  int mx = argmax_dot(C, normal), mn = argmax_dot(C, -normal);
  if (side(mx) < 0) return {};
  if (side(mn) >= 0) return {0, n, nullopt, nullopt};

  int l = mn, r = mx;
  if (r < l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) >= 0 ? r : l) = m;
  }
  int first = r;
  l = mx, r = mn;
  if (r <= l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) < 0 ? r : l) = m;
  }
  int last = r;
  if (last < first) last += n;
  int count = last - first;
  first %= n;
  ConvexPolygonCutResult ret{first, first + count, nullopt, nullopt};
  auto cross_point = [&](const Point& p, const Point& q) {
    return p + (q - p) * (cross(d, a - p) / cross(d, q - p));
  };
  if (side(first) > 0) ret.front = cross_point(C[(first + n - 1) % n], C[first]);
  int back = (ret.last - 1) % n;
  if (side(back) > 0) ret.back = cross_point(C[back], C[(back + 1) % n]);
  return ret;
}

Polygon convex_polygon_cut(const Polygon& C, const Point& a, const Point& b) {
  return convex_polygon_cut_info(C, a, b).to_polygon(C);
}

/**
 * @brief 多角形
 * @docs docs/geometry/polygon.md
 */

多角形 (geometry/polygon.hpp)

資料

Depends on

Required by

Verified with

Code

多角形
(geometry/polygon.hpp)