半平面交學習筆記
阿新 • • 發佈:2018-11-20
半平面交 - (S & I algorithm)
參考論文 演算法合集之《半平面交的新演算法及其實用價值》
問題簡介:
給出多個如 \(ax + by + c \ge 0\) 的限制( 接下來都以 \(ax+by+c \ge 0\) 為例) , 求解 \((x,y)\) 的集合
可以轉化為多個直線在平面上圍成的凸包
step1
將所有直線按角度排序,角度相同的保留下需要的一個(如圖)
step2
用一個雙端佇列儲存當前半平面交,每次通過判斷隊首與隊尾第一個交點是否滿足當前直線來更新
step3
先用隊尾判定隊首交點是否合法,再用隊首判斷隊尾交點是否合法
step4
現在佇列中的相鄰半平面的交點即為凸包的節點, 如果剩餘半平面數量小於3則無解
Code(POJ1279)
\\ 滿足Plane a的點為a.s->a.t的逆時針方向 #include <algorithm> #include <cmath> #include <cstdio> #include <iomanip> #include <iostream> #include <queue> #define sz q.size() using namespace std; inline char nc() { // return getchar(); static char buf[100000], * l = buf, * r = buf; if(l == r) r = (l = buf) + fread(buf, 1, 100000, stdin); if(l == r) return EOF; return *l++; } template<class T> void readin(T & x) { x = 0; int f = 1, ch = nc(); while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=nc();} while(ch>='0'&&ch<='9'){x=x*10-'0'+ch;ch=nc();} x *= f; } const double eps = 1e-9; int sign(double a) { if(fabs(a) < eps) return 0; return a < 0 ? -1 : 1; } struct Point { double x, y; Point() {} Point(double x, double y) : x(x), y(y) {} }; inline Point operator + (const Point & a, const Point & b) { return Point(a.x + b.x, a.y + b.y); } inline Point operator - (const Point & a, const Point & b) { return Point(a.x - b.x, a.y - b.y); } inline Point operator * (const Point & a, const double & b) { return Point(a.x * b, a.y * b); } inline Point operator / (const Point & a, const double & b) { return Point(a.x / b, a.y / b); } inline double cross(Point a, Point b) { return a.x * b.y - a.y * b.x; } inline double angle(Point a) { return atan2(a.y, a.x); } inline Point intersection(Point p0, Point p1, Point q0, Point q1) { return p0 + (p1 - p0) * cross(q1 - q0, p0 - q0) / cross(p1 - p0, q1 - q0); } inline double area(Point * p, int n) { double res = 0; p[n + 1] = p[1]; for(int i = 1; i <= n; i++) { res += cross(p[i], p[i + 1]); } return fabs(res / 2.0); } struct Plane { Point s, t; double ang; Plane() {} Plane(Point s, Point t) : s(s), t(t) { ang = angle(t - s); } }; int m; Point s[1505]; int comp(const Plane & a, const Plane & b) { double r = a.ang - b.ang; if(sign(r) != 0) return sign(r) == -1; return sign(cross(a.t - a.s, b.t - a.s)) == -1; } int judge(Plane p, Plane a, Plane b) { Point q = intersection(a.s, a.t, b.s, b.t); return sign(cross(p.t - p.s, q - p.s)) == -1; } bool SI(Plane * l, int n) { sort(l + 1, l + n + 1, comp); int t = 0; for(int i = 1; i <= n; i++) { if(i == 1 || sign(l[i].ang - l[i - 1].ang) != 0) { l[++t] = l[i]; } } n = t; if(n < 3) return false; deque<Plane> q; q.push_front(l[1]); q.push_front(l[2]); for(int i = 3; i <= n; i++) { while(sz > 1 && judge(l[i], q[0], q[1])) q.pop_front(); while(sz > 1 && judge(l[i], q[sz - 1], q[sz - 2])) q.pop_back(); q.push_front(l[i]); } while(sz > 1 && judge(q[sz - 1], q[0], q[1])) q.pop_front(); while(sz > 1 && judge(q[0], q[sz - 1], q[sz - 2])) q.pop_back(); if(q.size() < 3) return false; m = q.size(); q.push_back(q[0]); for(int i = 1; i <= m; i++) { s[i] = intersection(q[i - 1].s, q[i - 1].t, q[i].s, q[i].t); } return true; } int T, n; Plane l[1505]; Point p[1505]; double solve() { p[n + 1] = p[1]; for(int i = 1; i <= n; i++) { l[i] = Plane(p[i + 1], p[i]); } if(!SI(l, n)) return 0; return area(s, m); } int main() { // freopen("input.txt", "r", stdin); cout << fixed; readin(T); for(int kase = 1; kase <= T; kase++) { readin(n); for(int i = 1; i <= n; i++) { readin(p[i].x); readin(p[i].y); } cout << setprecision(2) << solve() << endl; } return 0; }