Dyd's Blog

He who has a strong enough why can bear almost any how.

凸包与半平面交

Posted on In
Symbols count in article: 8.4k Reading time ≈ 8 mins.

相似的算法

凸包与半平面交

问题引入

给定共面的 $n$ 个点，求一个周长最小的凸多边形，使得所有点都在这个多边形边上和内部，这个凸多边形被称为凸包

Andrew算法

解决这个问题有多个算法，这里介绍Andrew算法，其流程如下：

  1. 将所有点按 $x$ 为第一关键字排序
  2. 先考虑上凸包，具体地：
    • 维护一个栈，前两个点无条件入栈（第一个点定在凸包上）
    • 考虑栈顶元素和它的上一个点所在直线，如果当前点在直线上方，就弹出栈顶，加入当前点
    • 当前点变成下一个点
  3. 然后再反过来做一遍，维护下凸包，只需将上面的“直线上方”改成“直线下方”即可，注意一定要再扫一次第一个点
  4. 最后，注意特判所有点共线的情况

时间瓶颈在排序上，为 $O(n \log n)$

例题一

板子题

围住奶牛

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
#include <bits/stdc++.h>
#define DB double
#define STC static
using namespace std;
const int N = 1e5 + 5;

int n;
struct Point
{
    DB x, y;
    Point operator - (const Point &t) const
    {
        return {x - t.x, y - t.y};
    }
    bool operator < (const Point &t) const
    {
        return x == t.x ? y < t.y : x < t.x;
    }
} ;
DB dist(Point x, Point y)
{
    Point t = x - y;
    return sqrt(t.x * t.x + t.y * t.y);
}
DB cross(Point x, Point y)
{
    return x.x * y.y - y.x * x.y;
}
DB area(Point x, Point y, Point z)
{
    return cross(y - x, z - x);
}
DB andrew(Point x[])
{
    STC int stk[N];
    STC bool used[N];
    sort(x, x + n);
    int top = 0;
    for (int i = 0, t; i < n; ++i)
    {
        while (top >= 2 && area(x[stk[top - 1]], x[stk[top]], x[i]) <= 0)
        {
            // 凸包边界上的点即使被从栈中删掉，也不能删掉used上的标记
            if (!area(x[stk[top - 1]], x[stk[top]], x[i]))
                --top;
            else
                used[stk[top--]] = false;
        }
        stk[++top] = i;
        used[i] = true;
    }
    used[0] = false;
    for (int i = n - 1; i >= 0; --i)
    {
        if (used[i])
            continue;
        //为了应对所有点共线的情况,这里不能取<=
        while (top >= 2 && area(x[stk[top - 1]], x[stk[top]], x[i]) < 0)
            top -- ;
        stk[++top] = i;
    }
    DB res = 0;
    for (int i = 2; i <= top; ++i)
        res += dist(x[stk[i - 1]], x[stk[i]]);
    return res;
}
int main()
{
    STC Point q[N];
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ )
        scanf("%lf %lf", &q[i].x, &q[i].y);
    printf("%.2lf\n", andrew(q));
    return 0;
}

例题二

几乎板子题……

信用卡凸包

首先，如果信用卡是矩形，显然把所有顶点记下来求凸包即可

而对于四个角是圆的情况，每一个信用卡可以看成四个圆，我们现在要对这些圆求凸包

把凸包分成两部分，一部分是直线，一部分是曲线

曲线部分恰好是一个圆（凸多边形外交和为 $360$ ），而直线部分等于是求圆心的凸包

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
#include <bits/stdc++.h>
#define DB double
#define STC static
using namespace std;
const int N = 1e5 + 5;
const DB pi = acos(-1);
int n;
int dx[] = {1, 1, -1, -1}, dy[] = {1, -1, -1, 1};
struct Point
{
    DB x, y;
    Point operator - (const Point &t) const
    {
        return {x - t.x, y - t.y};
    }
    bool operator < (const Point &t) const
    {
        return x == t.x ? y < t.y : x < t.x;
    }
} ;
DB dist(Point x, Point y)
{
    Point t = x - y;
    return sqrt(t.x * t.x + t.y * t.y);
}
DB cross(Point x, Point y)
{
    return x.x * y.y - y.x * x.y;
}
DB area(Point x, Point y, Point z)
{
    return cross(y - x, z - x);
}
Point rotate(Point x, DB th)
{
    return {x.x * cos(th) + x.y * sin(th), -x.x * sin(th) + x.y * cos(th)};
}
DB andrew(Point x[], int len)
{
    STC int stk[N];
    STC bool used[N];
    sort(x, x + len);
    int top = 0;
    for (int i = 0, t; i < len; ++i)
    {
        while (top >= 2 && area(x[stk[top - 1]], x[stk[top]], x[i]) <= 0)
        {
            // 凸包边界上的点即使被从栈中删掉，也不能删掉used上的标记
            if (!area(x[stk[top - 1]], x[stk[top]], x[i]))
                --top;
            else
                used[stk[top--]] = false;
        }
        stk[++top] = i;
        used[i] = true;
    }
    used[0] = false;
    for (int i = len - 1; i >= 0; --i)
    {
        if (used[i])
            continue;
        //为了应对所有点共线的情况,这里不能取<=
        while (top >= 2 && area(x[stk[top - 1]], x[stk[top]], x[i]) < 0)
            top -- ;
        stk[++top] = i;
    }
    DB res = 0;
    for (int i = 2; i <= top; ++i)
        res += dist(x[stk[i - 1]], x[stk[i]]);
    return res;
}
int main()
{
    STC Point q[N];
    DB a, b, r;
    scanf("%d %lf %lf %lf", &n, &a, &b, &r);
    a = a / 2 - r, b = b / 2 - r;
    int cnt = 0;
    while (n--)
    {
        DB x, y, z;
        scanf("%lf %lf %lf", &x, &y, &z);
        for (int i = 0; i < 4; ++i)
        {
            Point t = rotate({dx[i] * b, dy[i] * a}, -z);
            q[cnt++] = {x + t.x, y + t.y};
        }
    }
    printf("%.2lf\n", andrew(q, cnt) + 2 * pi * r);
    return 0;
}

问题又一次引入

众所周知，一条直线将一个面划分成两半

给定共面的 $n$ 条直线，每条直线都将平面划分成两半，我们只取其中一半，将取的所有部分求交集，求最后的面积，这就叫半平面交

solve

为了简化问题，我们用向量表示直线（给直线赋予方向），固定我们取向量的左边，这样，对应两条角度相同的向量（就是直线），我们只保留左边那条，剩下的直线一定不存在角度相同的向量（即共向向量），然后：

  1. 将向量按照角度排序（这里使用 atan2(y, x) 函数，注意是 (y, x) 不是 (x, y) ）
  2. 用双端队列维护，按顺序扫描每一个向量，如果队头/队尾向量的交点在当前向量的右边，就删掉
  3. 最后用队头更新一下队尾，用队尾跟新一下队头

时间 $O(n \log n)$

例题一

板子

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
#include <bits/stdc++.h>
#define DB double
#define STC static
using namespace std;
const int N = 500 + 5;
const DB eps = 1e-8;
struct Point
{
	DB x, y;
	Point operator - (const Point &t) const
	{
		return {x - t.x, y - t.y};
	}
	
} ans[N];
struct Line
{
	Point st, ed;
	DB angle()
	{
		return atan2(ed.y - st.y, ed.x - st.x);
	}
} ;

int cmp(DB x, DB y)
{
	if (fabs(x - y) < eps)
		return 0;
	return x > y ? 1 : -1;
}
DB cross(Point x, Point y)
{
	return x.x * y.y - y.x * x.y;
}
DB area(Point x, Point y, Point z)
{
	return cross(y - x, z - x);
}
bool operator < (Line &x, Line &y)
{
	DB a = x.angle(), b = y.angle();
	if (!cmp(a, b))
		return area(x.st, x.ed, y.ed) < 0;
	return a < b;
}
Point jiao(Point p, Point u, Point q, Point v)
{
	DB t = cross(v, p - q) / cross(u, v);
	return {p.x + u.x * t, p.y + u.y * t};
}
Point jiao(Line x, Line y)
{
	return jiao(x.st, x.ed - x.st, y.st, y.ed - y.st);
}
bool on_r(Line &x, Line y, Line z)
{
	Point o = jiao(y, z);
	return cmp(area(x.st, x.ed, o), 0) <= 0; //这里取不取等依题意
}
DB halfp(Line x[], int len)
{
	sort(x, x + len);
	int hh = 0, tt = -1;
	STC int q[N];
	for (int i = 0; i < len; ++i)
	{
		if (i && !cmp(x[i].angle(), x[i - 1].angle()))
			continue;
		while (hh < tt && on_r(x[i], x[q[tt - 1]], x[q[tt]]))
			--tt;
		while (hh < tt && on_r(x[i], x[q[hh + 1]], x[q[hh]]))
			++hh;
		q[++tt] = i;
	}
	while (hh < tt && on_r(x[q[hh]], x[q[tt - 1]], x[q[tt]]))
		--tt;
	while (hh < tt && on_r(x[q[tt]], x[q[hh + 1]], x[q[hh]]))
		++hh;
	q[++tt] = q[hh];
	int k = 0;
    for (int i = hh; i < tt; ++i)
        ans[k++] = jiao(x[q[i]], x[q[i + 1]]);
   	DB res = 0;
    for (int i = 1; i + 1 < k; ++i)
        res += area(ans[0], ans[i], ans[i + 1]);
    return res / 2;
}
int main()
{
	int n, m, cnt = 0;
	STC Line l[N];
	STC Point pg[N];
	scanf("%d", &n);
	while (n--)
	{
		scanf("%d", &m);
		for (int i = 0; i < m; ++i)
			scanf("%lf %lf", &pg[i].x, &pg[i].y);
		for (int i = 0; i < m; ++i)
			l[cnt++] = {pg[i], pg[(i + 1) % m]};
	}
	printf("%.3lf\n", halfp(l, cnt));
	return 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
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
#include <bits/stdc++.h>
#define STC static
#define LD long double
#define PII pair<int, int>
using namespace std;
const int N = 10000 + 5;
const LD eps = 1e-18;
struct Point
{
	LD x, y;
	Point operator - (const Point &t) const
	{
		return {x - t.x, y - t.y};
	}
	
} ;
int ans[N];
struct Line
{
	Point st, ed;
	vector<int> id;
	LD angle()
	{
		return atan2(ed.y - st.y, ed.x - st.x);
	}
} ;
int cmp(LD x, LD y)
{
	if (fabs(x - y) < eps)
		return 0;
	return x > y ? 1 : -1;
}
LD cross(Point x, Point y)
{
	return x.x * y.y - y.x * x.y;
}
LD area(Point x, Point y, Point z)
{
	return cross(y - x, z - x);
}
bool operator < (Line &x, Line &y)
{
	LD a = x.angle(), b = y.angle();
	if (!cmp(a, b))
		return area(x.st, x.ed, y.ed) < 0;
	return a < b;
}
Point jiao(Point p, Point u, Point q, Point v)
{
	LD t = cross(v, p - q) / cross(u, v);
	return {p.x + u.x * t, p.y + u.y * t};
}
Point jiao(Line x, Line y)
{
	return jiao(x.st, x.ed - x.st, y.st, y.ed - y.st);
}
bool on_r(Line &x, Line y, Line z)
{
	Point o = jiao(y, z);
	return cmp(area(x.st, x.ed, o), 0) < 0;
}
void halfp(Line x[], int len)
{
	sort(x, x + len);
	int hh = 0, tt = -1;
	STC int q[N];
	for (int i = 0; i < len; ++i)
	{
		if (i && !cmp(x[i].angle(), x[i - 1].angle()))
			continue;
		while (hh < tt && on_r(x[i], x[q[tt - 1]], x[q[tt]]))
			--tt;
		while (hh < tt && on_r(x[i], x[q[hh + 1]], x[q[hh]]))
			++hh;
		q[++tt] = i;
	}
	while (hh < tt && on_r(x[q[hh]], x[q[tt - 1]], x[q[tt]]))
		--tt;
	while (hh < tt && on_r(x[q[tt]], x[q[hh + 1]], x[q[hh]]))
		++hh;
	q[++tt] = q[hh];
	int k = 0;
    for (int i = hh; i <= tt; ++i)
		for (auto id : x[q[i]].id)
        	ans[k++] = id;
   	sort(ans, ans + k);
	printf("%d\n", k);
    for (int i = 0; i < k; ++i)
        printf("%d ", ans[i]);
}
int main()
{
	map<PII, vector<int> > ids;
	int n, cnt = 0;
	STC Line l[N];
	STC int ki[N], vi[N];
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i)
		scanf("%d", &ki[i]);
	for (int i = 1; i <= n; ++i)
		scanf("%d", &vi[i]);
	for (int i = 1; i <= n; ++i)
		ids[{ki[i], vi[i]}].push_back(i);
	l[cnt++] = {{0, N}, {0, 0}};
	l[cnt++] = {{0, 0}, {N, 0}};
	for (auto i : ids) //简直毒瘤的map遍历
        l[cnt++] = {{0, i.first.first}, {1, i.first.first + i.first.second}, i.second};
	halfp(l, cnt);
	return 0;
}