数几何

半平面交

字数统计: 2k字   |   阅读时长: 8分
,

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

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的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