数几何

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

字数统计: 1.4k字   |   阅读时长: 6分
,

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

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的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++)

本题是洛谷上的题，直接用代码模板 Vector2、convex_hull 和 polygan_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