欢迎关注更多精彩
关注我，学习常用算法与数据结构，一题多解，降维打击。
本期话题：半平面求交
背景知识
学习资料
视频讲解
https://www.bilibili.com/video/BV1jL411C7Ct/?spm_id_from=333.1007.top_right_bar_window_history.content.click&vd_source=fb27f95f25902a2cc94d4d8e49f5f777
文本资料
https://oi-wiki.org//geometry/half-plane/
基本问题转化
在很多题目中，给定的线段是没有方向的。此时，我们需要先把所有的线段都转化成点加向量的方式。使得向量的左边为有效区域。这样就可以使用模板求解了。
要注意的问题
- 主要的问题是浮点型的判断大小问题。在排序和判断点与线的关系时都用到浮点型判断。有些题型会卡精度，能用整数判断尽量不要使用浮点判断。
- atan2计算比较耗时，可以事先保存。
代码模板
求多边形的核
题目链接：https://vjudge.net/problem/UVA-1571
多边形的核就是取核区域内任意一点，站在该可以观察到多边形内任意一点。
利用半平面求交可以得到多边形的核
#include<stdio.h>
#include<cmath>
#include <algorithm>
#include <vector>
#include <list>
#include <cstring>
#include <set>
using namespace std;
const double EPS = 1e-14;
const int N = 2e6 + 10;
int cmp(double d) {
if (abs(d) < EPS)return 0;
if (d > 0)return 1;
return -1;
}
class Point {
public:
double x, y;
int id;
Point() {}
Point(double a, double b) :x(a), y(b) {}
Point(const Point& p) :x(p.x), y(p.y), id(p.id) {}
void in() {
scanf("%lf %lf", &x, &y);
}
void out() {
printf("%.16f %.16f\n", x, y);
}
double dis() {
return sqrt(x * x + y * y);
}
double dis2() {
return x * x + y * y;
}
Point operator -() const {
return Point(-x, -y);
}
Point operator -(const Point& p) const {
return Point(x - p.x, y - p.y);
}
Point operator +(const Point& p) const {
return Point(x + p.x, y + p.y);
}
Point operator *(double d)const {
return Point(x * d, y * d);
}
Point operator /(double d)const {
return Point(x / d, y / d);
}
void operator -=(Point& p) {
x -= p.x;
y -= p.y;
}
void operator +=(Point& p) {
x += p.x;
y += p.y;
}
void operator *=(double d) {
x *= d;
y *= d;
}
void operator /=(double d) {
this ->operator*= (1 / d);
}
bool operator<(const Point& a) const {
return x < a.x || (abs(x - a.x) < EPS && y < a.y);
}
bool operator==(const Point& a) const {
return abs(x - a.x) < EPS && abs(y - a.y) < EPS;
}
};
// 向量操作
double cross(const Point& a, const Point& b) {
return a.x * b.y - a.y * b.x;
}
double dot(const Point& a, const Point& b) {
return a.x * b.x + a.y * b.y;
}
class Line {
public:
Point front, tail;
double ang;
int u, v;
Line() {}
Line(const Point& a, const Point& b) :front(a), tail(b) {
ang = atan2(front.y - tail.y, front.x - tail.x);
}
};
int cmp(const Line& a, const Line& b) {
//if (a.u == b.u && a.v == b.v)return 0;
return cmp(a.ang - b.ang);
}
// 点在直线哪一边>0 左边，<0边
double SideJudge(const Line& a, const Point& b) {
//return cmp(cross(a.front - a.tail, b - a.tail));
return cross(a.front - a.tail, b - a.tail);
}
int LineSort(const Line& a, const Line& b) {
int c = cmp(a, b);
if (c)return c < 0;
return cross(b.front - b.tail, a.front - b.tail) > 0;
}
/*
点p 到 p+r 表示线段1
点q 到 q+s 表示线段2
线段1 上1点用 p' = p+t*r (0<=t<=1)
线段2 上1点用 q' = q+u*s (0<=u<=1)
让两式相等求交点 p+t*r = q+u*s
两边都叉乘s
(p+t*r)Xs = (q+u*s)Xs
pXs + t*rXs = qXs
t = (q-p)Xs/(rXs)
同理，
u = (p-q)Xr/(sXr) -> u = (q-p)Xr/(rXs)
以下分4种情况：
1. 共线，sXr==0 && (q-p)Xr==0, 计算 (q-p)在r上的投影在r长度上的占比t0，
计算(q+s-p)在r上的投影在r长度上的占比t1，查看[t0, t1]是否与范围[0，1]有交集。
如果t0>t1, 则比较[t1, t0]是否与范围[0，1]有交集。
t0 = (q-p)*r/(r*r)
t1 = (q+s-p)*r/(r*r) = t0 + s · r / (r · r)
2. 平行sXr==0 && (q-p)Xr!=0
3. 0<=u<=1 && 0<=t<=1 有交点
4. 其他u, t不在0到范围内，没有交点。
*/
pair<double, double> intersection(const Point& q, const Point& s, const Point& p, const Point& r, bool &oneline) {
// 计算 (q-p)Xr
auto qpr = cross(q - p, r);
auto qps = cross(q - p, s);
auto rXs = cross(r, s);
if (cmp(rXs) == 0) {
oneline = true;
return { -1, -1 }; // 平行或共线
}
// 求解t, u
// t = (q-p)Xs/(rXs)
auto t = qps / rXs;
// u = (q-p)Xr/(rXs)
auto u = qpr / rXs;
return { u, t };
}
Point LineCross(const Line& a, const Line& b, bool &f) {
Point dira = a.front - a.tail;
Point dirb = b.front - b.tail;
bool oneline=false;
auto p = intersection(a.tail, dira, b.tail, dirb, oneline);
if (oneline)f = false;
return a.tail + dira * p.first;
}
class HalfPlane {
public:
vector<Line> lines;
vector<int> q;
vector<Point> t;
int len;
HalfPlane() {
lines.resize(N);
q.resize(N);
t.resize(N);
}
void reset() {
len = 0;
}
void addLine(const Line& a) {
lines[len++] = a;
}
bool run() {
sort(lines.begin(), lines.begin() + len, LineSort);
int l = -1, r = 0;
q[0] = 0;
for (int i = 1; i < len; ++i) {
if (cmp(lines[i], lines[i - 1]) == 0)continue;
while (r - l > 1 && SideJudge(lines[i], t[r]) < 0)r--;
while (r - l > 1 && SideJudge(lines[i], t[l + 2]) < 0)l++;
q[++r] = i;
bool f=true;
t[r] = LineCross(lines[q[r]], lines[q[r - 1]], f);
}
while (r - l > 1 && SideJudge(lines[q[l + 1]], t[r]) < 0)r--;
//if (r - l > 1) {
// bool f = true;
// t[r + 1] = LineCross(lines[q[l + 1]], lines[q[r]], f);
// r++;
// if (!f)r -= 2;
//}
统计交点
//l++;
//vector<Point> ans(r - l);
//for (int i = 0; i < ans.size(); ++i) {
// ans[i] = t[i + l + 1];
//}
return r-l>2;
}
};
Point oiPs[N * 2];
pair<int, int> ori[N * 2];
HalfPlane hp;
int bigDevid(int a, int b) {
for (int i = max(abs(a), abs(b)); i >= 1; i--) {
if (a % i == 0 && b % i == 0)return i;
}
return 1;
}
void solve() {
int n, m = 1;
//FILE* fp = fopen("ans.txt", "w");
while (scanf("%d", &n) != EOF && n) {
int a, b;
for (int i = 0; i < n; ++i) {
scanf("%d %d", &a, &b);
oiPs[i] = Point(a, b);
ori[i] = { a,b };
}
oiPs[n] = oiPs[0];
ori[n] = ori[0];
hp.reset();
for (int i = 0; i < n; ++i) {
hp.addLine(Line(oiPs[i+1], oiPs[i]));
hp.lines[i].u = ori[i+1].first - ori[i].first;
hp.lines[i].v = ori[i+1].second - ori[i].second;
int bd = bigDevid(hp.lines[i].u, hp.lines[i].v);
hp.lines[i].u /= bd;
hp.lines[i].v /= bd;
}
auto ps = hp.run();
if (ps)puts("1");
else puts("0");
m++;
}
}
int main() {
solve();
return 0;
}
/*
4
0 0
0 1
1 1
1 0
8
0 0
3 0
4 3
2 2
3 4
4 4
4 5
0 5
0
8
0 0
0 1
1 1
1 2
0 2
0 3
3 3
3 0
*/
练习一
链接：https://www.luogu.com.cn/problem/P4196
求多个凸多边形的交面积。
对每条边进行半平面求交，再利用三角形求多边形面积。
#include<stdio.h>
#include<cmath>
#include <algorithm>
#include <vector>
#include <list>
#include <cstring>
#include <set>
using namespace std;
const double EPS = 1e-14;
const int N = 2e6 + 10;
int cmp(double d) {
if (abs(d) < EPS)return 0;
if (d > 0)return 1;
return -1;
}
class Point {
public:
double x, y;
int id;
Point() {}
Point(double a, double b) :x(a), y(b) {}
Point(const Point& p) :x(p.x), y(p.y), id(p.id) {}
void in() {
scanf("%lf %lf", &x, &y);
}
void out() {
printf("%.16f %.16f\n", x, y);
}
double dis() {
return sqrt(x * x + y * y);
}
double dis2() {
return x * x + y * y;
}
Point operator -() const {
return Point(-x, -y);
}
Point operator -(const Point& p) const {
return Point(x - p.x, y - p.y);
}
Point operator +(const Point& p) const {
return Point(x + p.x, y + p.y);
}
Point operator *(double d)const {
return Point(x * d, y * d);
}
Point operator /(double d)const {
return Point(x / d, y / d);
}
void operator -=(Point& p) {
x -= p.x;
y -= p.y;
}
void operator +=(Point& p) {
x += p.x;
y += p.y;
}
void operator *=(double d) {
x *= d;
y *= d;
}
void operator /=(double d) {
this ->operator*= (1 / d);
}
bool operator<(const Point& a) const {
return x < a.x || (abs(x - a.x) < EPS && y < a.y);
}
bool operator==(const Point& a) const {
return abs(x - a.x) < EPS && abs(y - a.y) < EPS;
}
};
// 向量操作
double cross(const Point& a, const Point& b) {
return a.x * b.y - a.y * b.x;
}
double dot(const Point& a, const Point& b) {
return a.x * b.x + a.y * b.y;
}
class Line {
public:
Point front, tail;
double ang;
int u, v;
Line() {}
Line(const Point& a, const Point& b) :front(a), tail(b) {
ang = atan2(front.y - tail.y, front.x - tail.x);
}
};
int cmp(const Line& a, const Line& b) {
//if (a.u == b.u && a.v == b.v)return 0;
return cmp(a.ang - b.ang);
}
// 点在直线哪一边>0 左边，<0边
double SideJudge(const Line& a, const Point& b) {
//return cmp(cross(a.front - a.tail, b - a.tail));
return cross(a.front - a.tail, b - a.tail);
}
int LineSort(const Line& a, const Line& b) {
int c = cmp(a, b);
if (c)return c < 0;
return cross(b.front - b.tail, a.front - b.tail) > 0;
}
/*
点p 到 p+r 表示线段1
点q 到 q+s 表示线段2
线段1 上1点用 p' = p+t*r (0<=t<=1)
线段2 上1点用 q' = q+u*s (0<=u<=1)
让两式相等求交点 p+t*r = q+u*s
两边都叉乘s
(p+t*r)Xs = (q+u*s)Xs
pXs + t*rXs = qXs
t = (q-p)Xs/(rXs)
同理，
u = (p-q)Xr/(sXr) -> u = (q-p)Xr/(rXs)
以下分4种情况：
1. 共线，sXr==0 && (q-p)Xr==0, 计算 (q-p)在r上的投影在r长度上的占比t0，
计算(q+s-p)在r上的投影在r长度上的占比t1，查看[t0, t1]是否与范围[0，1]有交集。
如果t0>t1, 则比较[t1, t0]是否与范围[0，1]有交集。
t0 = (q-p)*r/(r*r)
t1 = (q+s-p)*r/(r*r) = t0 + s · r / (r · r)
2. 平行sXr==0 && (q-p)Xr!=0
3. 0<=u<=1 && 0<=t<=1 有交点
4. 其他u, t不在0到范围内，没有交点。
*/
pair<double, double> intersection(const Point& q, const Point& s, const Point& p, const Point& r, bool &oneline) {
// 计算 (q-p)Xr
auto qpr = cross(q - p, r);
auto qps = cross(q - p, s);
auto rXs = cross(r, s);
if (cmp(rXs) == 0) {
oneline = true;
return { -1, -1 }; // 平行或共线
}
// 求解t, u
// t = (q-p)Xs/(rXs)
auto t = qps / rXs;
// u = (q-p)Xr/(rXs)
auto u = qpr / rXs;
return { u, t };
}
Point LineCross(const Line& a, const Line& b, bool &f) {
Point dira = a.front - a.tail;
Point dirb = b.front - b.tail;
bool oneline=false;
auto p = intersection(a.tail, dira, b.tail, dirb, oneline);
if (oneline)f = false;
return a.tail + dira * p.first;
}
class HalfPlane {
public:
vector<Line> lines;
void addLine(const Line& a) {
lines.push_back(a);
}
vector<Point> run() {
sort(lines.begin(), lines.end(), LineSort);
vector<int> q(lines.size() + 10);
vector<Point> t(lines.size() + 10);
int l = -1, r = 0;
q[0] = 0;
for (int i = 1; i < lines.size(); ++i) {
if (cmp(lines[i], lines[i - 1]) == 0)continue;
while (r - l > 1 && SideJudge(lines[i], t[r]) < 0)r--;
while (r - l > 1 && SideJudge(lines[i], t[l + 2]) < 0)l++;
q[++r] = i;
bool f = true;
t[r] = LineCross(lines[q[r]], lines[q[r - 1]], f);
}
while (r - l > 1 && SideJudge(lines[q[l + 1]], t[r]) < 0)r--;
if (r - l > 1) {
bool f = true;
t[r + 1] = LineCross(lines[q[l + 1]], lines[q[r]], f);
r++;
}
// 统计交点
l++;
vector<Point> ans(r - l);
for (int i = 0; i < ans.size(); ++i) {
ans[i] = t[i + l + 1];
}
return ans;
}
};
Point oiPs[N];
void solve() {
int n, m;
scanf("%d", &n);
HalfPlane hp;
int a, b;
while (n--) {
scanf("%d", &m);
for (int i = 0; i < m; ++i) {
scanf("%d%d", &a, &b);
oiPs[i].x = a;
oiPs[i].y = b;
}
oiPs[m] = oiPs[0];
for (int i = 0; i < m; ++i) {
hp.addLine(Line(oiPs[i + 1], oiPs[i]));
}
}
auto keyPoints = hp.run();
double ans = 0;
for (int i = 2; i < keyPoints.size(); ++i) {
ans += cross(keyPoints[i - 1] - keyPoints[0], keyPoints[i] - keyPoints[0]);
}
printf("%.3f\n", ans / 2);
}
int main() {
solve();
return 0;
}
/*
3
3
-1 2
-2 1
-1 1
3
1 1
2 1
1 2
3
1 1
3 0
2 2
*/
本人码农，希望通过自己的分享，让大家更容易学懂计算机知识。创作不易，帮忙点击公众号的链接。
