半平面交笔记 Half Plane Intersection

发表于2020-08-08更新于2023-09-15字数统计4.1k阅读时长40分

Half Plane Intersection

参考资料

半平面

一条直线和直线的一侧。半平面是一个点集，因此是一条直线和直线的一侧构成的点集。当包含直线时，称为闭半平面；当不包含直线时，称为开半平面。
解析式一般为 Ax+By+C≥0Ax + By + C \ge 0。
在计算几何中用向量表示，整个题统一以向量的左侧或右侧为半平面。

半平面交

半平面交是指多个半平面的交集。因为半平面是点集，所以点集的交集仍然是点集。在平面直角坐标系围成一个区域(可以理解为向量集中每一个向量的左侧交)。
这就很像普通的线性规划问题了，得到的半平面交就是线性规划中的可行域。一般情况下半平面交是有限的，经常考察面积等问题的解决。

多边形的核

如果一个点集中的点与多边形上任意一点的连线与多边形没有其他交点，那么这个点集被称为多边形的核。
把多边形的每条边看成是首尾相连的向量，那么这些向量在多边形内部方向的半平面交就是多边形的核。

解法 - S&I算法

极角排序

利用 atan2(y, x) 函数将向量(线段)按极角排序，排序时，若遇到贡献向量(且方向相同)，则取靠近可行域的一个(即靠近左侧的向量)。判断方法为取一个向量的起点与另一个向量的起点、终点构成向量，通过叉积符号判断。

极角排序
1
2
3
4
5
6
7
8
9
10
11
/*cross product of two Vec2 vector (starting point is zero)*/
friend inline double operator*(const Vec2& a, const Vec2& b) {
    return a.x * b.y - a.y * b.x;
}
//s是Seg类的起点，t为终点
friend inline bool operator<(const Seg& a, const Seg& b) {//比较极角
    // theta1 != theta2 按theta排序
    if (abs(a.theta - b.theta) > EPS) return a.theta < b.theta;
    // theta1 == theta2 判断向量a在b的哪边，令最靠左的排在最左边
    return (b.s - a.s) * (b.t - a.s) > EPS;
}

维护单调队列

因为半平面交是一个凸多边形，所以需要维护一个凸壳。因为后来加入边的只可能会影响最开始加入的或最后加入的边（此时凸壳连通），只需要删除队首和队尾的元素，所以需要用单调队列。
具体维护方法见 oi-wiki.org

得到半平面交

如果半平面交是一个凸 nn 边形，最后在交点数组里会得到 nn 个点。我们再把它们首尾相连，就是一个统一方向（顺或逆时针）的多边形。
此时就可以用三角剖分求面积了。（求面积是最基础的考法）
偶尔会出现半平面交不存在或面积为 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
struct Vec2 {  //点向量类
    Vec2() = default;
    Vec2(double _x, double _y) : x(_x), y(_y) {}
    /*Vec2 vector a - Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }
    /*Vec2 vector a + Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*cross product of two Vec2 vector (starting point is zero)*/
    friend inline double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;
    }
    /*Vec2 * num operation @return Vec2(x, y)*/
    friend inline Vec2 operator*(const Vec2& a, double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*input coordinate of Vec2 through istream*/
    inline friend istream& operator>>(istream& is, Vec2& p) {
        is >> p.x >> p.y;
        return is;
    }
    /*@param x x坐标 @param y y坐标*/
    double x = 0, y = 0;
} vertex[MAXN], its[MAXN]/*intersected points*/;

线段类(极角排序，对应直线的交点)

线段类
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
struct Seg {  //线段类 起点s, 终点t
    Seg() = default;
    Seg(Vec2 _s, Vec2 _t) : s(_s), t(_t) { theta = atan2((t - s).y, (t - s).x); }
    /*comparison of two Seg @return bool*/
    friend inline bool operator<(const Seg& a, const Seg& b) {//比较极角
        // theta1 != theta2 按theta排序
        if (abs(a.theta - b.theta) > EPS) return a.theta < b.theta;
        // theta1 == theta2 判断向量a在b的哪边，令最靠左的排在最左边
        return (b.s - a.s) * (b.t - a.s) > EPS;
    }
    /*intersection of two Seg @return Vec2(x, y)*/
    inline friend Vec2 intersect(const Seg& a, const Seg& b) {
        double ratio = ((b.t - b.s) * (a.s - b.s)) / ((a.t - a.s) * (b.t - b.s));
        return a.s + (a.t - a.s) * ratio;
    }
    /*@param s 起点 @param t 终点*/
    Vec2 s, t;
    /*@param theta 极角*/
    double theta = 1e9;
} seg[MAXN], q[MAXN]/*deque 队列中存放线段*/;

半平面交

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
//求向量左侧的半平面交 
//队列中线段范围 (l, r]，若半平面交不闭合，则范围为 (l, r)
//第一个交点是 its[l + 2]，因为 seg[l + 2] 为第二条线段
inline void halfPlaneIntersection() {
    sort(seg + 1, seg + tot + 1);
    int l = 0, r = 0;
    for (int i = 1; i <= tot; i++) {
        if (abs(seg[i].theta - seg[i - 1].theta) > EPS) {
            while (r - l > 1 && (seg[i].t - its[r]) * (seg[i].s - its[r]) > EPS)
                --r;
            while (r - l > 1 && (seg[i].t - its[l + 2]) * (seg[i].s - its[l + 2]) > EPS)
                ++l;
            q[++r] = seg[i];
            if (r - l > 1) its[r] = intersect(q[r], q[r - 1]); //求新交点
        }
    }
    while (r - l > 1 && (q[l + 1].t - its[r]) * (q[l + 1].s - its[r]) > EPS)
        --r; //删除多余元素
    its[r + 1] = intersect(q[l + 1], q[r]), ++r; //终点线段与起点线段相交

    double ans(0);
    for (int i = l + 2; i < r; i++) { // l + 2 是起始交点 三角剖分求面积
        ans += (its[i] * its[i + 1]) / 2;
    }
    ans += (its[r] * its[l + 2]) / 2;
    cout << setprecision(3) << fixed << ans << endl;
}

例题

A CQOI2006 凸多边形

Description

模板题

Code

CQOI2006 凸多边形
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
#include <bits/stdc++.h>
#define PI acos(-1)

using namespace std;

constexpr double EPS = 1e-8;
constexpr int MAXN = 1e5 + 5;

struct Vec2 {  //点向量类
    Vec2() = default;
    Vec2(double _x, double _y) : x(_x), y(_y) {}
    /*Vec2 vector a - Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }
    /*Vec2 vector a + Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*cross product of two Vec2 vector (starting point is zero)*/
    friend inline double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;
    }
    /*Vec2 * num operation @return Vec2(x, y)*/
    friend inline Vec2 operator*(const Vec2& a, double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*input coordinate of Vec2 through istream*/
    inline friend istream& operator>>(istream& is, Vec2& p) {
        is >> p.x >> p.y;
        return is;
    }
    /*@param x x坐标 @param y y坐标*/
    double x = 0, y = 0;
} vertex[MAXN], its[MAXN]/*intersected points*/;

struct Seg {  //线段类 起点s, 终点t
    Seg() = default;
    Seg(Vec2 _s, Vec2 _t) : s(_s), t(_t) { theta = atan2((t - s).y, (t - s).x); }
    /*comparison of two Seg @return bool*/
    friend inline bool operator<(const Seg& a, const Seg& b) {//比较极角
        // theta1 != theta2 按theta排序
        if (abs(a.theta - b.theta) > EPS) return a.theta < b.theta;
        // theta1 == theta2 判断向量a在b的哪边，令最靠左的排在最左边
        return (b.s - a.s) * (b.t - a.s) > EPS;
    }
    /*intersection of two Seg @return Vec2(x, y)*/
    inline friend Vec2 intersect(const Seg& a, const Seg& b) {
        double ratio = ((b.t - b.s) * (a.s - b.s)) / ((a.t - a.s) * (b.t - b.s));
        return a.s + (a.t - a.s) * ratio;
    }
    /*@param s 起点 @param t 终点*/
    Vec2 s, t;
    /*@param theta 极角*/
    double theta = 1e9;
} seg[MAXN], q[MAXN]/*deque 队列中存放线段*/;
int n, cntVer, tot;
/*求向量左侧的半平面交 队列范围(l, r]*/
inline void halfPlaneIntersection() {
    sort(seg + 1, seg + tot + 1);
    int l = 0, r = 0;
    for (int i = 1; i <= tot; i++) {
        if (abs(seg[i].theta - seg[i - 1].theta) > EPS) {
            while (r - l > 1 && (seg[i].t - its[r]) * (seg[i].s - its[r]) > EPS)
                --r;
            while (r - l > 1 && (seg[i].t - its[l + 2]) * (seg[i].s - its[l + 2]) > EPS)
                ++l;
            q[++r] = seg[i];
            if (r - l > 1) its[r] = intersect(q[r], q[r - 1]); //求新交点
        }
    }
    while (r - l > 1 && (q[l + 1].t - its[r]) * (q[l + 1].s - its[r]) > EPS)
        --r; //删除多余元素
    its[r + 1] = intersect(q[l + 1], q[r]), ++r; //终点线段与起点线段相交

    double ans(0);
    for (int i = l + 2; i < r; i++) { // l + 2 是起始交点 三角剖分求面积
        ans += (its[i] * its[i + 1]) / 2;
    }
    ans += (its[r] * its[l + 2]) / 2;
    cout << setprecision(3) << fixed << ans << endl;
}

int main() {
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> cntVer;
        for (int j = 1; j <= cntVer; j++) cin >> vertex[j];
        for (int j = 1; j < cntVer; j++)
            seg[++tot] = Seg(vertex[j], vertex[j + 1]);
        seg[++tot] = Seg(vertex[cntVer], vertex[1]);
    }
    halfPlaneIntersection();
    return 0;
}

Description

UVA1304 Art Gallery

Solution

求凸多边形的核的面积，即求边的半平面交的面积。有的时候需要注意给出点的顺逆时针方向。

C ZJOI2008 瞭望塔

Description

题目链接
在山的轮廓线上选取一个位置建瞭望塔，使得在瞭望塔最顶端能看到山的每一个角落，求瞭望塔高的最小值。

Solution

先求出轮廓左侧半平面交(非闭合)，可以证明瞭望塔高取最小值时一定在 xix_i 或半平面交的交点所对应的 xx 处。
对于每个 xix_i，做垂直于水平方向的直线与半平面交的边界相交，求出 deltaYdeltaY。
对于半平面交的交点，同样做垂线，与轮廓线相交，求出 deltaYdeltaY。
注意：半平面交没有闭合，线段的取值范围为(l, r)。

Code

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

using namespace std;

constexpr double EPS = 1e-8;
constexpr int MAXN = 1e5 + 5;

struct Vec2 {  //点向量类
    Vec2() = default;
    Vec2(double _x, double _y) : x(_x), y(_y) {}
    /*Vec2 vector a - Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }
    /*Vec2 vector a + Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*cross product of two Vec2 vector (starting point is zero)*/
    friend inline double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;
    }
    /*Vec2 * num operation @return Vec2(x, y)*/
    friend inline Vec2 operator*(const Vec2& a, double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*input coordinate of Vec2 through istream*/
    inline friend istream& operator>>(istream& is, Vec2& p) {
        is >> p.x >> p.y;
        return is;
    }
    /*@param x x坐标 @param y y坐标*/
    double x = 0, y = 0;
} vertex[MAXN], its[MAXN];

struct Seg {  //线段类 起点s, 终点t
    Seg() = default;
    Seg(Vec2 _s, Vec2 _t) : s(_s), t(_t) {
        theta = atan2((t - s).y, (t - s).x);
    }
    /*comparison of two Seg @return bool*/
    friend inline bool operator<(const Seg& a, const Seg& b) {  //比较极角
        // theta1 != theta2 按theta排序
        if (abs(a.theta - b.theta) > EPS) return a.theta < b.theta;
        // theta1 == theta2 判断向量a在b的哪边，令最靠左的排在最左边
        return (b.s - a.s) * (b.t - a.s) > EPS;
    }
    /*intersection of two Seg @return Vec2(x, y)*/
    inline friend Vec2 intersect(const Seg& a, const Seg& b) {
        double ratio = ((b.t - b.s) * (a.s - b.s)) / ((a.t - a.s) * (b.t - b.s));
        return a.s + (a.t - a.s) * ratio;
    }
    /*@param s 起点 @param t 终点*/
    Vec2 s, t;
    /*@param theta 极角*/
    double theta = 1e9;
} seg[MAXN], q[MAXN];
int n, cntVer, tot, x[MAXN], y[MAXN];

inline void getHeight(int l, int r, double res = 1e15) {
    for (int i = 1; i <= n; i++) {
        double h = 0;
        for (int j = l + 1; j < r; j++) {
            auto tmpIts = intersect(q[j], Seg(Vec2(x[i], -1), Vec2(x[i], y[i])));
            h = max(h, tmpIts.y);
        }
        res = min(res, h - y[i]);
    }
    for (int i = 1; i < n; i++)
        for (int j = l + 2; j < r; j++) {
            if (its[j].x >= x[i] && its[j].x <= x[i + 1]) {
                auto tmpIts =
                    intersect(Seg(Vec2(its[j].x, -1), Vec2(its[j].x, its[j].y)),
                              Seg(Vec2(x[i], y[i]), Vec2(x[i + 1], y[i + 1])));
                res = min(res, its[j].y - tmpIts.y);
            }
        }
    cout << setprecision(3) << fixed << res << endl;
}
/*求向量左侧的半平面交 队列范围(l, r]*/
inline void halfPlaneIntersection() {
    sort(seg + 1, seg + tot + 1);
    int l = 0, r = 0;
    for (int i = 1; i <= tot; i++) {
        if (abs(seg[i].theta - seg[i - 1].theta) > EPS) {
            while (r - l > 1 && (seg[i].t - its[r]) * (seg[i].s - its[r]) > EPS)
                --r;
            while (r - l > 1 &&
                   (seg[i].t - its[l + 2]) * (seg[i].s - its[l + 2]) > EPS)
                ++l;
            q[++r] = seg[i];
            if (r - l > 1) its[r] = intersect(q[r], q[r - 1]);  //求新交点
        }
    }
    while (r - l > 1 && (q[l + 1].t - its[r]) * (q[l + 1].s - its[r]) > EPS)
        --r;                                      //删除多余元素
    its[r + 1] = intersect(q[l + 1], q[r]), ++r;  //终点线段与起点线段相交

    getHeight(l, r);
}

int main() {
    cin >> n;
    for (int i = 1; i <= n; i++) cin >> x[i];
    for (int i = 1; i <= n; i++) cin >> y[i];
    for (int i = 1; i < n; i++)
        seg[++tot] = Seg(Vec2(x[i], y[i]), Vec2(x[i + 1], y[i + 1]));
    halfPlaneIntersection();
    return 0;
}

D UVA1396 Most Distant Point from the Sea

Description

题目链接
求凸多边形的最大内切圆半径。

Solution

将多边形的边向内平移，当半平面交面积为 0 或不存在时，即所平移的距离为最大内切圆的半径。接下来二分移动的距离即可。

Code

Most Distant Point from the Sea
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
#include <bits/stdc++.h>
#define PI acos(-1)

using namespace std;
//求多边形最大内切圆半径 二分 + 半平面交
constexpr long double EPS = 1e-11;
constexpr int MAXN = 300;

struct Vec2 {  //点向量类
    Vec2() = default;
    Vec2(long double _x, long double _y) : x(_x), y(_y) {}
    /*Vec2 vector a - Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator-(const Vec2& a, const Vec2& b) {
        return Vec2(a.x - b.x, a.y - b.y);
    }
    /*Vec2 vector a + Vec2 vector b @return Vec2(x, y)*/
    friend inline Vec2 operator+(const Vec2& a, const Vec2& b) {
        return Vec2(a.x + b.x, a.y + b.y);
    }
    /*cross product of two Vec2 vector (starting point is zero)*/
    friend inline long double operator*(const Vec2& a, const Vec2& b) {
        return a.x * b.y - a.y * b.x;
    }
    /*Vec2 * num operation @return Vec2(x, y)*/
    friend inline Vec2 operator*(const Vec2& a, long double num) {
        return Vec2(a.x * num, a.y * num);
    }
    /*input coordinate of Vec2 through istream*/
    inline friend istream& operator>>(istream& is, Vec2& p) {
        is >> p.x >> p.y;
        return is;
    }
    /*@param x x坐标 @param y y坐标*/
    long double x = 0, y = 0;
} vertex[MAXN], its[MAXN]/*intersected points*/;

struct Seg {  //线段类 起点s, 终点t
    Seg() = default;
    Seg(Vec2 _s, Vec2 _t) : s(_s), t(_t) { theta = atan2((t - s).y, (t - s).x); }
    /*comparison of two Seg @return bool*/
    friend inline bool operator<(const Seg& a, const Seg& b) {//比较极角
        // theta1 != theta2 按theta排序
        if (abs(a.theta - b.theta) > EPS) return a.theta < b.theta;
        // theta1 == theta2 判断向量a在b的哪边，令最靠左的排在最左边
        return (b.s - a.s) * (b.t - a.s) > EPS;
    }
    /*intersection of two Seg @return Vec2(x, y)*/
    inline friend Vec2 intersect(const Seg& a, const Seg& b) {
        long double ratio = ((b.t - b.s) * (a.s - b.s)) / ((a.t - a.s) * (b.t - b.s));
        return a.s + (a.t - a.s) * ratio;
    }
    /*@param s 起点 @param t 终点*/
    Vec2 s, t;
    /*@param theta 极角*/
    long double theta = 1e9;
} seg[MAXN], sg[MAXN], q[MAXN]/*deque 队列中存放线段*/;
int n, tot;
/*求向量左侧的半平面交 队列范围[l, r]*/
inline bool halfPlaneIntersection() {
    sort(sg + 1, sg + tot + 1);
    int l = 0, r = 0;
    for (int i = 1; i <= tot; i++) {
        if (abs(sg[i].theta - sg[i - 1].theta) > EPS) {
            while (r - l > 1 && (sg[i].t - its[r]) * (sg[i].s - its[r]) > EPS)
                --r;
            while (r - l > 1 && (sg[i].t - its[l + 2]) * (sg[i].s - its[l + 2]) > EPS)
                ++l;
            q[++r] = sg[i];
            if (r - l > 1) its[r] = intersect(q[r], q[r - 1]); //求新交点
        }
    }
    while (r - l > 1 && (q[l + 1].t - its[r]) * (q[l + 1].s - its[r]) > EPS)
        --r; //删除多余元素
    its[r + 1] = intersect(q[l + 1], q[r]), ++r; //终点线段与起点线段相交

    long double ans(0);
    for (int i = l + 2; i < r; i++) { // l + 2 是起始交点 三角剖分求面积
        ans += (its[i] * its[i + 1]) / 2;
    }
    ans += (its[r] * its[l + 2]) / 2;
    return ans > EPS;
}

inline void translate(long double mid) { //平移 mid
    for (int i = 1; i <= tot; i++) {
        long double theta = seg[i].theta + PI / 2;
        long double deltaX = mid * cos(theta);
        long double deltaY = mid * sin(theta);
        auto delta = Vec2(deltaX, deltaY);
        sg[i] = Seg(seg[i].s + delta, seg[i].t + delta);
    }
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    while (cin >> n && n) {
        tot = 0;
        for (int i = 1; i <= n; i++) cin >> vertex[i];
        for (int i = 1; i <= n; i++)
            seg[++tot] = Seg(vertex[i], vertex[i % n + 1]);
        long double l = 0, r = 1e5;
        while (l < r - EPS) {
            long double mid = (l + r) / 2;
            translate(mid);
            if (halfPlaneIntersection()) l = mid;
            else r = mid;
        }
        cout << setprecision(8) << fixed << l << endl;
    }
    return 0;
}