二维凸包 & 旋转卡壳

二维凸包 & 旋转卡壳

参考：

1. 模板

cpp <convex_hull>
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
#include <bits/stdc++.h>
template <class T>
struct Vec2 {
    T x, y;
    Vec2() : Vec2(0, 0) {}
    Vec2(T x, T y) : x(x), y(y) {}
    Vec2(const Vec2& rhs) = default;
    Vec2& operator=(const Vec2& rhs) = default;
    Vec2& operator+=(const Vec2& rhs) { return x += rhs.x, y += rhs.y, *this; }
    Vec2& operator-=(const Vec2& rhs) { return x -= rhs.x, y -= rhs.y, *this; }
    Vec2& operator*=(T rhs) { return x *= rhs, y *= rhs, *this; }
    Vec2& operator/=(T rhs) { return x /= rhs, y /= rhs, *this; }
    Vec2 operator+(const Vec2& rhs) const { return Vec2(*this) += rhs; }
    Vec2 operator-(const Vec2& rhs) const { return Vec2(*this) -= rhs; }
    Vec2 operator*(T rhs) const { return Vec2(*this) *= rhs; }
    Vec2 operator/(T rhs) const { return Vec2(*this) /= rhs; }
    auto operator<=>(const Vec2& rhs) const = default;
    T cross(const Vec2& rhs) const { return x * rhs.y - y * rhs.x; }
    T dot(const Vec2& rhs) const { return x * rhs.x + y * rhs.y; }
    T norm2() const { return dot(*this); }
    T norm() const { return std::sqrt(norm2()); }
    T operator[](int i) const { return i == 0 ? x : y; }
    T& operator[](int i) { return i == 0 ? x : y; }
    friend Vec2 operator*(T lhs, const Vec2& rhs) { return rhs * lhs; }
    friend std::istream& operator>>(std::istream& is, Vec2& rhs) { return is >> rhs.x >> rhs.y; }
    friend std::ostream& operator<<(std::ostream& os, const Vec2& rhs) {
        return os << '(' << rhs.x << ',' << rhs.y << ')';
    }
    Vec2 foot_of_perpendicular(const Vec2& p1, const Vec2& p2) const {
        assert(p1 != p2);
        auto v = p2 - p1;
        return p1 + v * v.dot(*this - p1) / v.norm2();
    }
};
// Returns a list of points on the convex hull in counter-clockwise order.
template <class T>
std::vector<int> convex_hull(std::vector<Vec2<T>>& points, bool contain_points_on_edge = false) {
    int N = points.size();
    if (N == 1) return {0};
    std::sort(points.begin(), points.end());
    std::vector<int> stk(2 * N);
    std::vector<bool> vi(N);
    int k = 0;
    for (int i = 0; i < N; ++i) {
        while (k > 1) {
            auto cross = (points[stk[k - 1]] - points[stk[k - 2]]).cross(points[i] - points[stk[k - 1]]);
            if (contain_points_on_edge ? cross >= 0 : cross > 0) break;
            vi[stk[--k]] = false;
        }
        vi[i] = true;
        stk[k++] = i;
    }
    for (int i = N - 2, t = k + 1; i >= 0; --i) {
        if (i && vi[i]) continue;
        while (k >= t) {
            auto cross = (points[stk[k - 1]] - points[stk[k - 2]]).cross(points[i] - points[stk[k - 1]]);
            if (contain_points_on_edge ? cross >= 0 : cross > 0) break;
            vi[stk[--k]] = false;
        }
        vi[i] = true;
        stk[k++] = i;
    }
    stk.resize(k - 1);
    return stk;
}
// f(i, j, k): j == (i + 1) % chull.size();
// k is the farthest point from the line (points[i], points[j]);
// should be used with `contain_points_on_edge = false`
template <class T, class F>
void rotating_calipers(const std::vector<Vec2<T>>& points, const std::vector<int>& chull, F&& f,
                       bool stop_at_first = false) {
    int N = chull.size();
    assert(N >= 3);
    int k = 2;
    for (int i = 0; i < N; ++i) {
        int j = i + 1 == N ? 0 : i + 1;
        for (;;) {
            int k_ = k + 1 == N ? 0 : k + 1;
            auto old_cross = (points[chull[k]] - points[chull[j]]).cross(points[chull[i]] - points[chull[j]]);
            auto new_cross = (points[chull[k_]] - points[chull[j]]).cross(points[chull[i]] - points[chull[j]]);
            if (stop_at_first ? new_cross <= old_cross : new_cross < old_cross) break;
            k = k_;
        }
        f(chull[i], chull[j], chull[k]);
    }
}
// f(i, j, l, r, k): j == (i + 1) % chull.size();
// k is the farthest point from the line (points[i], points[j]);
// l, r are the left and right most points from the line (points[i], points[j]);
// should be used with `contain_points_on_edge = false`
template <class T, class F>
void rotating_calipers2(const std::vector<Vec2<T>>& points, const std::vector<int>& chull, F&& f,
                        bool stop_at_first = false) {
    int N = chull.size();
    assert(N >= 3);
    int l, r = 1, k = 2;
    for (int i = 0; i < N; ++i) {
        int j = i + 1 == N ? 0 : i + 1;
        for (;;) {
            int k_ = k + 1 == N ? 0 : k + 1;
            auto old_cross = (points[chull[k]] - points[chull[j]]).cross(points[chull[i]] - points[chull[j]]);
            auto new_cross = (points[chull[k_]] - points[chull[j]]).cross(points[chull[i]] - points[chull[j]]);
            if (stop_at_first ? new_cross <= old_cross : new_cross < old_cross) break;
            k = k_;
        }
        for (;;) {
            int r_ = r + 1 == N ? 0 : r + 1;
            auto old_dot = (points[chull[r]] - points[chull[i]]).dot(points[chull[j]] - points[chull[i]]);
            auto new_dot = (points[chull[r_]] - points[chull[i]]).dot(points[chull[j]] - points[chull[i]]);
            if (stop_at_first ? new_dot <= old_dot : new_dot < old_dot) break;
            r = r_;
        }
        if (i == 0) {
            l = r;
            if (stop_at_first) {
                for (;;) {
                    int l_ = l + 1 == N ? 0 : l + 1;
                    auto old_dot = (points[chull[l]] - points[chull[i]]).dot(points[chull[j]] - points[chull[i]]);
                    auto new_dot = (points[chull[l_]] - points[chull[i]]).dot(points[chull[j]] - points[chull[i]]);
                    if (new_dot < old_dot) break;
                    l = l_;
                }
            }
        }
        for (;;) {
            int l_ = l + 1 == N ? 0 : l + 1;
            auto old_dot = (points[chull[l]] - points[chull[j]]).dot(points[chull[i]] - points[chull[j]]);
            auto new_dot = (points[chull[l_]] - points[chull[j]]).dot(points[chull[i]] - points[chull[j]]);
            if (stop_at_first ? new_dot <= old_dot : new_dot < old_dot) break;
            l = l_;
        }
        f(chull[i], chull[j], chull[l], chull[r], chull[k]);
    }
}

2. 例题

2.1 洛谷P2742 二维凸包

给 N 个点，求这些点的凸包的周长。

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include <bits/stdc++.h>
#include <convex_hull>
using namespace std;
int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    int N;
    cin >> N;
    vector<Vec2<double>> points(N);
    for (auto& p : points) cin >> p;
    auto chull = convex_hull(points, false);
    double res = 0;
    for (int i = 0; i < chull.size(); ++i) res += (points[chull[i]] - points[chull[(i + 1) % chull.size()]]).norm();
    cout << fixed << setprecision(2) << res << '\n';
}

2.2 洛谷P1452 求凸包上两点间最大距离

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <bits/stdc++.h>
#include <convex_hull>
using namespace std;
int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    using ll = long long;
    int N;
    cin >> N;
    vector<Vec2<ll>> points(N);
    for (auto& p : points) cin >> p;
    auto chull = convex_hull(points);
    ll max_dist2 = 0;
    if (chull.size() == 2)
        max_dist2 = (points[chull[0]] - points[chull[1]]).norm2();
    else {
        auto f = [&](int i, int j, int k) {
            max_dist2 = max({max_dist2, (points[i] - points[k]).norm2(), (points[j] - points[k]).norm2()});
        };
        rotating_calipers(points, chull, f);
    }
    cout << max_dist2;
}

2.3 洛谷P3187 最小矩形覆盖

cpp
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
#include <bits/stdc++.h>
#include <convex_hull>
using namespace std;
int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    using ld = long double;
    int N;
    cin >> N;
    vector<Vec2<ld>> points(N);
    for (auto& p : points) cin >> p;
    auto chull = convex_hull(points);
    ld min_area = numeric_limits<ld>::infinity();
    array<Vec2<ld>, 4> res;
    auto f = [&](int i_, int j_, int l_, int r_, int k_) {
        auto &i = points[i_], &j = points[j_], &l = points[l_], &r = points[r_], &k = points[k_];
        auto a = l.foot_of_perpendicular(i, j), b = r.foot_of_perpendicular(i, j);
        auto area = (b - a).cross(k - a);
        if (area > min_area) return;
        min_area = area;
        auto h = k - k.foot_of_perpendicular(i, j);
        auto c = b + h, d = a + h;
        res = {a, b, c, d};
    };
    rotating_calipers2(points, chull, f);
    cout << fixed << setprecision(5) << min_area << '\n';
    for (auto& p : res) cout << p.x << ' ' << p.y << '\n';
}
算法 > 计算几何
#算法 #计算几何 #二维凸包 #旋转卡壳
二维凸包 & 旋转卡壳
https://blog.fredbill.eu.org/2024/01/07/算法/计算几何/二维凸包 & 旋转卡壳/
作者
FredBill
发布于
2024年1月7日
许可协议