凸包的直径,旋转卡尺算法

  |  

摘要: 本文介绍计算几何中求凸包的直径的算法:旋转卡尺

【对算法,数学,计算机感兴趣的同学,欢迎关注我哈,阅读更多原创文章】
我的网站:潮汐朝夕的生活实验室
我的公众号:算法题刷刷
我的知乎:潮汐朝夕
我的github:FennelDumplings
我的leetcode:FennelDumplings


各位好,在文章 向量的实现 中,我们实现了一个二维向量的代码模板,主要是 Vector2 这个类。

此后基于二维向量的代码模板,我们解决了计算几何中关于直线和线段极角多边形中的常见问题,以及力扣中的几个相关题目。

在文章 凸包 中,我们学习了凸包的算法,本文进一步看一下怎么用旋转卡尺算法求凸包的直径。

利用凸包上的一些性质,以及旋转卡尺算法,还可以高效求解多边形上的一些问题:例如多边形直径,多边形宽度,最小面积矩形覆盖等。

凸包的直径

凸多边形的直径是指被凸多边形完全包含的最长线段的长度。

朴素算法

枚举每对顶点,找出其中距离最远的一对。

旋转卡壳算法

首先将多边形夹入到两条平行的直线之间。然后沿着多边形的边旋转直线并测得顶点之间的距离。两条直线始终以多边形的两个顶点 a, b 为中心旋转,过程如下图:

旋转卡壳求多边形直径的过程

1
2
3
4
5
首先找出最左端和最右端的点,然后分别向两个顶点加垂线(向量)
左边的记为 a,右边的记为 b 这两个向量方向始终相反。只需要保存第一条的方向向量。

求出当前直线线到下一个顶点之间的角度(共两个)
选择角度更小的直线进行旋转。更新直线方向并将旋转点更换成下一个点。

代码模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
double polygan_diameter(const vector<Vector2>& p)
{
// p 以逆时针给出凸多边形的顶点
int n = p.size();
// 找最左端和最右端的点
int left = min_element(p.begin(), p.end()) - p.begin();
int right = max_element(p.begin(), p.end()) - p.begin();
// 分别向 p[left] 和 p[right] 加垂线,两条垂涎有相反方向
// 标出 a 的方向
Vector2 calipersA(0, 1); // 卡尺的方向向量
double ans = (p[right] - p[left]).norm();
// toNext[i] := p[i] 到下一个顶点方向的单位向量
vector<Vector2> toNext(n);
for(int i = 0; i < n; ++i)
{
// normalize() 返回方向相同的单位向量
toNext[i] = (p[(i + 1) % n] - p[i]).normalize();
}
// a, b 分别表示两条直线正在以哪个顶点为中心旋转
int a = left, b = right;
// 一直执行到两条直线旋转半圈后方向互换为止
while(a != right || b != left)
{
// 查看 a, b 哪个到下一个顶点的角度更小
double cosThetaA = calipersA.dot(toNext[a]);
double cosThetaB = -calipersA.dot(toNext[b]);
if(cosThetaA > cosThetaB)
{
// ThetaA < ThetaB
// [0, PI] 的余弦值递减,不用计算实际角度就可以比较角度的大小关系
calipersA = toNext[a];
a = (a + 1) % n;
}
else
{
calipersA = toNext[b] * (-1);
b = (b + 1) % n;
}
ans = max(ans, (p[a] - p[b]).norm());
}
return ans;
}

模板题:【模板】旋转卡壳

给定平面上 n 个点,求凸包直径。

1
2
3
4
5
6
7
8
9
10
输入格式
第一行一个正整数 n。
接下来 n 行,每行两个整数 ,x,y,表示一个点的坐标。

输出格式
输出一行一个整数,表示答案的平方。

数据范围
2 <= n <= 50000
∣x∣,∣y∣<=1e4

输入输出样例
输入
4
0 0
0 1
1 1
1 0
输出
2

代码 (C++)

本题是洛谷上的题,直接用代码模板 Vector2convex_hullpolygan_diameter 即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
// 求凸包
vector<Vector2> convex_hull(vector<Vector2>& p)
{
int n = p.size();
sort(p.begin(), p.end());
vector<Vector2> result(n * 2);
int k = 0; // 当前的凸包点个数,即 result.size()
// 构造凸包的下侧
for(int i = 0; i < n; ++i)
{
while(k > 1 && (result[k - 1] - result[k - 2]).cross(p[i] - result[k - 1]) < 0)
--k;
result[k++] = p[i];
}
// 构造凸包的上侧
int t = k;
for(int i = n - 2; i >= 0; --i)
{
while(k > t && (result[k - 1] - result[k - 2]).cross(p[i] - result[k - 1]) < 0)
--k;
result[k++] = p[i];
}
result.resize(k - 1);
return result;
}

double polygan_diameter(const vector<Vector2>& p);

int main()
{
int n;
// fstream fin("test/P1452_4.in");
cin >> n;
// fin >> n;
vector<Vector2> p;
for(int i = 0; i < n; ++i)
{
int x, y;
cin >> x >> y;
// fin >> x >> y;
p.emplace_back(x, y);
}
// 处理退化情况:所有点共线
bool flag = true;
for(int i = 0; i < n; ++i)
{
int j = (i + 1) % n;
int k = (i + 2) % n;
if(fabs(ccw(p[i], p[j], p[k])) > EPS)
{
flag = false;
break;
}
}
if(flag)
{
// 所有点共线
Vector2 P = p[0];
double ans1 = 0; // PiP 的部分(PPi.x > 0)
double ans2 = 0; // PPi 的部分(PPi.x < 0)
for(int i = 1; i < n; ++i)
{
Vector2 PPi = p[i] - P;
if(fabs(PPi.x) < EPS)
continue;
else if(PPi.x > EPS)
ans1 = max(ans1, PPi.norm());
else
ans2 = max(ans2, PPi.norm());
}
int ans = round((ans1 + ans2) * (ans1 + ans2));
cout << ans << endl;
return 0;
}
vector<Vector2> ch = convex_hull(p);
double ans = polygan_diameter(ch);
cout << (int)round(ans * ans) << endl;
}

题目

下面是 acwing 上两道旋转卡尺的两道题目,比较难,感兴趣的可以做一下。


Share