mydcwfy's Blog

凸包与旋转卡壳

Posted on
Symbols count in article: 13k Reading time ≈ 28 mins.

前置知识：计算几何。

1. 凸包

直接理解为我们用一条橡皮筋围住这些点。

我们就不再介绍 Javis 算法，直接介绍两个时间复杂度更优的 Graham 和 Andrew 算法。

1）Graham 算法

我们首先选择最左下方的点，注意到这个点一定是凸包上的点。假设这个点为 A。

接着，我们按照 A 点与这些点的连线与水平线的夹角排序。注意到如果夹角相同，将最远的点放在最前面。

接着，我们扫描整个数组，同时维护一个栈，记录当前的凸包里的点。

我们要加入一个点的时候，如果当前的栈里的点会和这个点形成一个顺时针的转角的话，我们就将栈顶弹出。

看到这个图，我们现在要将红线上面的点加入，容易发现红线下面这个点是一定不需要的，我们直接弹出即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
bool cmp(const Point &t1, const Point &t2)
{
    return (t1 - p[1]) * (t2 - p[1]) > eps || (fabs((t1 - p[1]) * (t2 - p[1])) < eps && dist(p[1], t1) > dist(p[1], t2));
}

int t = 1;
for (int i = 2; i <= n; ++ i)
    if (p[i].y < p[t].y || (fabs(p[i].y - p[t].y) < eps && p[i].x < p[t].x)) t = i;
swap(p[1], p[t]);
sort(p + 2, p + n + 1, cmp);
stk[++ top] = p[1], stk[++ top] = p[2];
for (int i = 3; i <= n; ++ i)
{
    while (top > 1 && (stk[top] - stk[top - 1]) * (p[i] - stk[top]) < -eps) top --;
    stk[++ top] = p[i];
}
stk[++ top] = p[1];//第一个点进入两次，便于计算周长面积

2）Andrew 算法

我们不再按照夹角排序，直接使用 x 坐标排序，结果又如何呢？

我们这么遍历，发现只会走到一半的凸壳。

于是，我们再会过来遍历一次，就可以把另一半的凸壳遍历到了！

这个算法避免的一些复杂的夹角的计算，常数略小。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
bool cmp(const Point &p1, const Point &p2){return (p1.x == p2.x) ? (p1.y < p2.y) : p1.x < p2.x;}

void Andrew(int n)
{
	sort(p + 1, p + n + 1, cmp);
	for (int i = 1; i <= n; ++ i)
	{
		while (top > 1 && (p[stk[top]] - p[stk[top - 1]]) * (p[i] - p[stk[top]]) <= 0)
			if ((p[stk[top]] - p[stk[top - 1]]) * (p[i] - p[stk[top]]) < 0) usd[stk[top --]] = false;
			else top --;
		usd[stk[++ top] = i] = true;
	}
	usd[1] = false;
	for (int i = n; i; -- i)
	{
		if (usd[i]) continue;
		while (top > 1 && (p[stk[top]] - p[stk[top - 1]]) * (p[i] - p[stk[top]]) <= 0)
			if ((p[stk[top]] - p[stk[top - 1]]) * (p[i] - p[stk[top]]) < 0) usd[stk[top --]] = false;
			else top --;
		usd[stk[++ top] = i] = true;
	}
}

3）例题

T1：信用卡凸包

题目传送门 Luogu

题目传送门 AcWing

我们可以发现，一个凸多边形的外角是 $360^{\circ}$。

那么，圆形一定贡献的是一个周角的大小。剩下的，我们发现可以向内平移到圆心的位置。按圆心凸包即可。

可以发现，黑直线的长度就是红线的长度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
int main()
{
	int n;
	scanf("%d %lf %lf %lf", &n, &a, &b, &R);
	a -= 2 * R, b -= 2 * R;
	double x, y, th;
	for (int i = 1; i <= n; ++ i)
	{
		scanf("%lf %lf %lf", &x, &y, &th);
		x += eps, y += eps, th += eps;
		p[(i << 2) - 3] = (Point){x, y} + rotate({b / 2, a / 2}, th);
		p[(i << 2) - 2] = (Point){x, y} + rotate({-b / 2, a / 2}, th);
		p[(i << 2) - 1] = (Point){x, y} + rotate({-b / 2, -a / 2}, th);
		p[(i << 2)] = (Point){x, y} + rotate({b / 2, -a / 2}, th);
	}
	Andrew(n << 2);
	double res = 0;
	for (int i = 1; i < top; ++ i)
		res += dist(p[stk[i]], p[stk[i + 1]]);
	printf("%.2lf", res + 2 * acos(-1) * R);
	return 0;
}

2. 半平面交

我们介绍一种办法（我似乎并不知道是什么算法），能在 $O(n \log n)$ 的时间求出围住的凸多边形（注意好像不能判断无解的情况）。

为了方便，我们定义半平面存储为一条有向直线，在这条直线逆时针（可以理解为左边）的部分即为半平面。

1）算法流程

我们首先按照每条线与 $x$ 轴的夹角排序，如果相同，说明是平行的，我们按照从左至右的顺序。

然后，我们顺次插入每一个半平面。我们维护一个如此的双端队列来表示当前的半平面。

接着，我们插入的时候，看一下有没有在当前半平面外面的交点，也就是判断有没有在直线顺时针的点，有的话就删除。

比如蓝点在新加入的蓝线右边，所以肯定会导致最后加入的黑线被弹出队尾。

注意有可能线的方向不定导致会弹出队首，所以前后都要判断。

最后的时候，我们再用队首来尝试弹出队尾，队尾尝试弹出队首。

注意我们求夹角的时候，可以使用 atan2(x, y)，这样可以防止出现 nan 的情况，因为他会判断 $y$ 是否等于 0。

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
Point get_point_(Point p1, Point k1, Point p2, Point k2)
{
	Point u = p1 - p2;
	double t = (k2 * u) / (k1 * k2);
	return {t * k1.x + p1.x, t * k1.y + p1.y};
}

Point get_point(const Line &l1, const Line &l2){return get_point_(l1.st, l1.ed - l1.st, l2.st, l2.ed - l2.st);}

bool on_right(const Line &a, const Line &b, const Line &c)
{
	Point t = get_point(b, c);
	return sgn((a.ed - a.st) * (t - a.st)) <= 0;
}

void half_plane(Line l[], int n)
{
	sort(l + 1, l + n + 1, lcmp);
	hh = 1;
	for (int i = 1; i <= n; ++ i)
	{
		if (i > 1 && !cmp(angle(l[i]), angle(l[i - 1]))) continue;
		while (hh < tt && on_right(l[i], l[q[tt - 1]], l[q[tt]])) tt --;
		while (hh < tt && on_right(l[i], l[q[hh]], l[q[hh + 1]])) hh ++;
		q[++ tt] = i;
	}
	while (hh < tt && on_right(l[q[hh]], l[q[tt - 1]], l[q[tt]])) tt --;
	while (hh < tt && on_right(l[q[tt]], l[q[hh]], l[q[hh + 1]])) hh ++;
	q[++ tt] = q[hh];
}

2）例题

T2：[JLOI 2013]赛车

题目传送门 Luogu

半平面交的模板题，注意用 long double，并且把精度调高一点。

详细代码 
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
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <vector>
#include <map>
using namespace std;

typedef pair<int, int> PII;
typedef long long LL;
typedef long double LD;
const LD eps = 1e-18;
const int N = 2e4 + 10;

struct Point{
	LD x, y;
	Point operator +(Point t)const{
		return {x + t.x, y + t.y};
	}
	Point operator -(Point t)const{
		return {x - t.x, y - t.y};
	}
	LD operator *(Point t)const{
		return x * t.y - y * t.x;
	}
};
struct Line{
	Point st, ed;
	vector<int> id;
}l[N];
int q[N], hh, tt;
int ki[N], vi[N];

auto angle = [](const Line &l){return atan2(l.ed.y - l.st.y, l.ed.x - l.st.x);};
auto sgn = [](LD x){return (fabs(x) < eps) ? 0 : (x > 0 ? 1 : -1);};
auto cmp = [](LD x, LD y){return sgn(x - y);};
auto lcmp = [](const Line &l1, const Line &l2){
	LD A = angle(l1), B = angle(l2);
	if (cmp(A, B) == 0) return sgn((l1.ed - l1.st) * (l2.ed - l1.st)) < 0;
	return A < B;
};

Point get_point_(Point p1, Point k1, Point p2, Point k2)
{
	Point del = p1 - p2;
	LD t = (k2 * del) / (k1 * k2);
	return {t * k1.x + p1.x, t * k1.y + p1.y};
}

Point get_point(const Line &l1, const Line &l2){return get_point_(l1.st, l1.ed - l1.st, l2.st, l2.ed - l2.st);}

bool on_right(const Line &a, const Line &b, const Line &c)
{
	Point t = get_point(b, c);
	return sgn((t - a.st) * (a.ed - a.st)) > 0;
}

void half_plane(int n)
{
	sort(l + 1, l + n + 1, lcmp);
	hh = 1, tt = 0;
	for (int i = 1; i <= n; ++ i)
	{
		if (i > 1 && !cmp(angle(l[i]), angle(l[i - 1]))) continue;
		while (hh < tt && on_right(l[i], l[q[tt - 1]], l[q[tt]])) tt --;
		while (hh < tt && on_right(l[i], l[q[hh]], l[q[hh + 1]])) hh ++;
		q[++ tt] = i;
	}
	while (hh < tt && on_right(l[q[tt]], l[q[hh]], l[q[hh + 1]])) hh ++;
	while (hh < tt && on_right(l[q[hh]], l[q[tt - 1]], l[q[tt]])) tt --;
	q[++ tt] = q[hh];
	vector<int> ans;
	for (int i = hh; i < tt; ++ i)
		for (auto x : l[q[i]].id) ans.push_back(x);
	sort(ans.begin(), ans.end());
	cout << ans.size() << endl;
	for (auto x : ans) cout << x << ' ';
}

int main()
{
	map<PII, vector<int> > ids;
	int cnt = 0, n;
	cin >> n;
	l[++ cnt] = {{0, 0}, {10000, 0}, vector<int>(0)};
	l[++ cnt] = {{0, 10000}, {0, 0}, vector<int>(0)};
	for (int i = 1; i <= n; ++ i) cin >> ki[i];
	for (int i = 1; i <= n; ++ i) cin >> vi[i];
	for (int i = 1; i <= n; ++ i)
		ids[{ki[i], vi[i]}].push_back(i);
	for (const auto &c : ids)
		l[++ cnt] = {{0, c.first.first}, {1, c.first.first + c.first.second}, c.second};
	half_plane(cnt);
	return 0;
}

3. 三维凸包

1）三维空间向量

a. 加减、数乘、模长

与二维向量相同，不再赘述。

b. 点乘

得到是一个数。

$(x1, y1, z1) \cdot (x2, y2, z2) = x1x2 + y1y2 + z1z2$。

注意满足 $|A||B|\cos = A \cdot B$。

c. 叉乘

得到不是数，是一个行列式的结果。

注意 $i, j, k$ 是空间单位向量，得到的明显也是一个三位向量。

展开行列式，可以得到：$(x1, y1, z1) \times (x2, y2, z2) = (y1z2 - y2z1, x2z1 - x1z2, x1y2 - x2y1)$。

d. 多面体欧拉定理

点数 - 棱数 + 面数 = 2。

e. 平面的法向量

法向量是指垂直于

任意取两个向量

2）三维凸包 - 增量法

其实是一个暴力算法，时间复杂度为 $O(n ^ 2)$。

我们首先找到任意不共面的 4 个点，使得可以形成一个凸包。

然后，我们尝试加入每一个点。

如果它在凸包里，则跳过；否则我们考虑加入。

假设这个点是一个太阳，那么一定会有一些面是“白天”，还有一些棱是“晨昏线”。

我们将”白天“全部删除，并且将”晨昏线“与新加的点相连，形成新的凸包。

具体看代码。

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
int Convex_3d(int n)
{
	static bitset<N> g[N];
	int cnt = 0, ncnt = 0;
	pl[++ cnt] = {1, 2, 3}, pl[++ cnt] = {3, 2, 1};
	for (int i = 4; i <= n; ++ i)
	{
		ncnt = 0;
		for (int j = 1; j <= cnt; ++ j){
			bool t = pl[j].above(p[i]);
			if (!t) np[++ ncnt] = pl[j];
			for (int k = 0; k < 3; ++ k)
				g[pl[j].v[k]][pl[j].v[(k + 1) % 3]] = t;
		}
		for (int j = 1; j <= cnt; ++ j)
			for (int k = 0; k < 3; ++ k)
			{
				int a = pl[j].v[k], b = pl[j].v[(k + 1) % 3];
				if (g[a][b] && !g[b][a]) np[++ ncnt] = {a, b, i};
			}
		cnt = ncnt;
		for (int j = 1; j <= cnt; ++ j) pl[j] = np[j];
	}
	return cnt;
}

4. 旋转卡壳

其实不是一种模板或者算法，而是一种思想或者是做题的技巧。

1）思想

我们定义对踵点为任意两条平行的直线，向中间靠拢的时候碰到的点。

很明显，我们的对踵点的个数是 $O(n)$ 级别的，所以我们并不需要枚举每一条平行的直线。我们可以考虑直接枚举与凸包边平行的边。

下面假设我们要寻找直径，也就是最长的线段。

我们先确定 $p(i)$ 为对踵点的一端，我们怎样才能寻找到与之相对的对踵点呢？

很明显，我们可以暴力枚举，但时间复杂度是 $O(n ^ 2)$，不优。

我们可以发现，随着我们的 $p(i)$ 一直都是顺时针（或者逆时针）旋转的，$p(j)$ 也一定是顺时针旋转的。

所以，我们可以使用双指针算法，可以做到 $O(n)$，瓶颈在前面的凸包，时间复杂度为 $O(n \log n)$。

怎样判断哪个点是最远的呢？我们回顾叉积的定义，发现是三角形的面积。我们再拉一个 $p(i + 1)$ 过来，就可以求面积了，同时这条直线就是 $p(i), p(i + 1)$。

具体看代码。

1
2
3
4
5
for (int i = 0, j = 2; i < n; ++ i)
    {
        while ((s[i + 1] - s[i]) * (s[j] - s[i]) < (s[i + 1] - s[i]) * (s[j + 1] - s[i])) j = (j + 1) % n;
        ans = max(ans, max(dist2(s[i], s[j]), dist2(s[i + 1], s[j])));
    }

2）例题

T1：求最远点距离

题目传送门 Luogu

题目传送门 AcWing

刚才讲到的，直接上代码。

详细代码 
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
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cstdio>
using namespace std;

typedef long long LL;

const int N = 5e4 + 10;
struct Point{
    int x, y;
    Point operator -(Point t)const{
        return {x - t.x, y - t.y};
    }
    LL operator *(Point t)const{
        return (LL)x * t.y - (LL)y * t.x;
    }
}p[N];
int stk[N], top, n;

auto cmp = [](Point t1, Point t2){
    return t1.x == t2.x ? t1.y < t2.y : t1.x < t2.x;
};

LL dist2(Point p1, Point p2){
    return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y);
}

void Andrew()
{
    static bool usd[N];
    sort(p + 1, p + n + 1, cmp);
    for (int i = 1; i <= n; ++ i)
    {
        while (top >= 2 && (p[stk[top]] - p[stk[top - 1]]) * (p[i] - p[stk[top]]) <= 0)
            usd[stk[top --]] = false;
        usd[stk[++ top] = i] = true;
    }
    usd[1] = false;
    for (int i = n; i; -- i)
    {
        if (usd[i]) continue;
        while (top >= 2 && (p[stk[top]] - p[stk[top - 1]]) * (p[i] - p[stk[top]]) < 0)
            usd[stk[top --]] = false;
        usd[stk[++ top] = i] = true;
    }
    top --;
}

LL get_dist()
{
    if (top <= 2) return dist2(p[stk[1]], p[stk[2]]);
    LL ans = 0;
    for (int i = 1, j = 2; i <= top; ++ i)
    {
        while ((p[stk[i + 1]] - p[stk[i]]) * (p[stk[j]] - p[stk[i]]) < (p[stk[i + 1]] - p[stk[i]]) * (p[stk[j + 1]] - p[stk[i]]))
            j = j % top + 1;
        ans = max(ans, max(dist2(p[stk[i + 1]], p[stk[j]]), dist2(p[stk[i]], p[stk[j]])));
    }
    return ans;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++ i) cin >> p[i].x >> p[i].y;
    Andrew();
    cout << get_dist() << endl;
    return 0;
}

T2：最小矩形覆盖

题目传送门 Luogu

题目传送门 AcWing

先搁着 qwq

详细代码 
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
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;

const int N = 3e5 + 10;
const double INF = 1e20, eps = 1e-12, Pi = acos(-1);

struct Point{
    double x, y;
    Point operator +(Point t)const{
        return {x + t.x, y + t.y};
    }
    Point operator -(Point t)const{
        return {x - t.x, y - t.y};
    }
    Point operator *(double t)const{
        return {x * t, y * t};
    }
    Point operator /(double t)const{
        return {x / t, y / t};
    }
    double operator *(Point t)const{
        return x * t.y - y * t.x;
    }
    double operator &(Point t)const{
        return x * t.x + y * t.y;
    }
}p[N], ans[4];
int stk[N], top, n;
double mxa = INF;

int sgn(double x){
    return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1);
}
int cmp(double x, double y){return sgn(x - y);}
double len(Point a){return sqrt(a & a);}
double area(Point a, Point b, Point c){return (b - a) * (c - a);}
Point normal(Point t){return t / len(t);}
Point rotate(Point a, double th){return {a.x * cos(th) + a.y * sin(th), a.y * cos(th) - a.x * sin(th)};}
bool pcmp(Point a, Point b)
{
    return a.x == b.x ? a.y < b.y : a.x < b.x;
}

double project(Point a, Point b, Point c)
{
    return ((b - a) & (c - a)) / len(b - a);
}

void Andrew()
{
    static bool usd[N];
    sort(p + 1, p + n + 1, pcmp);
    for (int i = 1; i <= n; ++ i)
    {
        while (top >= 2 && area(p[stk[top - 1]], p[stk[top]], p[i]) >= 0)
            usd[stk[top --]] = false;
        usd[stk[++ top] = i] = true;
    }
    usd[1] = false;
    for (int i = n; i; -- i)
    {
        if (usd[i]) continue;
        while (top >= 2 && area(p[stk[top - 1]], p[stk[top]], p[i]) >= 0)
            usd[stk[top --]] = false;
        usd[stk[++ top] = i] = true;
    }
    reverse(stk + 1, stk + top + 1);
    top --;
}

void find_ju()
{
    for (int i = 1, j = 3, k = 2, l = 3; i <= top; ++ i)
    {
        Point d = p[stk[i]], e = p[stk[i + 1]];
        while (cmp(area(d, e, p[stk[j]]), area(d, e, p[stk[j + 1]])) < 0)
            j = j % top + 1;
        while (cmp(project(d, e, p[stk[k]]), project(d, e, p[stk[k + 1]])) < 0)
            k = k % top + 1;
        if (i == 1) l = j;
        while (cmp(project(d, e, p[stk[l]]), project(d, e, p[stk[l + 1]])) > 0)
            l = l % top + 1;
        Point x = p[stk[j]], y = p[stk[k]], z = p[stk[l]];
        double h = area(d, e, x) / len(e - d),
               w = ((y - z) & (e - d)) / len(e - d);
        // cout << i << ' ' << h << ' ' << w << endl;
        if (h * w < mxa){
            mxa = h * w;
            ans[0] = d + normal(e - d) * project(d, e, y);
            ans[3] = d + normal(e - d) * project(d, e, z);
            Point t = normal(rotate(e - d, -Pi / 2));
            ans[1] = ans[0] + t * h;
            ans[2] = ans[3] + t * h;
        }
    }
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++ i) cin >> p[i].x >> p[i].y;
    Andrew();
    // cout << top << endl, exit(0);
    find_ju();
    printf("%.5lf\n", mxa);
    int k = 0;
    for (int i = 1; i < 4; ++ i)
        if (cmp(ans[i].y, ans[k].y) < 0 || !cmp(ans[i].y, ans[k].y) && cmp(ans[i].x, ans[k].x) < 0)
            k = i;
    for (int i = 0; i < 4; ++ i, k ++)
    {
        if (k == 4) k = 0;
        double x = ans[k].x + eps, y = ans[k].y + eps;
        if (fabs(x) < eps) x = 0;
        if (fabs(y) < eps) y = 0;
        printf("%.5lf %.5lf\n", x, y);
    }
    return 0;
}