算是四个比较重要的计算几何知识点吧,写在这里备忘。
PS:本文章默认读者了解最基本的计算几何知识,所以不提供任何模板。
凸包:
凸包的概念:给定$n(n≥3)$个点,求能把这些点包含在内的面积最小的多边形。
一般我们使用Andrew算法求解凸包,时间复杂度为$O(nlogn)$。
Andrew算法的思路很简单,就是正反做两次扫描,分别构造下凸包和上凸包。
考虑构造下凸包的过程,我们维护一个栈,假设我们现在已经有3个点$p_1, p_2, p_3$,现在发现$p_4p_3$对于$p_3p_2$是右拐弯的,那说明$p_3$不在凸包上,于是将$p_3$出栈,继续对$p_2$进行检查。
构造上凸包也是同理,每个点最多进出一次栈,所以扫描的复杂度是$O(n)$,由于我们还要对所有的点进行一次升序排序($x$坐标为第一关键字,$y$坐标为第二关键字),所以总的复杂度为$O(nlogn)$。
模板题:Luogu P2742
1 #include <bits/stdc++.h> 2 using namespace std; 3 4 const double EPS = 1e-8; 5 const int MAXN = 10010; 6 7 struct Point 8 { 9 double x, y; 10 Point(double X = 0.0, double Y = 0.0){x = X, y = Y;} 11 Point operator + (Point B) {return Point(x + B.x, y + B.y);} 12 Point operator - (Point B) {return Point(x - B.x, y - B.y);} 13 bool operator == (Point B) {return fabs(x - B.x) < EPS && fabs(y - B.y) < EPS;} 14 }; 15 typedef Point Vector; 16 inline double Cross(Vector A, Vector B) {return A.x * B.y - A.y * B.x;} 17 inline double Dis(Point A, Point B) {return hypot(A.x - B.x, A.y - B.y);} 18 inline bool cmp(const Point &A, const Point &B) {return fabs(A.x - B.x) < EPS ? A.y - B.y < EPS : A.x - B.x < EPS;} 19 20 int n, cnt; double res; 21 Point p[MAXN], ans[MAXN]; 22 23 void Convex() 24 { 25 sort(p + 1, p + 1 + n, cmp); 26 n = unique(p + 1, p + 1 + n) - p - 1; 27 for(int i = 1; i <= n; ++i) 28 { 29 while(cnt > 1 && Cross(ans[cnt] - ans[cnt - 1], p[i] - ans[cnt - 1]) < EPS) --cnt; 30 ans[++cnt] = p[i]; 31 } 32 int j = cnt; 33 for(int i = n - 1; i; --i) 34 { 35 while(cnt > j && Cross(ans[cnt] - ans[cnt - 1], p[i] - ans[cnt - 1]) < EPS) --cnt; 36 ans[++cnt] = p[i]; 37 } 38 if(n > 1) --cnt; 39 } 40 41 int main() 42 { 43 while(scanf("%d", &n) != EOF && n) 44 { 45 for(int i = 1; i <= n; ++i) scanf("%lf %lf", &p[i].x, &p[i].y); 46 cnt = 0, res = 0.0; Convex(); 47 if(cnt == 2) res = Dis(ans[1], ans[2]); 48 else if(cnt >= 3) 49 { 50 res = Dis(ans[cnt], ans[1]); 51 for(int i = 1; i < cnt; ++i) res += Dis(ans[i], ans[i + 1]); 52 } 53 printf("%.2f\n", res); 54 } 55 return 0; 56 }
最近点对:
暴力是$O(n^2)$的,考虑分治(分治前先进行排序,排序方式参考凸包算法)。
当目前区间只有一个点或者两个点时,直接计算距离即可。
现在问题是怎样合并两个区间,假设两个区间分别为$S_1, S_2$。
令$S_1, S_2$中的最近距离分别为$d_1, d_2$。
那我们就在中间点$p[mid]$附近找到所有离它的距离小于$min(d_1, d_2)$的点,令这些点的集合为$tmp$。
然后我们按照$y$坐标对$tmp$中的点进行排序,之后暴力枚举$tmp$中的点对,并利用$y$坐标的单调性进行剪枝。
可以证明这样的复杂度是$O(nlogn)$的。
模板题:Luogu P1429
1 #include <bits/stdc++.h> 2 using namespace std; 3 4 const int MAXN = 200010; 5 const double EPS = 1e-8; 6 7 struct Point 8 { 9 double x, y; 10 }; 11 inline double Dis(Point A, Point B) {return hypot(A.x - B.x, A.y - B.y);} 12 inline bool cmp1(Point A, Point B) {return fabs(A.x - B.x) < EPS ? A.y - B.y < EPS : A.x - B.x < EPS;} 13 inline bool cmp2(Point A, Point B) {return A.y - B.y < EPS;} 14 15 int n; 16 Point p[MAXN], tmp[MAXN]; 17 18 double solve(int l, int r) 19 { 20 double dis = 1e20; 21 if(l == r) return dis; 22 if(l + 1 == r) return Dis(p[l], p[r]); 23 int mid = (l + r) >> 1; 24 double d1 = solve(l, mid); 25 double d2 = solve(mid + 1, r); 26 dis = min(d1, d2); 27 int cnt = 0; 28 for(int i = l; i <= r; ++i) 29 if(fabs(p[mid].x - p[i].x) <= dis) 30 tmp[++cnt] = p[i]; 31 sort(tmp + 1, tmp + 1 + cnt, cmp2); 32 for(int i = 1; i <= cnt; ++i) 33 for(int j = i + 1; j <= cnt; ++j) 34 { 35 if(tmp[j].y - tmp[i].y >= dis) break; 36 dis = min(dis, Dis(tmp[j], tmp[i])); 37 } 38 return dis; 39 } 40 41 int main() 42 { 43 scanf("%d", &n); 44 for(int i = 1; i <= n; ++i) scanf("%lf %lf", &p[i].x, &p[i].y); 45 sort(p + 1, p + 1 + n, cmp1); 46 printf("%.4f\n", solve(1, n)); 47 return 0; 48 }
旋转卡壳:
首先你要会凸包。
然后考虑一个问题,平面最远点对。
首先最远点对肯定是在凸包上,所以先求出凸包。
我们称两条平行线与凸包的交点为对踵点对,可以证明对踵点对的数量是$O(n)$的。
根据下图,容易发现对踵点具有一个特别的性质。
简而言之,就是这个三角形的面积变化是具有单峰函数的性质的。
那我们逆时针枚举每个点,计算这个点的对踵点,我们只要接着上一个点的对踵点逆时针枚举即可。
发现这样就相当于是用两条平行线卡着凸包的外壳逆时针转了一圈,这也是算法名称的来源。
模板题:POJ 2187
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cmath> 5 #include<algorithm> 6 #include<vector> 7 #include<set> 8 using namespace std; 9 10 const int MAXN = 50010; 11 12 struct Point 13 { 14 int x, y; 15 Point(int X = 0, int Y = 0){x = X, y = Y;} 16 Point operator + (Point B) {return Point(x + B.x, y + B.y);} 17 Point operator - (Point B) {return Point(x - B.x, y - B.y);} 18 }; 19 inline bool cmp(const Point &A, const Point &B) {return A.x == B.x ? A.y < B.y : A.x < B.x;} 20 inline int dis(Point A, Point B) {return (A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y);} 21 inline int Cross(Point A, Point B) {return A.x * B.y - A.y * B.x;} 22 23 int n, m, ans; 24 Point p[MAXN], q[MAXN]; 25 26 void Convex_hull() 27 { 28 sort(p, p + n, cmp); 29 for(int i = 0; i < n; ++i) 30 { 31 while(m > 1 && Cross(q[m - 1] - q[m - 2], p[i] - q[m - 1]) <= 0) --m; 32 q[m++] = p[i]; 33 } 34 for(int i = n - 2, j = m; ~i; --i) 35 { 36 while(m > j && Cross(q[m - 1] - q[m - 2], p[i] - q[m - 1]) <= 0) --m; 37 q[m++] = p[i]; 38 } 39 if(n > 1) --m; 40 } 41 42 int solve() 43 { 44 if(m == 2) return dis(q[0], q[1]); 45 for(int i = 0, j = 1; i < m; ++i) 46 { 47 while(abs(Cross(q[(i + 1) % m] - q[i], q[j] - q[i])) < 48 abs(Cross(q[(i + 1) % m] - q[i], q[(j + 1) % m] - q[i]))) j = (j + 1) % m; 49 ans = max(ans, dis(q[i], q[j])); 50 } 51 return ans; 52 } 53 54 int main() 55 { 56 scanf("%d", &n); 57 for(int i = 0; i < n; ++i) scanf("%d %d", &p[i].x, &p[i].y); 58 Convex_hull(); printf("%d\n", solve()); 59 return 0; 60 }
半平面交:
高中数学必修5?
我们用一条有向直线来表示一个半平面(直线的左边是半平面)
为了方便,我们可能会人为的添加几个半平面,以避免最后半平面交不闭合的情况。
下面给出一种复杂度为$O(nlogn)$的算法求解半平面。
模板题:Luogu P4196
1 #include<bits/stdc++.h> 2 using namespace std; 3 4 const int MAXN = 510; 5 const double EPS = 1e-8; 6 inline int sgn(double x) {return fabs(x) < EPS ? 0 : (x > 0 ? 1 : -1);} 7 8 struct Point 9 { 10 double x, y; 11 Point operator + (Point B) {return (Point){x + B.x, y + B.y};} 12 Point operator - (Point B) {return (Point){x - B.x, y - B.y};} 13 Point operator * (double k) {return (Point){x * k, y * k};} 14 Point operator / (double k) {return (Point){x / k, y / k};} 15 }; 16 typedef Point Vector; 17 inline double Cross(Point A, Point B) {return A.x * B.y - A.y * B.x;} 18 19 struct Line 20 { 21 Point p, v; 22 double ang; 23 }; 24 inline bool cmp(Line A, Line B) {return sgn(A.ang - B.ang) == 0 ? sgn(Cross(A.v - A.p, B.v - A.p)) > 0 : sgn(A.ang - B.ang) < 0;} 25 26 Point Cross_point(Line A, Line B) 27 { 28 Vector v1 = A.v - A.p, v2 = B.v - B.p, u = B.p - A.p; 29 double t = Cross(u, v1) / Cross(v1, v2); 30 return B.p + v2 * t; 31 } 32 33 inline bool OnRight(Line A, Line B, Line C) 34 { 35 Point p = Cross_point(A, B); 36 return sgn(Cross(C.v - C.p, p - C.p)) < 0; 37 } 38 39 Line a[MAXN], q[MAXN]; 40 Point b[MAXN]; 41 int n, tot, cnt; 42 43 void HPI() 44 { 45 sort(a + 1, a + 1 + tot, cmp); 46 for(int i = 1; i <= tot; ++i) 47 { 48 if(sgn(a[i].ang - a[i - 1].ang) != 0) ++cnt; 49 a[cnt] = a[i]; 50 } 51 int l = 1, r = 0; 52 q[++r] = a[1], q[++r] = a[2]; 53 for(int i = 3; i <= cnt; ++i) 54 { 55 while(l < r && OnRight(q[r - 1], q[r], a[i])) --r; 56 while(l < r && OnRight(q[l + 1], q[l], a[i])) ++l; 57 q[++r] = a[i]; 58 } 59 while(l < r && OnRight(q[r - 1], q[r], q[l])) --r; 60 while(l < r && OnRight(q[l + 1], q[l], q[r])) ++l; 61 q[r + 1] = q[l]; tot = 0; 62 for(int i = l; i <= r; ++i) b[++tot] = Cross_point(q[i], q[i + 1]); 63 } 64 65 int main() 66 { 67 scanf("%d", &n); 68 for(int i = 1; i <= n; ++i) 69 { 70 int k; scanf("%d", &k); 71 for(int j = 1; j <= k; ++j) scanf("%lf %lf", &b[j].x, &b[j].y); 72 b[k + 1] = b[1]; 73 for(int j = 1; j <= k; ++j) a[++tot].p = b[j], a[tot].v = b[j + 1]; 74 } 75 for(int i = 1; i <= tot; ++i) a[i].ang = atan2(a[i].v.y - a[i].p.y, a[i].v.x - a[i].p.x); 76 HPI(); b[tot + 1] = b[1]; double ans = 0.0; 77 if(tot > 2) for(int i = 1; i <= tot; ++i) ans += Cross(b[i], b[i + 1]); 78 printf("%.3f\n", fabs(ans) / 2.0); 79 return 0; 80 }
原文:https://www.cnblogs.com/Aegir/p/11311843.html