This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub dnx04/algo
#include "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
#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