模板-計算幾何
阿新 • • 發佈:2019-02-13
struct Point { double x, y; Point(double x = 0, double y = 0): x(x), y(y) {} }; typedef Point Vector; Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); } Vector operator - (Vector A, Vector B) { return Vector(A.x-B.x, A.y-B.y); } Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); } Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); } bool operator < (const Vector& a, const Vector& b) { return a.x < b.x || (a.x == b.x && a.y < b.y); } const double eps = 1e-10; int dcmp(double x) { if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; } bool operator == (const Vector& a, const Vector& b) { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; } double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; } double Length(Vector A) { return sqrt(Dot(A, A)); } double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); } double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; } double Area2(Point A, Point B, Point C) { return Cross(B-A, C-A); } // 旋轉公式 用 (length*cos(theta), length*sin(theta)) 表示向量即可推出 Vector Rotate(Vector A, double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad)); } // 法向量 Vector Normal(Vector A) { double L = Length(A); return Vector(-A.y/L, A.x/L); } // P(x1, y1), v(x2, y2), Q(x3, y3), w(x4, y4), P-Q=u(x5, y5). // P + v*t1 = Q + w*t2 // 用虛數表示如下 // x1+y1*i + t1*x2+t1*y2*i = x3+y3*i + t2*x4+t2*y4*i // x1 + t1*x2 + (y1+t1*y2)*i = x3 + t2*x4 + (y3+t2*y4)*i // => /x1+t1*x2 = x3+t2*x4 -> x5 + t1*x2 = t2*x4 -> x5*y4 + t1*x2*y4 = t2*x4*y4 (1) // \y1+t1*y2 = y3+t2*y4 -> y5 + t1*y2 = t2*y4 -> y5*x4 + t1*y2*x4 = t2*y4*x4 (2) // (1)-(2) => (x5*y4-y5*x4) + t1(x2*y4-y2*x4) = 0 // => t1 = (x4*y5-x5*y4) / (x2*y4-y2*x4) // = Cross(w, u) / Cross(v, w) Vector GetLineIntersection(Point P, Vector v, Point Q, Vector w) { Vector u = P-Q; double t = Cross(w, u) / Cross(v, w); return P+v*t; } double DistanceToLine(Point P, Point A, Point B) { Vector v1 = B-A, v2 = P-A; return fabs(Cross(v1, v2)) / Length(v1); } double DistanceToSegment(Point P, Point A, Point B) { if(A == B) return Length(P-A); // 兩點 Vector v1 = B-A, v2 = P-A, v3 = P-B; if(dcmp(Dot(v1, v2)) < 0) return Length(v2); if(dcmp(Dot(v1, v3)) > 0) return Length(v3); return fabs(Cross(v1, v2) / Length(v1)); } Point GetLineProjection(Point P, Point A, Point B) { Vector v = B-A; return A+v*(Dot(v, P-A) / Dot(v, v)); } bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2) { double c1=Cross(a2-a1,b1-a1), c2=Cross(a2-a1,b2-a1), c3=Cross(b2-b1,a1-b1), c4=Cross(b2-b1,a2-b1); return (dcmp(c1)*dcmp(c2) < 0) && (dcmp(c3)*dcmp(c4) < 0); } bool OnSegment(Point P, Point a1, Point a2) { return dcmp(Cross(a1-P, a2-P)) == 0 && dcmp(Dot(a1-P, a2-P)) < 0; } double PolygonArea(Point* P, int n) { double area = 0; for(int i = 1; i < n-1; i++) area += Cross(P[i]-P[0], P[i+1]-P[0]); return area/2; // 有向面積 } struct Line { Point p; Vector v; Point point(double a) { return p + v*a; } }; const double PI = acos(double(-1)); struct Circle { Point c; double r; Circle(Point c, double r):c(c),r(r) {} Point point(double a) { return Point(c.x + cos(a)*r, c.y + sin(a)*r); } }; int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector& sol) { double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y; double e = a*a + c*c, f = 2 * (a*b+c*d), g = b*b + d*d - C.r*C.r; double delta = f*f - 4*e*g; if(dcmp(delta) < 0) return 0; if(dcmp(delta) == 0) { t1 = t2 = -f / (2*e); sol.push_back(L.point(t1)); return 1; } t1 = (-f - sqrt(delta)) / (2*e); sol.push_back(L.point(t1)); t2 = (-f + sqrt(delta)) / (2*e); sol.push_back(L.point(t2)); return 2; } double angle(Vector v) { return atan2(v.y, v.x); } int getCircleCircleIntersection(Circle C1, Circle C2, vector& sol) { double d = Length(C1.c - C2.c); if(dcmp(d) == 0) { if(dcmp(C1.r-C2.r) == 0) return -1; return 0; } if(dcmp(C1.r+C2.r-d) < 0) return 0; if(dcmp(fabs(C1.r-C2.r)-d > 0)) return 0; double a = angle(C2.c-C1.c); double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d)); Point p1 = C1.point(a-da), p2 = C1.point(a+da); sol.push_back(p1); if(p1 == p2) return 1; sol.push_back(p2); return 2; } int getTangents(Point p, Circle C, Vector* v) { Vector u = C.c-p; double dist = Length(u); if(dist < C.r) return 0; else if(dcmp(dist-C.r) == 0) { v[0] = Rotate(u, PI/2); return 1; } else { double ang = asin(C.r / dist); v[0] = Rotate(u, -ang); v[1] = Rotate(u, +ang); return 2; } } // 可用 dcmp 比較 // 如果不希望在凸包邊上有輸入點, 可以將 <= 換為 < int ConvexHull(Point* p, Point* ch, int n) { sort(p, p+n); // x 為第一關鍵字, y 為第二關鍵字 int m = 0; for(int i = 0; i < n; i++) { while(m > 1 && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--; ch[m++] = p[i]; } int k = m; for(int i = n-2; i >= 0; i--) { while(m > k && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--; ch[m++] = p[i]; } if(n > 1) m--; return m; } bool OnLeft(const Line& L, const Point& p) { return Cross(L.v, p-L.P) > 0; } vector HalfplaneIntersection(vector L) { int n = L.size(); sort(L.begin(), L.end()); // 按極角排序 int first, last; // 雙端佇列的第一個元素和最後一個元素的下標 vector p(n); // p[i]為q[i]和q[i+1]的交點 vector q(n); // 雙端佇列 vector ans; // 結果 q[first=last=0] = L[0]; // 雙端佇列初始化為只有一個半平面L[0] for(int i = 1; i < n; i++) { while(first < last && !OnLeft(L[i], p[last-1])) last--; while(first < last && !OnLeft(L[i], p[first])) first++; q[++last] = L[i]; if(fabs(Cross(q[last].v, q[last-1].v)) < eps) { // 兩向量平行且同向,取內側的一個 last--; if(OnLeft(q[last], L[i].P)) q[last] = L[i]; } if(first < last) p[last-1] = GetLineIntersection(q[last-1], q[last]); } while(first < last && !OnLeft(q[first], p[last-1])) last--; // 刪除無用平面 if(last - first <= 1) return ans; // 空集 p[last] = GetLineIntersection(q[last], q[first]); // 計算首尾兩個半平面的交點 // 從deque複製到輸出中 for(int i = first; i <= last; i++) ans.push_back(p[i]); return ans; } int diameter2(vector& points) { vector p = ConvexHull(points); int n = p.size(); if(n == 1) return 0; if(n == 2) return Dist2(p[0], p[1]); p.push_back(p[0]); // 免得取模 int ans = 0; for(int u = 0, v = 1; u < n; u++) { // 一條直線貼住邊p[u]-p[u+1] while(1) { // 當 Area(p[u], p[u+1], p[v+1]) <= Area(p[u], p[u+1], p[v])時停止旋轉 // 即 Cross(p[u+1]-p[u], p[v+1]-p[u]) - Cross(p[u+1]-p[u], p[v]-p[u]) <= 0 // 根據 Cross(A,B) - Cross(A,C) = Cross(A,B-C) // 化簡得 Cross(p[u+1]-p[u], p[v+1]-p[v]) <= 0 int diff = Cross(p[u+1]-p[u], p[v+1]-p[v]); if(diff <= 0) { ans = max(ans, Dist2(p[u], p[v])); // u和v是對踵點 if(diff == 0) ans = max(ans, Dist2(p[u], p[v+1])); // diff == 0時u和v+1也是對踵點 break; } v = (v + 1) % n; } } return ans; }