旋转卡壳笔记 Rotating Calipers

发表于2020-08-08更新于2023-09-15字数统计3.8k阅读时长39分

Rotating Calipers

参考资料

凸多边形的切线

如果一条直线与凸多边形有交点，并且整个凸多边形都在这条直线的一侧，那么这条直线就是该凸多边形的一条切线。

对踵点

如果过凸多边形上两点作一对平行线，使得整个多边形都在这两条线之间，那么这两个点被称为一对对踵点。

凸多边形的直径

即凸多边形上任意两个点之间距离的最大值。直径一定会在对踵点中产生，如果两个点不是对踵点，那么两个点中一定可以让一个点向另一个点的对踵点方向移动使得距离更大。并且点与点之间的距离可以体现为线与线之间的距离，在非对踵点之间构造平行线，一定没有在对踵点构造平行线优，这一点可以通过平移看出。

旋转卡壳


通过求出每一条边的对踵点，即可求出凸多边形的直径。通过上面的图可以看出，边的对踵点同时是这两个点的对踵点，因此并没有遗漏。

具体来讲，通过向量叉积求出边与顶点形成三角形的面积，此三角形的面积也就反映了点到直线的距离。由于面积是单峰的，那么对单独一条边的对踵点，我们可以用三分来求，对于所有边，我们可以用 two−pointertwo-pointer 来 O(n)O(n) 求出。因为随着边在逆时针枚举，它的对踵点 TT 也在逆时针移动。

步骤

首先求出凸包，目的是把所有点按逆时针排序。然后枚举逆时针边，定义 TT 为当前的对踵点，若 TT 逆时针变动能增大到形成的三角形面积改变 TT，直到继续改变会导致距离减小为止。每次求出边（设为 ABAB）的对踵点 TiT_i 后，分别用 ∣ATi∣,∣BTi∣|AT_i|, |BT_i| 来更新答案，即凸多边形的直径。
时间复杂度 O(n)O(n)。

Code

Rotating Calipers
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*Vec2类中operator*重载为了向量叉积 norm为向量取模 cov为凸壳上的点*/
inline long double rotatingCalipers(long double res = 0) { //旋转卡壳
    convex();
    int up = 1, left = 1, right = 1;
    if (top == 3) return normSquare(cov[1] - cov[2]);
    for (int i = 1; i < top; i++) {
        while (true) { //寻找对踵点 two polong doubleers
            auto s = cov[i], t = cov[i + 1];
            auto x1 = cov[up], x2 = cov[up + 1];
            long double s1 = abs((s - x1) * (t - x1));
            long double s2 = abs((s - x2) * (t - x2));
            if (s1 < s2 + EPS) up = up % (top - 1) + 1; //h1 <= h2
            else break;
        }
        res = max(res, normSquare(cov[i] - cov[up]));
        res = max(res, normSquare(cov[i + 1] - cov[up]));
    }   
    return res;
}

下附向量类

向量类
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
struct Vec2 {
    Vec2() = default;
    Vec2(long double _x, long double _y) : x(_x), y(_y) {}

    friend inline Vec2 intersect(const Vec2&, const Vec2&, const Vec2&, const Vec2&);

    friend inline bool operator<(const Vec2& a, const Vec2& b) {
        return a.x == b.x ? a.y < b.y : a.x < b.x;
    }

    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }

    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*@return cross product of two Vec2*/
    friend inline long double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
    /*@return 向量与数乘后的向量*/
    friend inline Vec2 operator*(const Vec2& a, long double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*@param a 2D vector @return norm of the vector*/
    friend inline long double norm(const Vec2& a) {
        return sqrt(a.x * a.x + a.y * a.y);
    }
    /*@param a 2D vector @return square norm of the vector*/
    friend inline long double normSquare(const Vec2& a) {
        return a.x * a.x + a.y * a.y;
    }
    long double x = 0, y = 0;
} vec[MAXN];

/*@return intersect point of two Seg*/
inline Vec2 intersect(const Vec2& s1, const Vec2& t1,
                        const Vec2& s2, const Vec2& t2) {
    long double ratio = ((s1 - s2) * (t2 - s2)) / ((t2 - s2) * (t1 - s1));
    return s1 + (t1 - s1) * ratio;
}

例题

A Beauty Contest G

Description

模板题

Solution

模板题。不过求的是模长的平方。

Code

Beauty Contest G
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include <bits/stdc++.h>
#define PI acos(-1)

using namespace std;
//小数旋转卡壳
constexpr int MAXN = 2e5 + 5;
constexpr long double EPS = 1e-10;

struct Vec2 {
    Vec2() = default;
    Vec2(long double _x, long double _y) : x(_x), y(_y) {}

    friend inline Vec2 intersect(const Vec2&, const Vec2&, const Vec2&, const Vec2&);

    friend inline bool operator<(const Vec2& a, const Vec2& b) {
        return a.x == b.x ? a.y < b.y : a.x < b.x;
    }

    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }

    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*@return cross product of two Vec2*/
    friend inline long double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
    /*@return 向量与数乘后的向量*/
    friend inline Vec2 operator*(const Vec2& a, long double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*@param a 2D vector @return norm of the vector*/
    friend inline long double norm(const Vec2& a) {
        return sqrt(a.x * a.x + a.y * a.y);
    }
    /*@param a 2D vector @return square norm of the vector*/
    friend inline long double normSquare(const Vec2& a) {
        return a.x * a.x + a.y * a.y;
    }

    long double x = 0, y = 0;
} vec[MAXN], cov[MAXN];

/*@return intersect point of two Seg*/
inline Vec2 intersect(const Vec2& s1, const Vec2& t1,
                        const Vec2& s2, const Vec2& t2) {
    long double ratio = ((s1 - s2) * (t2 - s2)) / ((t2 - s2) * (t1 - s1));
    return s1 + (t1 - s1) * ratio;
}

int n, top;

inline void convex() {
    top = 0;
    sort(vec + 1, vec + n + 1);
    for (int i = 1; i <= n; i++) {  //下凸壳
        while (top > 1 && (cov[top] - cov[top - 1]) * (vec[i] - cov[top]) <= 0)
            cov[top--] = Vec2(0, 0);
        cov[++top] = vec[i];
    }
    int up = top; //上凸壳起点
    for (int i = n - 1; i > 0; i--) {  //上凸壳
        while (top > up && (cov[top] - cov[top - 1]) * (vec[i] - cov[top]) <= 0)
            cov[top--] = Vec2(0, 0);
        cov[++top] = vec[i];
    }
}

inline long double rotatingCalipers(long double res = 0) { //旋转卡壳
    convex();
    int up = 1, left = 1, right = 1;
    if (top == 3) return normSquare(cov[1] - cov[2]);
    for (int i = 1; i < top; i++) {
        while (true) { //寻找对踵点 two polong doubleers
            auto s = cov[i], t = cov[i + 1];
            auto x1 = cov[up], x2 = cov[up + 1];
            long double s1 = abs((s - x1) * (t - x1));
            long double s2 = abs((s - x2) * (t - x2));
            if (s1 < s2 + EPS) up = up % (top - 1) + 1; //h1 <= h2
            else break;
        }
        res = max(res, normSquare(cov[i] - cov[up]));
        res = max(res, normSquare(cov[i + 1] - cov[up]));
    }   
    return res;
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> vec[i].x >> vec[i].y;
    cout << setprecision(0) << fixed << rotatingCalipers() << endl; //返回最远点对距离的平方
    return 0;
}

B HNOI2007 最小矩形覆盖

Description

题目链接
给定一些点的坐标，要求求能够覆盖所有点的最小面积的矩形，输出所求矩形的面积和四个顶点坐标。

Solution

首先最小矩形覆盖下的矩形的一条边一定在凸包上。与旋转卡壳求凸多边形直径类似，利用单峰性求出距离直线最远的顶点 upup，最靠右的点 rightright (通过向量点 >0>0 积判断)，最靠左的点 leftleft (通过向量点积 <0<0 判断)。通过这三个点和选取的凸包上的一条边可以求出矩形四个顶点，进而求出最小矩形覆盖面积。
顶点求法：由 up,left,rightup, left, right 对应的矩形边以及选取的直线 ABAB 相交可得(实现起来略微复杂，两线段对应直线的交点见 intersect()intersect() 函数)。

Code

最小矩形覆盖
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include <bits/stdc++.h>
#define PI acos(-1)

using namespace std;

constexpr int MAXN = 2e5 + 5;
constexpr long double EPS = 1e-10;

struct Vec2 {
    Vec2() = default;
    Vec2(long double _x, long double _y) : x(_x), y(_y) {}

    friend inline bool operator<(const Vec2& a, const Vec2& b) {
        return a.x == b.x ? a.y < b.y : a.x < b.x;
    }

    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }

    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*@return cross product of two Vec2*/
    friend inline long double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
    /*@return 向量与数乘后的向量*/
    friend inline Vec2 operator*(const Vec2& a, long double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*@param a 2D vector @return norm of the vector*/
    friend inline long double norm(const Vec2& a) {
        return sqrt(a.x * a.x + a.y * a.y);
    }
    /*@param a 2D vector @return square norm of the vector*/
    friend inline long double normSquare(const Vec2& a) {
        return a.x * a.x + a.y * a.y;
    }
    /*@return dot product of two vectors*/
    friend inline long double dot(const Vec2& a, const Vec2& b) {
        return a.x * b.x + a.y * b.y;
    }
    long double x = 0, y = 0;
} vec[MAXN], ans[4];

/*@return intersect point of two Seg*/
inline Vec2 intersect(const Vec2& s1, const Vec2& t1, const Vec2& s2,
                      const Vec2& t2) {
    long double ratio = ((s1 - s2) * (t2 - s2)) / ((t2 - s2) * (t1 - s1));
    return s1 + (t1 - s1) * ratio;
}

int cnt;
int stk[MAXN], used[MAXN], top;

inline void convex() {  //求凸包 旧版写法
    sort(vec + 1, vec + cnt + 1);
    top = 0, stk[++top] = 1;
    for (int i = 2; i <= cnt; i++) {  //下凸壳
        while (top > 1 &&
               (vec[stk[top]] - vec[stk[top - 1]]) * (vec[i] - vec[stk[top]]) <=
                   EPS)
            used[stk[top--]] = false;
        used[i] = true, stk[++top] = i;
    }
    long double tmp = top;
    for (int i = cnt - 1; i > 0; i--) {  //上凸壳
        if (!used[i]) {
            while (top > tmp && (vec[stk[top]] - vec[stk[top - 1]]) *
                                        (vec[i] - vec[stk[top]]) <=
                                    EPS)
                used[stk[top--]] = false;
            used[i] = true, stk[++top] = i;
        }
    }
}

inline void rotatingCalipers(long double res = 1e18) { //旋转卡壳
    convex();
    int up = 1, left = 3, right = 1;
    for (int i = 1; i < top; i++) {
        auto s = vec[stk[i]], t = vec[stk[i + 1]];
        while (true) {  //寻找对踵点 two polong doubleers
            auto x1 = vec[stk[up]], x2 = vec[stk[up + 1]];
            long double s1 = abs((s - x1) * (t - x1));
            long double s2 = abs((s - x2) * (t - x2));
            if (s1 < s2 + EPS) up = up % (top - 1) + 1;  // h1 <= h2
            else break;
        }
        while (true) {  //右边最远的点
            auto x1 = vec[stk[right]], x2 = vec[stk[right + 1]];
            long double s1 = dot(x1 - s, t - s);
            long double s2 = dot(x2 - s, t - s);
            if (s1 < s2 + EPS) right = right % (top - 1) + 1;  // r1 <= r2
            else break;
        }
        while (true) {  //左边最远的点
            auto x1 = vec[stk[left]], x2 = vec[stk[left + 1]];
            long double s1 = dot(x1 - s, t - s);
            long double s2 = dot(x2 - s, t - s);
            if (s1 > EPS) left = left % (top - 1) + 1;
            else if (s2 < s1 + EPS) left = left % (top - 1) + 1;  // l2 <= l1
            else break;
        }
        long double theta = atan2(t.y - s.y, t.x - s.x);
        auto p1 = intersect(
            s, t, vec[stk[left]],
            vec[stk[left]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        auto p2 = intersect(
            s, t, vec[stk[right]],
            vec[stk[right]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        auto p3 = intersect(
            vec[stk[up]], vec[stk[up]] + Vec2(cos(theta), sin(theta)),
            vec[stk[right]],
            vec[stk[right]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        auto p4 = intersect(
            vec[stk[up]], vec[stk[up]] + Vec2(cos(theta), sin(theta)),
            vec[stk[left]],
            vec[stk[left]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        long double area = (p3 - p2) * (p1 - p2);
        if (area < res - EPS) { // area < res
            res = area;
            ans[0] = p1, ans[1] = p2, ans[2] = p3, ans[3] = p4;
        }
    }
    cout << res << endl;
}

inline void print() {
    int st;
    long double minY = 1e18, minX = 1e18;
    for (int i = 0; i < 4; i++) {
        if (ans[i].y < minY - EPS) {
            st = i, minY = ans[i].y, minX = ans[i].y;
        } else if (abs(ans[i].y - minY) < EPS) { //相等
            if (ans[i].x < minX - EPS) { //选x小的
                st = i, minY = ans[i].y, minX = ans[i].x;
            }
        }
    }
    for (int i = 0; i < 4; i++) {
        double x = ans[(st + i) % 4].x, y = ans[(st + i) % 4].y;
        if (abs(x) < EPS) x = .0;
        if (abs(y) < EPS) y = .0;
        cout << x << " " << y << endl;
    }
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cout << setprecision(5) << fixed;
    cin >> cnt;
    for (int i = 1; i <= cnt; i++) cin >> vec[i].x >> vec[i].y;
    rotatingCalipers(); //旋转卡壳求出矩形
    print(); //输出
    return 0;
}

C UVA10173 Smallest Bounding Rectangle

Description

题目链接
多组数据求最小矩形覆盖的面积。

Solution

与 B HNOI2007 最小矩形覆盖做法相同，只需输出面积。
注意多组数据 usedused 数组清零。

Code

UVA10173 Smallest Bounding Rectangle
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#include <bits/stdc++.h>
#define PI acos(-1)

using namespace std;

constexpr int MAXN = 2e5 + 5;
constexpr long double EPS = 1e-11;

struct Vec2 {
    Vec2() = default;
    Vec2(long double _x, long double _y) : x(_x), y(_y) {}

    friend inline bool operator<(const Vec2& a, const Vec2& b) {
        return a.x == b.x ? a.y < b.y : a.x < b.x;
    }

    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }

    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*@return cross product of two Vec2*/
    friend inline long double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
    /*@return 向量与数乘后的向量*/
    friend inline Vec2 operator*(const Vec2& a, long double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*@param a 2D vector @return norm of the vector*/
    friend inline long double norm(const Vec2& a) {
        return sqrt(a.x * a.x + a.y * a.y);
    }
    /*@param a 2D vector @return square norm of the vector*/
    friend inline long double normSquare(const Vec2& a) {
        return a.x * a.x + a.y * a.y;
    }
    /*@return dot product of two vectors*/
    friend inline long double dot(const Vec2& a, const Vec2& b) {
        return a.x * b.x + a.y * b.y;
    }
    long double x = 0, y = 0;
} vec[MAXN];

/*@return intersect point of two Seg*/
inline Vec2 intersect(const Vec2& s1, const Vec2& t1, const Vec2& s2,
                      const Vec2& t2) {
    long double ratio = ((s1 - s2) * (t2 - s2)) / ((t2 - s2) * (t1 - s1));
    return s1 + (t1 - s1) * ratio;
}

int cnt;
int stk[MAXN], used[MAXN], top;

inline void convex() {  //求凸包
    sort(vec + 1, vec + cnt + 1);
    top = 0, stk[++top] = 1;
    memset(used, 0, sizeof(used));
    for (int i = 2; i <= cnt; i++) {  //下凸壳
        while (top > 1 &&
               (vec[stk[top]] - vec[stk[top - 1]]) * (vec[i] - vec[stk[top]]) < EPS)
            used[stk[top--]] = false;
        used[i] = true, stk[++top] = i;
    }
    long double tmp = top;
    for (int i = cnt - 1; i > 0; i--) {  //上凸壳
        if (!used[i]) {
            while (top > tmp && (vec[stk[top]] - vec[stk[top - 1]]) *
                                        (vec[i] - vec[stk[top]]) < EPS)
                used[stk[top--]] = false;
            used[i] = true, stk[++top] = i;
        }
    }
}

inline void rotatingCalipers(long double res = 1e18) { //旋转卡壳
    convex();
    int up = 2, left = 2, right = 1;
    if (top <= 3) { cout << .0 << endl; return; }
    for (int i = 1; i < top; i++) {
        auto s = vec[stk[i]], t = vec[stk[i + 1]];
        while (true) {  //寻找对踵点 two polong doubleers
            auto x1 = vec[stk[up]], x2 = vec[stk[up + 1]];
            long double s1 = abs((s - x1) * (t - x1));
            long double s2 = abs((s - x2) * (t - x2));
            if (s1 < s2 + EPS) up = up % (top - 1) + 1;  // h1 <= h2
            else break;
        }
        while (true) {  //右边最远的点
            auto x1 = vec[stk[right]], x2 = vec[stk[right + 1]];
            long double s1 = dot(x1 - s, t - s);
            long double s2 = dot(x2 - s, t - s);
            if (s1 < s2 + EPS) right = right % (top - 1) + 1;  // r1 <= r2
            else break;
        }
        while (true) {  //左边最远的点
            auto x1 = vec[stk[left]], x2 = vec[stk[left + 1]];
            long double s1 = dot(x1 - s, t - s);
            long double s2 = dot(x2 - s, t - s);
            if (s1 > EPS) left = left % (top - 1) + 1;
            else if (s2 < s1 + EPS) left = left % (top - 1) + 1;  // l2 <= l1
            else break;
        }
        long double theta = atan2(t.y - s.y, t.x - s.x);
        auto p1 = intersect( //直线相交求顶点
            s, t, vec[stk[left]],
            vec[stk[left]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        auto p2 = intersect( //直线相交求顶点
            s, t, vec[stk[right]],
            vec[stk[right]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        auto p3 = intersect( //直线相交求顶点
            vec[stk[up]], vec[stk[up]] + Vec2(cos(theta), sin(theta)),
            vec[stk[right]],
            vec[stk[right]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        auto p4 = intersect( //直线相交求顶点
            vec[stk[up]], vec[stk[up]] + Vec2(cos(theta), sin(theta)),
            vec[stk[left]],
            vec[stk[left]] + Vec2(cos(theta + PI / 2), sin(theta + PI / 2)));
        long double area = (p3 - p2) * (p1 - p2);
        res = min(res, area);
    }
    cout << res << endl;
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cout << setprecision(4) << fixed;
    while (cin >> cnt && cnt) {
        for (int i = 1; i <= cnt; i++) cin >> vec[i].x >> vec[i].y;
        rotatingCalipers(); //旋转卡壳求出矩形
    }
    return 0;
}