ACM_Notebook_new
This documentation is automatically generated by online-judge-tools/verification-helper
Geometry/polygon.h
- View this file on GitHub
- Last update: 2023-02-02 12:34:59+08:00
- 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
Verified with
- Geometry/tests/aizu_cgl_3_a_polygon_area.test.cpp
- Geometry/tests/aizu_cgl_3_b_polygon_is_convex.test.cpp
- Geometry/tests/aizu_cgl_3_c_polygon_in_polygon.test.cpp
- Geometry/tests/aizu_cgl_4_a_polygon_convex_hull.test.cpp
- Geometry/tests/aizu_cgl_4_b_polygon_convex_diameter.test.cpp
- Geometry/tests/aizu_cgl_4_c_polygon_convex_cut.test.cpp
- Geometry/tests/polygon_in_convex.test.cpp
- Geometry/tests/z_polygon_area.test.cpp
- Geometry/tests/z_polygon_convexhull.test.cpp
- Geometry/tests/z_polygon_is_convex.test.cpp
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.h"
// 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