geometry/polygon.hpp

View this file on GitHub
Last update: 2023-08-13 15:15:35+07:00
Include: #include "geometry/polygon.hpp"
Link: https://cses.fi/problemset/task/2195
Link: https://open.kattis.com/problems/convexhull
Link: https://open.kattis.com/problems/cuttingpolygon
Link: https://open.kattis.com/problems/pointinpolygon

Code

// Polygon: convex hull, in polygon, convex diameter ..
// Functions with param vector<P<T>> works with both integers and floating
// points Functions with param Polygon works with floating points only.

typedef vector<Point> Polygon;

// Convex Hull:
// If minimum point --> #define REMOVE_REDUNDANT
// Known issues:
// - Max. point does not work when some points are the same.
// Tested:
// - (min points) https://open.kattis.com/problems/convexhull
// - (max points) https://cses.fi/problemset/task/2195

// #define REMOVE_REDUNDANT
template <typename T>
T area2(P<T> a, P<T> b, P<T> c) {
  return a % b + b % c + c % a;
}

template <typename T>
void ConvexHull(vector<P<T>> &pts) {
  sort(pts.begin(), pts.end());
  pts.erase(unique(pts.begin(), pts.end()), pts.end());
  vector<P<T>> up, dn;
  for (int i = 0; i < (int)pts.size(); i++) {
#ifdef REMOVE_REDUNDANT
    while (up.size() > 1 && area2(up[up.size() - 2], up.back(), pts[i]) >= 0)
      up.pop_back();
    while (dn.size() > 1 && area2(dn[dn.size() - 2], dn.back(), pts[i]) <= 0)
      dn.pop_back();
#else
    while (up.size() > 1 && area2(up[up.size() - 2], up.back(), pts[i]) > 0)
      up.pop_back();
    while (dn.size() > 1 && area2(dn[dn.size() - 2], dn.back(), pts[i]) < 0)
      dn.pop_back();
#endif
    up.push_back(pts[i]);
    dn.push_back(pts[i]);
  }
  pts = dn;
  for (int i = (int)up.size() - 2; i >= 1; i--) pts.push_back(up[i]);

#ifdef REMOVE_REDUNDANT
  if (pts.size() <= 2) return;
  dn.clear();
  dn.push_back(pts[0]);
  dn.push_back(pts[1]);
  for (int i = 2; i < (int)pts.size(); i++) {
    if (onSegment(dn[dn.size() - 2], pts[i], dn.back())) dn.pop_back();
    dn.push_back(pts[i]);
  }
  if (dn.size() >= 3 && onSegment(dn.back(), dn[1], dn[0])) {
    dn[0] = dn.back();
    dn.pop_back();
  }
  pts = dn;
#endif
}

// Area, perimeter, centroid
template <typename T>
T signed_area2(vector<P<T>> p) {
  T area(0);
  for (int i = 0; i < (int)p.size(); i++) {
    area += p[i] % p[(i + 1) % p.size()];
  }
  return area;
}
double area(const Polygon &p) { return std::abs(signed_area2(p) / 2.0); }
Point centroid(Polygon p) {
  Point c(0, 0);
  double scale = 3.0 * signed_area2(p);
  for (int i = 0; i < (int)p.size(); i++) {
    int j = (i + 1) % p.size();
    c = c + (p[i] + p[j]) * (p[i].x * p[j].y - p[j].x * p[i].y);
  }
  return c / scale;
}
double perimeter(Polygon p) {
  double res = 0;
  for (int i = 0; i < (int)p.size(); ++i) {
    int j = (i + 1) % p.size();
    res += (p[i] - p[j]).len();
  }
  return res;
}

// Is convex: checks if polygon is convex.
// Return true for degen points (all vertices are collinear)
template <typename T>
bool is_convex(const vector<P<T>> &P) {
  int sz = (int)P.size();
  int min_ccw = 2, max_ccw = -2;
  for (int i = 0; i < sz; i++) {
    int c = ccw(P[i], P[(i + 1) % sz], P[(i + 2) % sz]);
    min_ccw = min(min_ccw, c);
    max_ccw = max(max_ccw, c);
  }
  return min_ccw * max_ccw >= 0;
}

// Inside polygon: O(N). Works with any polygon
// Tested:
// - https://open.kattis.com/problems/pointinpolygon
// - https://open.kattis.com/problems/cuttingpolygon
enum PolygonLocation { OUT, ON, IN };
PolygonLocation in_polygon(const Polygon &p, Point q) {
  if ((int)p.size() == 0) return PolygonLocation::OUT;

  // Check if point is on edge.
  int n = p.size();
  REP(i, n) {
    int j = (i + 1) % n;
    Point u = p[i], v = p[j];

    if (u > v) swap(u, v);

    if (ccw(u, v, q) == 0 && u <= q && q <= v) return PolygonLocation::ON;
  }

  // Check if point is strictly inside.
  int c = 0;
  for (int i = 0; i < n; i++) {
    int j = (i + 1) % n;
    if (((p[i].y <= q.y && q.y < p[j].y) || (p[j].y <= q.y && q.y < p[i].y)) &&
        q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) /
                           (double)(p[j].y - p[i].y)) {
      c = !c;
    }
  }
  return c ? PolygonLocation::IN : PolygonLocation::OUT;
}

// Check point in convex polygon, O(logN)
#define Det(a, b, c)                           \
  ((double)(b.x - a.x) * (double)(c.y - a.y) - \
   (double)(b.y - a.y) * (c.x - a.x))
PolygonLocation in_convex(vector<Point> &l, Point p) {
  if (l.empty()) return PolygonLocation::OUT;
  if (l.size() <= 2) {
    return onSegment(l[0], l[1 % l.size()], p) ? PolygonLocation::ON
                                               : PolygonLocation::OUT;
  }

  int a = 1, b = l.size() - 1, c;
  if (Det(l[0], l[a], l[b]) > 0) swap(a, b);

  if (onSegment(l[0], l[a], p)) return ON;
  if (onSegment(l[0], l[b], p)) return ON;

  if (Det(l[0], l[a], p) > 0 || Det(l[0], l[b], p) < 0) return OUT;
  while (abs(a - b) > 1) {
    c = (a + b) / 2;
    if (Det(l[0], l[c], p) > 0)
      b = c;
    else
      a = c;
  }
  int t = cmp(Det(l[a], l[b], p), 0);
  return (t == 0) ? ON : (t < 0) ? IN : OUT;
}

// Cut a polygon with a line. Returns half on left-hand-side.
// To return the other half, reverse the direction of Line l (by negating l.a,
// l.b) The line must be formed using 2 points
Polygon polygon_cut(const Polygon &P, Line l) {
  Polygon Q;
  for (int i = 0; i < (int)P.size(); ++i) {
    Point A = P[i], B = (i == ((int)P.size()) - 1) ? P[0] : P[i + 1];
    if (ccw(l.A, l.B, A) != -1) Q.push_back(A);
    if (ccw(l.A, l.B, A) * ccw(l.A, l.B, B) < 0) {
      Point p;
      areIntersect(Line(A, B), l, p);
      Q.push_back(p);
    }
  }
  return Q;
}

// Find intersection of 2 convex polygons
// Helper method
bool intersect_1pt(Point a, Point b, Point c, Point d, Point &r) {
  double D = (b - a) % (d - c);
  if (cmp(D, 0) == 0) return false;
  double t = ((c - a) % (d - c)) / D;
  double s = -((a - c) % (b - a)) / D;
  r = a + (b - a) * t;
  return cmp(t, 0) >= 0 && cmp(t, 1) <= 0 && cmp(s, 0) >= 0 && cmp(s, 1) <= 0;
}
Polygon convex_intersect(Polygon P, Polygon Q) {
  const int n = P.size(), m = Q.size();
  int a = 0, b = 0, aa = 0, ba = 0;
  enum { Pin, Qin, Unknown } in = Unknown;
  Polygon R;
  do {
    int a1 = (a + n - 1) % n, b1 = (b + m - 1) % m;
    double C = (P[a] - P[a1]) % (Q[b] - Q[b1]);
    double A = (P[a1] - Q[b]) % (P[a] - Q[b]);
    double B = (Q[b1] - P[a]) % (Q[b] - P[a]);
    Point r;
    if (intersect_1pt(P[a1], P[a], Q[b1], Q[b], r)) {
      if (in == Unknown) aa = ba = 0;
      R.push_back(r);
      in = B > 0 ? Pin : A > 0 ? Qin : in;
    }
    if (C == 0 && B == 0 && A == 0) {
      if (in == Pin) {
        b = (b + 1) % m;
        ++ba;
      } else {
        a = (a + 1) % m;
        ++aa;
      }
    } else if (C >= 0) {
      if (A > 0) {
        if (in == Pin) R.push_back(P[a]);
        a = (a + 1) % n;
        ++aa;
      } else {
        if (in == Qin) R.push_back(Q[b]);
        b = (b + 1) % m;
        ++ba;
      }
    } else {
      if (B > 0) {
        if (in == Qin) R.push_back(Q[b]);
        b = (b + 1) % m;
        ++ba;
      } else {
        if (in == Pin) R.push_back(P[a]);
        a = (a + 1) % n;
        ++aa;
      }
    }
  } while ((aa < n || ba < m) && aa < 2 * n && ba < 2 * m);
  if (in == Unknown) {
    if (in_convex(Q, P[0])) return P;
    if (in_convex(P, Q[0])) return Q;
  }
  return R;
}

// Find the diameter of polygon.
// Returns:
// <diameter, <ids of 2 points>>
pair<double, pair<int, int>> convex_diameter(const Polygon &p) {
  const int n = p.size();
  int is = 0, js = 0;
  for (int i = 1; i < n; ++i) {
    if (p[i].y > p[is].y) is = i;
    if (p[i].y < p[js].y) js = i;
  }
  double maxd = (p[is] - p[js]).len();
  int i, maxi, j, maxj;
  i = maxi = is;
  j = maxj = js;
  do {
    int ii = (i + 1) % n, jj = (j + 1) % n;
    if ((p[ii] - p[i]) % (p[jj] - p[j]) >= 0)
      j = jj;
    else
      i = ii;
    if ((p[i] - p[j]).len() > maxd) {
      maxd = (p[i] - p[j]).len();
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return {maxd, std::minmax(maxi, maxj)}; /* farthest pair is (maxi, maxj). */
}

// Pick theorem
// Given non-intersecting polygon.
// S = area
// I = number of integer points strictly Inside
// B = number of points on sides of polygon
// S = I + B/2 - 1

#line 1 "geometry/polygon.hpp"
// Polygon: convex hull, in polygon, convex diameter ..
// Functions with param vector<P<T>> works with both integers and floating
// points Functions with param Polygon works with floating points only.

typedef vector<Point> Polygon;

// Convex Hull:
// If minimum point --> #define REMOVE_REDUNDANT
// Known issues:
// - Max. point does not work when some points are the same.
// Tested:
// - (min points) https://open.kattis.com/problems/convexhull
// - (max points) https://cses.fi/problemset/task/2195

// #define REMOVE_REDUNDANT
template <typename T>
T area2(P<T> a, P<T> b, P<T> c) {
  return a % b + b % c + c % a;
}

template <typename T>
void ConvexHull(vector<P<T>> &pts) {
  sort(pts.begin(), pts.end());
  pts.erase(unique(pts.begin(), pts.end()), pts.end());
  vector<P<T>> up, dn;
  for (int i = 0; i < (int)pts.size(); i++) {
#ifdef REMOVE_REDUNDANT
    while (up.size() > 1 && area2(up[up.size() - 2], up.back(), pts[i]) >= 0)
      up.pop_back();
    while (dn.size() > 1 && area2(dn[dn.size() - 2], dn.back(), pts[i]) <= 0)
      dn.pop_back();
#else
    while (up.size() > 1 && area2(up[up.size() - 2], up.back(), pts[i]) > 0)
      up.pop_back();
    while (dn.size() > 1 && area2(dn[dn.size() - 2], dn.back(), pts[i]) < 0)
      dn.pop_back();
#endif
    up.push_back(pts[i]);
    dn.push_back(pts[i]);
  }
  pts = dn;
  for (int i = (int)up.size() - 2; i >= 1; i--) pts.push_back(up[i]);

#ifdef REMOVE_REDUNDANT
  if (pts.size() <= 2) return;
  dn.clear();
  dn.push_back(pts[0]);
  dn.push_back(pts[1]);
  for (int i = 2; i < (int)pts.size(); i++) {
    if (onSegment(dn[dn.size() - 2], pts[i], dn.back())) dn.pop_back();
    dn.push_back(pts[i]);
  }
  if (dn.size() >= 3 && onSegment(dn.back(), dn[1], dn[0])) {
    dn[0] = dn.back();
    dn.pop_back();
  }
  pts = dn;
#endif
}

// Area, perimeter, centroid
template <typename T>
T signed_area2(vector<P<T>> p) {
  T area(0);
  for (int i = 0; i < (int)p.size(); i++) {
    area += p[i] % p[(i + 1) % p.size()];
  }
  return area;
}
double area(const Polygon &p) { return std::abs(signed_area2(p) / 2.0); }
Point centroid(Polygon p) {
  Point c(0, 0);
  double scale = 3.0 * signed_area2(p);
  for (int i = 0; i < (int)p.size(); i++) {
    int j = (i + 1) % p.size();
    c = c + (p[i] + p[j]) * (p[i].x * p[j].y - p[j].x * p[i].y);
  }
  return c / scale;
}
double perimeter(Polygon p) {
  double res = 0;
  for (int i = 0; i < (int)p.size(); ++i) {
    int j = (i + 1) % p.size();
    res += (p[i] - p[j]).len();
  }
  return res;
}

// Is convex: checks if polygon is convex.
// Return true for degen points (all vertices are collinear)
template <typename T>
bool is_convex(const vector<P<T>> &P) {
  int sz = (int)P.size();
  int min_ccw = 2, max_ccw = -2;
  for (int i = 0; i < sz; i++) {
    int c = ccw(P[i], P[(i + 1) % sz], P[(i + 2) % sz]);
    min_ccw = min(min_ccw, c);
    max_ccw = max(max_ccw, c);
  }
  return min_ccw * max_ccw >= 0;
}

// Inside polygon: O(N). Works with any polygon
// Tested:
// - https://open.kattis.com/problems/pointinpolygon
// - https://open.kattis.com/problems/cuttingpolygon
enum PolygonLocation { OUT, ON, IN };
PolygonLocation in_polygon(const Polygon &p, Point q) {
  if ((int)p.size() == 0) return PolygonLocation::OUT;

  // Check if point is on edge.
  int n = p.size();
  REP(i, n) {
    int j = (i + 1) % n;
    Point u = p[i], v = p[j];

    if (u > v) swap(u, v);

    if (ccw(u, v, q) == 0 && u <= q && q <= v) return PolygonLocation::ON;
  }

  // Check if point is strictly inside.
  int c = 0;
  for (int i = 0; i < n; i++) {
    int j = (i + 1) % n;
    if (((p[i].y <= q.y && q.y < p[j].y) || (p[j].y <= q.y && q.y < p[i].y)) &&
        q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) /
                           (double)(p[j].y - p[i].y)) {
      c = !c;
    }
  }
  return c ? PolygonLocation::IN : PolygonLocation::OUT;
}

// Check point in convex polygon, O(logN)
#define Det(a, b, c)                           \
  ((double)(b.x - a.x) * (double)(c.y - a.y) - \
   (double)(b.y - a.y) * (c.x - a.x))
PolygonLocation in_convex(vector<Point> &l, Point p) {
  if (l.empty()) return PolygonLocation::OUT;
  if (l.size() <= 2) {
    return onSegment(l[0], l[1 % l.size()], p) ? PolygonLocation::ON
                                               : PolygonLocation::OUT;
  }

  int a = 1, b = l.size() - 1, c;
  if (Det(l[0], l[a], l[b]) > 0) swap(a, b);

  if (onSegment(l[0], l[a], p)) return ON;
  if (onSegment(l[0], l[b], p)) return ON;

  if (Det(l[0], l[a], p) > 0 || Det(l[0], l[b], p) < 0) return OUT;
  while (abs(a - b) > 1) {
    c = (a + b) / 2;
    if (Det(l[0], l[c], p) > 0)
      b = c;
    else
      a = c;
  }
  int t = cmp(Det(l[a], l[b], p), 0);
  return (t == 0) ? ON : (t < 0) ? IN : OUT;
}

// Cut a polygon with a line. Returns half on left-hand-side.
// To return the other half, reverse the direction of Line l (by negating l.a,
// l.b) The line must be formed using 2 points
Polygon polygon_cut(const Polygon &P, Line l) {
  Polygon Q;
  for (int i = 0; i < (int)P.size(); ++i) {
    Point A = P[i], B = (i == ((int)P.size()) - 1) ? P[0] : P[i + 1];
    if (ccw(l.A, l.B, A) != -1) Q.push_back(A);
    if (ccw(l.A, l.B, A) * ccw(l.A, l.B, B) < 0) {
      Point p;
      areIntersect(Line(A, B), l, p);
      Q.push_back(p);
    }
  }
  return Q;
}

// Find intersection of 2 convex polygons
// Helper method
bool intersect_1pt(Point a, Point b, Point c, Point d, Point &r) {
  double D = (b - a) % (d - c);
  if (cmp(D, 0) == 0) return false;
  double t = ((c - a) % (d - c)) / D;
  double s = -((a - c) % (b - a)) / D;
  r = a + (b - a) * t;
  return cmp(t, 0) >= 0 && cmp(t, 1) <= 0 && cmp(s, 0) >= 0 && cmp(s, 1) <= 0;
}
Polygon convex_intersect(Polygon P, Polygon Q) {
  const int n = P.size(), m = Q.size();
  int a = 0, b = 0, aa = 0, ba = 0;
  enum { Pin, Qin, Unknown } in = Unknown;
  Polygon R;
  do {
    int a1 = (a + n - 1) % n, b1 = (b + m - 1) % m;
    double C = (P[a] - P[a1]) % (Q[b] - Q[b1]);
    double A = (P[a1] - Q[b]) % (P[a] - Q[b]);
    double B = (Q[b1] - P[a]) % (Q[b] - P[a]);
    Point r;
    if (intersect_1pt(P[a1], P[a], Q[b1], Q[b], r)) {
      if (in == Unknown) aa = ba = 0;
      R.push_back(r);
      in = B > 0 ? Pin : A > 0 ? Qin : in;
    }
    if (C == 0 && B == 0 && A == 0) {
      if (in == Pin) {
        b = (b + 1) % m;
        ++ba;
      } else {
        a = (a + 1) % m;
        ++aa;
      }
    } else if (C >= 0) {
      if (A > 0) {
        if (in == Pin) R.push_back(P[a]);
        a = (a + 1) % n;
        ++aa;
      } else {
        if (in == Qin) R.push_back(Q[b]);
        b = (b + 1) % m;
        ++ba;
      }
    } else {
      if (B > 0) {
        if (in == Qin) R.push_back(Q[b]);
        b = (b + 1) % m;
        ++ba;
      } else {
        if (in == Pin) R.push_back(P[a]);
        a = (a + 1) % n;
        ++aa;
      }
    }
  } while ((aa < n || ba < m) && aa < 2 * n && ba < 2 * m);
  if (in == Unknown) {
    if (in_convex(Q, P[0])) return P;
    if (in_convex(P, Q[0])) return Q;
  }
  return R;
}

// Find the diameter of polygon.
// Returns:
// <diameter, <ids of 2 points>>
pair<double, pair<int, int>> convex_diameter(const Polygon &p) {
  const int n = p.size();
  int is = 0, js = 0;
  for (int i = 1; i < n; ++i) {
    if (p[i].y > p[is].y) is = i;
    if (p[i].y < p[js].y) js = i;
  }
  double maxd = (p[is] - p[js]).len();
  int i, maxi, j, maxj;
  i = maxi = is;
  j = maxj = js;
  do {
    int ii = (i + 1) % n, jj = (j + 1) % n;
    if ((p[ii] - p[i]) % (p[jj] - p[j]) >= 0)
      j = jj;
    else
      i = ii;
    if ((p[i] - p[j]).len() > maxd) {
      maxd = (p[i] - p[j]).len();
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return {maxd, std::minmax(maxi, maxj)}; /* farthest pair is (maxi, maxj). */
}

// Pick theorem
// Given non-intersecting polygon.
// S = area
// I = number of integer points strictly Inside
// B = number of points on sides of polygon
// S = I + B/2 - 1