模板-计算几何
計算幾何算法相關模版, 可能有錯誤, 省選前持續更正中
重要的不是模版內容, 而是提供算法的實現思路.
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; }vectorHalfplaneIntersection(vectorL) {int n = L.size();sort(L.begin(), L.end()); // 按極角排序int first, last; // 雙端隊列的第一個元素和最后一個元素的下標vectorp(n); // p[i]為q[i]和q[i+1]的交點 vectorq(n); // 雙端隊列 vectorans; // 結果 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) { vectorp = 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; }
總結
