凸包笔记 Convex Hull

发表于2020-07-14更新于2023-09-15字数统计2.4k阅读时长26分

Convex Hull

参考资料
凸包常用的求法有 Graham 扫描法和 Andrew 算法。
Andrew 算法求凸包:
运用向量外积以及单调栈，单调栈例题。
向量类(重载了一些运算符)

Vector
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
struct Vec2 {
    Vec2() = default;
    Vec2(int _x, int _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 int operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x; //×
    }
    int x = 0, y = 0;
} vec[MAXN], cov[MAXN];

Andrew 算法求凸包

Andrew
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
inline double convex(double res = 0) {
    top = 0, up = 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];
    }
    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];
    }
    for (int i = 1; i < top; i++) {  //凸包中的所有点 求周长
        Vec2 tmp = cov[i + 1] - cov[i];
        res += sqrt(pow(tmp.x, 2) + pow(tmp.y, 2));
    }
    return res;
}

例题

A UVA1303 Wall

Problem

题目链接

Solution

先求出给定点的凸包，再将凸包向外推 LL。
顶点处产生圆弧，各顶点圆弧对应圆周角总和为 360°360°(运用多边形内角和可推出)。
答案为凸包长度加上一个半径为L的圆的周长。

Code

Wall
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
#include <bits/stdc++.h>
#define PI acos(-1)

using namespace std;

constexpr int MAXN = 1e3 + 3;

int testCase, n, L, top, up;

struct Vec2 {
    int x, y;

    Vec2() : x(0), y(0) {}
    Vec2(int _x, int _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 int operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x; //×
    }    
} vec[MAXN], cov[MAXN];

inline double convex(double res = 0) {
    top = 0, up = 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];
    }
    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];
    }
    for (int i = 1; i < top; i++) {  //凸包中的所有点
        Vec2 tmp = cov[i + 1] - cov[i];
        res += sqrt(pow(tmp.x, 2) + pow(tmp.y, 2));
    }
    return res;
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cin >> testCase;
    while (testCase--) {
        cin >> n >> L;
        for (int i = 1; i <= n; i++) cin >> vec[i].x >> vec[i].y;
        double res = convex() + PI * 2 * L;
        cout << setprecision(0) << fixed << res << endl;
        if (testCase) cout << endl;
    }
    return 0;
}

B USACO5.1 圈奶牛Fencing the Cows

Problem

题目链接
Andrew 算法模板。

C UVA11626 Convex Hull

Problem

题目链接
Andrew 算法模板。

D UVA811 The Fortified Forest

Problem

题目链接

Solution

N≤15N \le 15，考虑dfs枚举哪些树被砍掉，每个组合跑一次 AndrewAndrew 算法进行比较(注意排序前需拷贝)。

Code

The Fortified Forest
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
#include <bits/stdc++.h>
//旧版凸包写法
using namespace std;

constexpr int MAXN = 20;

struct Vector {
    long long x = 0, y = 0, value = 0, len = 0;

    Vector() = default;
    Vector(long long _x, long long _y, long long _v, long long _l)
        : x(_x), y(_y), value(_v), len(_l) {}

    friend inline istream& operator>>(istream& is, Vector& obj) {
        is >> obj.x >> obj.y >> obj.value >> obj.len;
        return is;
    }

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

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

    friend inline long long operator*(const Vector& a, const Vector& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
} vec[MAXN], vec2[MAXN], vec3[MAXN];

int n, testCase, ansSize, tmpSize;
double ansLen, ansValue;
vector<int> vecTmp, vecAns;

inline void convex() {
    int stk[MAXN] = {0}, used[MAXN] = {0}, top(0);
    double res(0), value(0);
    for (int i = 1; i <= tmpSize; i++) vec2[i] = vec3[i];
    sort(vec2 + 1, vec2 + 1 + tmpSize);
    for (auto &i : vecTmp) value += vec[i].value, res += vec[i].len;
    stk[++top] = 1;
    for (int i = 2; i <= tmpSize; i++) {  //下凸壳
        while (top > 1 &&
               (vec2[stk[top]] - vec2[stk[top - 1]]) * (vec2[i] - vec2[stk[top]]) <= 0)
            used[stk[top--]] = false;
        used[i] = true;
        stk[++top] = i;
    }
    int tmp = top;
    for (int i = tmpSize - 1; i > 0; i--) {  //上凸壳
        if (!used[i]) {
            while (top > tmp && (vec2[stk[top]] - vec2[stk[top - 1]]) *
                                        (vec2[i] - vec2[stk[top]]) <= 0)
                used[stk[top--]] = false;
            used[i] = true;
            stk[++top] = i;
        }
    }
    for (int i = 1; i < top; i++) {
        Vector tmp = vec2[stk[i + 1]] - vec2[stk[i]];
        res -= sqrt(tmp.x * tmp.x + tmp.y * tmp.y);
    }

    if (res >= 0 && (value < ansValue || (value == ansValue && ansSize > vecTmp.size()))) {
        ansLen = res;
        ansValue = value;
        ansSize = vecTmp.size();
        vecAns = vecTmp;
    }
}

void dfs(int cur) {
    if (cur == n + 1) {
        convex();
        return;
    }
    vecTmp.emplace_back(cur);
    dfs(cur + 1); //chopped
    vecTmp.pop_back();

    vec3[++tmpSize] = vec[cur];
    dfs(cur + 1); //not chopped
    --tmpSize;
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    int cnt(0);
    while (cin >> n && n) {
        cnt++;
        if (cnt > 1) cout << endl;
        ansLen = 1e15, ansValue = 1e15, ansSize = 20, tmpSize  = 0;
        vecTmp.clear(), vecAns.clear(), vecTmp.clear();
        for (int i = 1; i <= n; i++) cin >> vec[i];
        dfs(1);
        cout.unsetf(ios_base::fixed);
        cout << "Forest " << cnt << endl;
        cout << "Cut these trees:";
        for (auto &i : vecAns) cout << " " << i;
        cout << endl;
        cout << "Extra wood: " << setprecision(2) << fixed << ansLen << endl;
    }
    return 0;
}

E SHOI2012 信用卡凸包

Problem

题目链接

Solution

将矩形向内缩 1/4 圆半径 r，对其旋转后的顶点求凸包，再外推，加上一个半径为r的圆。

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
#include <bits/stdc++.h>
//旧版凸包写法
using namespace std;

constexpr int MAXN = 1e5 + 5;

struct Vector {
    double x = .0, y = .0;

    Vector() = default;
    Vector(double _x, double _y) : x(_x), y(_y) {}

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

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

    friend inline double operator*(const Vector& a, const Vector& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
} vec[MAXN];

int cnt;

inline double convex() {
    int stk[MAXN] = {0}, used[MAXN] = {0}, top(0);
    memset(used, 0, sizeof(used));
    double res(.0);
    sort(vec + 1, vec + cnt + 1);
    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]]) <= 0)
            used[stk[top--]] = false;
        used[i] = true;
        stk[++top] = i;
    }
    int 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]]) <= 0)
                used[stk[top--]] = false;
            used[i] = true;
            stk[++top] = i;
        }
    }
    for (int i = 1; i < top; i++) { //凸包中的所有点
        Vector tmp = vec[stk[i + 1]] - vec[stk[i]];
        res += sqrt(pow(tmp.x, 2) + pow(tmp.y , 2));
    }   
    return res;
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    int n;
    double a, b, r;
    cin >> n >> a >> b >> r;
    a = a / 2 - r, b = b / 2 - r;
    for (int i = 1; i <= n; i++) {
        double x, y, theta;
        cin >> x >> y >> theta;
        double len = sqrt(pow(a, 2) + pow(b, 2));
        double alpha1 = atan2(a, b) + theta;
        double alpha2 = atan2(a, -b) + theta;
        double alpha3 = atan2(-a, -b) + theta;
        double alpha4 = atan2(-a, b) + theta;
        vec[++cnt] = Vector(x + len * cos(alpha1), y + len * sin(alpha1));
        vec[++cnt] = Vector(x + len * cos(alpha2), y + len * sin(alpha2));
        vec[++cnt] = Vector(x + len * cos(alpha3), y + len * sin(alpha3));
        vec[++cnt] = Vector(x + len * cos(alpha4), y + len * sin(alpha4));
    }
    cout << setprecision(2) << fixed << convex() + acos(-1) * r * 2.0 << endl;
    return 0;
}

F CERC2016 凸轮廓线 Convex Contour

Problem

题目链接

Solution

将正方形和矩形的顶点计入，圆不能仅将上下两个点计入求凸包。
对于圆，可割圆法变为正1e4边形，将顶点计入求凸包。

Code

Convex Contour
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
#include <bits/stdc++.h>

#define PI acos(-1)

using namespace std;

constexpr int MAXN = 2e5 + 5;
constexpr int M = 1e4;

struct Vec2 {
    Vec2() = default;
    Vec2(double _x, 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 double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;  //×
    }
    double x = .0, y = .0;
} vec[MAXN], cov[MAXN];

int cnt;
int top(0), up;

inline double convex(double res = 0) {
    sort(vec + 1, vec + cnt + 1);
    for (int i = 1; i <= cnt; i++) {  //下凸壳
        while (top > 1 && (cov[top] - cov[top - 1]) * (vec[i] - cov[top]) <= 0)
            cov[top--] = Vec2(0, 0);
        cov[++top] = vec[i];
    }
    up = top; //上凸壳起点
    for (int i = cnt - 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];
    }
    for (int i = 1; i < top; i++) {  //凸包中的所有点
        Vec2 tmp = cov[i + 1] - cov[i];
        res += sqrt(pow(tmp.x, 2) + pow(tmp.y, 2));
    }
    return res;
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    double n, st(0);
    string s;
    cin >> n >> s;
    for (int i = 0; i < s.size(); i++) {
        if (s[i] == 'S') {
            vec[++cnt] = Vec2(st, 0);
            vec[++cnt] = Vec2(st, 2);
            vec[++cnt] = Vec2(st + 2, 0);
            vec[++cnt] = Vec2(st + 2, 2);
        } else if (s[i] == 'T') {
            vec[++cnt] = Vec2(st, 0);
            vec[++cnt] = Vec2(st + 2, 0);
            vec[++cnt] = Vec2(st + 1, sqrt(3.0));
        } else {
            for (int j = 1; j <= M; j++) {
                double theta = PI * 2 * j / M;
                vec[++cnt] = Vec2(st + 1 + cos(theta), 1 + sin(theta));
            }
        }
        st += 2;
    }
    cout << setprecision(14) << fixed << convex() / 2 << endl;
    return 0;
}