半平面交

  |  

摘要: 本文介绍计算几何中半平面交算法

【对算法,数学,计算机感兴趣的同学,欢迎关注我哈,阅读更多原创文章】
我的网站:潮汐朝夕的生活实验室
我的公众号:算法题刷刷
我的知乎:潮汐朝夕
我的github:FennelDumplings
我的leetcode:FennelDumplings


各位好,在文章 向量的实现 中,我们实现了一个二维向量的代码模板,主要是 Vector2 这个类。

此后基于二维向量的代码模板,我们解决了计算几何中关于直线和线段极角多边形中的常见问题,以及力扣中的几个相关题目。

在文章 凸包凸包的直径,旋转卡尺算法 中,我们学习了凸包的算法,本文进一步看一下怎么用旋转卡尺算法求凸包的直径。

本文我们学习半平面交算法,主要的应用是求多边形的核、求可以放进多边形的圆的最大半径等。首先复习一下所需的有向直线的代码模板,然后介绍半平面交的定义,以及用半平面交求多边形的核的算法和代码模板。


有向直线

这里简要复习一下有向直线的代码模板,算法细节参考文章 直线和线段

点和方向向量的形式 $\vec{a} + p \cdot \vec{b}$ 的有向直线的代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
struct Line
{
// line: a+pb, a 是点, b 是方向向量
Vector2 a;
Vector2 b;
double deg; // 极角
Line(Vector2 a, Vector2 b):a(a),b(b)
{
deg = b.polar();
}
bool operator<(const Line& l) const
{
// 方向向量的极角序
return deg < l.deg;
}
};

判断一个点是否在有向直线的左侧,代码如下:

1
2
3
4
5
bool on_left(const Line& l, const Vector2& p)
{
// 点 p 是否在有向直线 l 左侧
return l.b.cross(p - l.a) > EPS;
}

求两个有向直线的交点。需要注意平行的情况,代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
// 直线与直线相交
bool line_intersection(const Line& l1, const Line& l2, Vector2& x)
{
double det = l1.b.cross(l2.b); // 两个方向向量的外积
// 平行或重叠返回false
if(fabs(det) < EPS)
return false;
// 两条有向直线的交点
Vector2 u = l1.a - l2.a;
double t = l2.b.cross(u) / det;
x = l1.a + l1.b * t;
return true;
}

半平面交

定义

  • 半平面:一条直线把二维平面分成两个平面

  • 半平面交:定义一个半平面为一个向量的左侧,半平面交即为若干个半平面的交集。

半平面交的结果

肯定是凸的:

  1. 凸多边形
  2. 无界,因为有可能若干半平面没有形成封闭
  3. 直线,线段,点,空

半平面交的应用

  1. 求解一个区域,可以看到给定图形的各个角落。(多边形的核)
  2. 求可以放进多边形的圆的最大半径。

算法 (求多边形的核)

选取一个正方向,一般是逆时针。核在有向直线的左边。

多边形 vector<Vector2> p,p 中为多边形的点,逆时针给出。以顶点 p[i] 为起点,逆时针方向的方向向量为 (p[(i + 1) % n] - p[i]).normalize()

1
2
3
Line lines(n);

lines[i] = Line(p[i], (p[(i + 1) % n] - p[i]).normalize());

将有向直线按极角排序 [0, 2PI),其中有向直线用点和方向向量的形式: $\vec{a} + p \cdot \vec{b}$,其中 $\vec{a}$ 为 p[i],$\vec{b}$ 为 (p[(i + 1) % n] - p[i]).normalize()

按顺序遍历有向直线,取左边区域,删右边区域,最后剩下的区域为半平面交。

半平面交通过有向直线下标集合 L的方式维护

最终留下的直线集合的直线数目如果至少有三条,则有核,核为全部直线围成的凸包。

求这些直线围成的凸包面积:先按顺序求出交点,再用叉积求这些交点围成的多边形面积。

  • 取左边区域,删右边区域的具体做法
1
2
3
4
5
6
7
8
1. 以逆时针为正方向,建边。当输入方向不确定时,用叉乘面积的正负判定输入的方向。
2. 对线段根据极角排序
3. 去除极角相同的情况下,位置在右边的边
4. 用双端队列储存线段集合 deq,遍历所有线段
5. 判断该线段加入后对半平面交的影响,对双端队列的头部和尾部进行判断,因为线段加入是有序的
如果某条线段对于新的半平面交没有影响,则从队列中剔除掉
没有影响:队列头部(尾部)两条直线的交点在新加入直线的右侧
6. 最后剩下的线段集合 L,即使最后要求的半平面交。

代码模板

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
vector<Vector2> half_plane_intersection(vector<Line>& lines)
{
// 极角序
sort(lines.begin(), lines.end());
// 去除极角相同时右边的边
vector<Line> tmp;
for(int i = 0; i < (int)lines.size(); ++i)
{
if(tmp.empty() || fabs(tmp.back().deg - lines[i].deg) > EPS)
{
tmp.push_back(lines[i]);
}
else if(on_left(tmp.back(), lines[i].a))
{
tmp.pop_back();
tmp.push_back(lines[i]);
}
}
lines.swap(tmp);
int N = lines.size();
// cout << N << endl;
// 有向直线集合表示半平面,初始只有半平面 {0}
vector<Line> deq(N); // 数组形式双端队列
// deq_p[i] := lines[i] 与 lines[(i + 1) % N] 的交点
vector<Vector2> deq_p(N);
int left = 0, right = 0;
deq[0] = lines[0];
for(int i = 1; i < N; ++i)
{
// cout << lines[i].a.x << " " << lines[i].a.y << " : " << lines[i].deg << endl;
// 判断该线段加入后对半平面交的影响,对双端队列的头部和尾部进行判断,因为线段加入是有序的
while(left < right && !on_left(lines[i], deq_p[right - 1]))
--right;
while(left < right && !on_left(lines[i], deq_p[left]))
++left;
deq[++right] = lines[i];
if(left < right)
line_intersection(deq[right - 1], deq[right], deq_p[right - 1]);
}
while(left < right && !on_left(deq[left], deq_p[right - 1]))
--right;
if(right - left <= 1)
return {};
// 首位两个半平面的交
line_intersection(deq[right], deq[left], deq_p[right]);
vector<Vector2> poly(right - left + 1);
for(int i = 0; i <= right - left; ++i)
poly[i] = deq_p[left + i];
return poly;
}

题目

以下题目的代码中的用到的代码模板如下:

  • Vector2
  • Line
  • on_left
  • line_intersection
  • half_plane_intersection
  • area

P2283 多边形

给定一个多边形(不一定是凸多边形)展厅,展厅中有些区域可以看到展厅内的所有地方,求这些区域的面积。

代码 (C++)

这是洛谷的一道题,比较直接,就是求多边形的核的面积。

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
int main()
{
int n;
cin >> n;
vector<Vector2> A; // 点
vector<Line> lines; // 边
A.assign(n, Vector2());
for(int j = 0; j < n; ++j)
{
int x, y;
cin >> x >> y;
A[j].x = x;
A[j].y = y;
}
// 如果有向面积为负,则点是顺时针给的,翻转成逆时针
if(area(A) < -EPS)
reverse(A.begin(), A.end());
// 建边,逆时针方向
for(int j = 0; j < n; ++j)
{
int k = (j + 1) % n;
lines.emplace_back(A[j], (A[k] - A[j]).normalize());
}
vector<Vector2> poly = half_plane_intersection(lines);
double ans = fabs(area(poly));
cout << std::fixed << std::setprecision(2);
cout << ans + 1e-9 << endl;
}

2803. 凸多边形

逆时针给出 n 个凸多边形的顶点坐标,求它们交的面积。

代码 (C++)

这是 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
int main()
{
int n;
cin >> n;
vector<Vector2> A; // 点
vector<Line> lines; // 边
for(int i = 0; i < n; ++i)
{
int m;
cin >> m;
A.assign(m, Vector2());
for(int j = 0; j < m; ++j)
{
int x, y;
cin >> x >> y;
A[j].x = x;
A[j].y = y;
}
// 如果有向面积为负,则点是顺时针给的,翻转成逆时针
if(area(A) < -EPS)
reverse(A.begin(), A.end());
// 建边,逆时针方向
for(int j = 0; j < m; ++j)
{
int k = (j + 1) % m;
lines.emplace_back(A[j], (A[k] - A[j]).normalize());
}
}
vector<Vector2> poly = half_plane_intersection(lines);
double ans = area(poly);
cout << std::fixed << std::setprecision(3);
cout << ans + 1e-9 << endl;
}

Share