算是四个比较重要的计算几何知识点吧,写在这里备忘。
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)$的算法求解半平面。
- 对所有半平面按极角排序。
- 将第一个半平面加入双端队列中。
- 逐个加入半平面,然后分三种情况进行处理。
- 情况1:如果队尾的两个半平面的交点在当前半平面的外面,那么删除队尾半平面。
- 情况2:如果队首的两个半平面的交点在当前半平面的外面,那么删除队首半平面。
- 情况3:如果队尾的两个半平面的交点在队首半平面的外面,那么删除队尾半平面。
- 如果遇到两个半平面的极角相同,那我们应该保留靠左边的那一个。
模板题: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 }