#include <bits/stdc++.h> using namespace std; //Starting from here code is imported from: https://github.com/bicsi/kactl/tree/master/content/geometry #define DOUBLE long double const double kInf = 2.0e5; struct Point { DOUBLE x, y; Point operator+(const Point& oth) const { return {x + oth.x, y + oth.y}; } Point operator-(const Point& oth) const { return {x - oth.x, y - oth.y}; } Point operator*(const DOUBLE d) const { return {x * d, y * d}; } Point operator/(const DOUBLE d) const { return {x / d, y / d}; } }; int sgn(DOUBLE d) { return max(d, -d) < 1e-5 ? 0 : d < 0 ? -1 : 1; } DOUBLE cross(Point a, Point b) { return a.x * b.y - a.y * b.x; } DOUBLE dot(Point a, Point b) { return a.x * b.x + a.y * b.y; } DOUBLE det(Point a, Point b, Point c) { return cross(b - a, c - a); } Point LineIntersection(Point a, Point b, Point p, Point q) { DOUBLE c1 = det(a, b, p), c2 = det(a, b, q); assert(sgn(c2 - c1)); return (q * c1 - p * c2) / (c1 - c2); } struct Angle { long long x, y; Angle(long long x, long long y) : x(x), y(y) {} Angle operator-(Angle a) const { return {x-a.x, y-a.y}; } int quad() const { if (y < 0) return (x >= 0) + 2; if (y > 0) return (x <= 0); return (x <= 0) * 2; } Point p() const { return {1. * x, 1. * y}; } Angle t90() const { return {-y, x}; } }; bool operator<(Angle a, Angle b) { return make_pair(a.quad(), a.y * b.x) < make_pair(b.quad(), a.x * b.y); } struct HalfplaneSet : multimap<Angle, Point> { using Iter = multimap<Angle, Point>::iterator; HalfplaneSet() { insert({{+1, 0}, {-kInf, -kInf}}); insert({{0, +1}, {+kInf, -kInf}}); insert({{-1, 0}, {+kInf, +kInf}}); insert({{0, -1}, {-kInf, +kInf}}); } Iter get_next(Iter it) { return (next(it) == end() ? begin() : next(it)); } Iter get_prev(Iter it) { return (it == begin() ? prev(end()) : prev(it)); } Iter fix(Iter it) { return it == end() ? begin() : it; } void Cut(Angle a, Angle b) { if (empty()) return; int old_size = size(); auto eval = [&](Iter it) { return sgn(det(a.p(), b.p(), it->second)); }; auto intersect = [&](Iter it) { return LineIntersection(a.p(), b.p(), it->second, it->first.p() + it->second); }; auto it = fix(lower_bound(b - a)); if (eval(it) >= 0) return; while (size() && eval(get_prev(it)) < 0) fix(erase(get_prev(it))); while (size() && eval(get_next(it)) < 0) it = fix(erase(it)); if (empty()) return; if (eval(get_next(it)) > 0) it->second = intersect(it); else it = fix(erase(it)); if (old_size <= 2) return; it = get_prev(it); auto p = intersect(it); insert(it, {b - a, p}); if (eval(it) == 0) erase(it); } DOUBLE Maximize(Angle c) { assert(!empty()); auto it = fix(lower_bound(c.t90())); return dot(it->second, c.p()); } void Check(Angle c) { if (empty()) return; DOUBLE ans = -2e18; for (auto it = begin(); it != end(); ++it) ans = max(ans, dot(it->second, c.p())); assert(sgn(Maximize(c) - ans) == 0); } void output() { cerr << "Polygon\n"; for (auto x : *this) { cerr << x.second.x << " " << x.second.y << '\n'; cerr << "dir: " << x.first.x << " " << x.first.y << '\n'; } cerr << "...\n"; } DOUBLE Area() { if (size() <= 2) return 0; DOUBLE ret = 0; for (auto it = begin(); it != end(); ++it) ret += cross(it->second, get_next(it)->second); return ret; } }; //End of imported code int main() { int n; cin >> n; vector<vector<Angle>> p(n, vector<Angle>()); for(int i = 0; i < n; i += 1) for(int j = 0; j < 2; j += 1) { int x, y; cin >> x >> y; p[i].push_back(Angle(x, y)); } HalfplaneSet hs; for(int i = 0; i < n; i += 1) for(int j = 0; j < n; j += 1) { if(i == j) continue; for(int ii = 0; ii < 2; ii += 1) for(int jj = 0; jj < 2; jj += 1) { bool can = true; for(int k = 0; k < n; k += 1) { if(det(p[i][ii].p(), p[j][jj].p(), p[k][0].p()) > 0 and det(p[i][ii].p(), p[j][jj].p(), p[k][1].p()) > 0) { can = false; } } if(can) { hs.Cut(p[j][jj], p[i][ii]); } } } cout.precision(20); cout << hs.Area() / 2 << "\n"; }
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 | #include <bits/stdc++.h> using namespace std; //Starting from here code is imported from: https://github.com/bicsi/kactl/tree/master/content/geometry #define DOUBLE long double const double kInf = 2.0e5; struct Point { DOUBLE x, y; Point operator+(const Point& oth) const { return {x + oth.x, y + oth.y}; } Point operator-(const Point& oth) const { return {x - oth.x, y - oth.y}; } Point operator*(const DOUBLE d) const { return {x * d, y * d}; } Point operator/(const DOUBLE d) const { return {x / d, y / d}; } }; int sgn(DOUBLE d) { return max(d, -d) < 1e-5 ? 0 : d < 0 ? -1 : 1; } DOUBLE cross(Point a, Point b) { return a.x * b.y - a.y * b.x; } DOUBLE dot(Point a, Point b) { return a.x * b.x + a.y * b.y; } DOUBLE det(Point a, Point b, Point c) { return cross(b - a, c - a); } Point LineIntersection(Point a, Point b, Point p, Point q) { DOUBLE c1 = det(a, b, p), c2 = det(a, b, q); assert(sgn(c2 - c1)); return (q * c1 - p * c2) / (c1 - c2); } struct Angle { long long x, y; Angle(long long x, long long y) : x(x), y(y) {} Angle operator-(Angle a) const { return {x-a.x, y-a.y}; } int quad() const { if (y < 0) return (x >= 0) + 2; if (y > 0) return (x <= 0); return (x <= 0) * 2; } Point p() const { return {1. * x, 1. * y}; } Angle t90() const { return {-y, x}; } }; bool operator<(Angle a, Angle b) { return make_pair(a.quad(), a.y * b.x) < make_pair(b.quad(), a.x * b.y); } struct HalfplaneSet : multimap<Angle, Point> { using Iter = multimap<Angle, Point>::iterator; HalfplaneSet() { insert({{+1, 0}, {-kInf, -kInf}}); insert({{0, +1}, {+kInf, -kInf}}); insert({{-1, 0}, {+kInf, +kInf}}); insert({{0, -1}, {-kInf, +kInf}}); } Iter get_next(Iter it) { return (next(it) == end() ? begin() : next(it)); } Iter get_prev(Iter it) { return (it == begin() ? prev(end()) : prev(it)); } Iter fix(Iter it) { return it == end() ? begin() : it; } void Cut(Angle a, Angle b) { if (empty()) return; int old_size = size(); auto eval = [&](Iter it) { return sgn(det(a.p(), b.p(), it->second)); }; auto intersect = [&](Iter it) { return LineIntersection(a.p(), b.p(), it->second, it->first.p() + it->second); }; auto it = fix(lower_bound(b - a)); if (eval(it) >= 0) return; while (size() && eval(get_prev(it)) < 0) fix(erase(get_prev(it))); while (size() && eval(get_next(it)) < 0) it = fix(erase(it)); if (empty()) return; if (eval(get_next(it)) > 0) it->second = intersect(it); else it = fix(erase(it)); if (old_size <= 2) return; it = get_prev(it); auto p = intersect(it); insert(it, {b - a, p}); if (eval(it) == 0) erase(it); } DOUBLE Maximize(Angle c) { assert(!empty()); auto it = fix(lower_bound(c.t90())); return dot(it->second, c.p()); } void Check(Angle c) { if (empty()) return; DOUBLE ans = -2e18; for (auto it = begin(); it != end(); ++it) ans = max(ans, dot(it->second, c.p())); assert(sgn(Maximize(c) - ans) == 0); } void output() { cerr << "Polygon\n"; for (auto x : *this) { cerr << x.second.x << " " << x.second.y << '\n'; cerr << "dir: " << x.first.x << " " << x.first.y << '\n'; } cerr << "...\n"; } DOUBLE Area() { if (size() <= 2) return 0; DOUBLE ret = 0; for (auto it = begin(); it != end(); ++it) ret += cross(it->second, get_next(it)->second); return ret; } }; //End of imported code int main() { int n; cin >> n; vector<vector<Angle>> p(n, vector<Angle>()); for(int i = 0; i < n; i += 1) for(int j = 0; j < 2; j += 1) { int x, y; cin >> x >> y; p[i].push_back(Angle(x, y)); } HalfplaneSet hs; for(int i = 0; i < n; i += 1) for(int j = 0; j < n; j += 1) { if(i == j) continue; for(int ii = 0; ii < 2; ii += 1) for(int jj = 0; jj < 2; jj += 1) { bool can = true; for(int k = 0; k < n; k += 1) { if(det(p[i][ii].p(), p[j][jj].p(), p[k][0].p()) > 0 and det(p[i][ii].p(), p[j][jj].p(), p[k][1].p()) > 0) { can = false; } } if(can) { hs.Cut(p[j][jj], p[i][ii]); } } } cout.precision(20); cout << hs.Area() / 2 << "\n"; } |