题面
逆时针顺序给出一些凸多边形,求它们交的面积
多边形个数2 <= n <= 10,每一个多边形边数3 <= mi <= 50
分析
半平面交其实就类似数学中的线性规划,从给出的一些域中找到所有域都满足的可行域。
数据量基本能做 O(n2) ,但是选择 O(n) 的做法
用凸多边形的逆时针顶点序来记录可行域的顶点。
维护时依据的是按极角序逐步加入有向直线,并且删去已经无效的直线(相当于删去可行域顶点)
图中取的都是有向线段左侧为可行域:
按极角序即1,2,3,4顺序加入,当加入4时候,发现半平面的上一个顶点 A 在4右侧,说明这是个无效顶点,删去,并且加入新端点 B
还有一种可能,是加入新线后,转了一圈,导致第一条直线无效而删除。
像这个,新加入4后,发现只有B下侧区域才有效,所以 A 也被删除了,但 A 是队首元素。
平行时,取靠内的(更强的约束条件)
所以发现要维护可行域多边形,需要用双端队列。
对于判断点在有向直线右侧, 需要用到外积知识:
A 和 u 是为了记录这个有向直线信息的,A 是直线上一点,u 是直线上的某个向量,从而 A+tu 可以表示直线上所有点。
在这里发现
u
⃗
×
A
P
⃗
>
0
\vec u × \vec{AP}>0
u×AP>0当且仅当 P 在线左侧,且A点只要在直线上可以随意取
这就完成了点在线左侧的判定
代码
总体上还是比较复杂
#include <iostream>
#include<cstdio>
#include<algorithm>
#include<string.h>
#include<string>
#include<math.h>
#include<vector>
#include<iomanip>
using namespace std;
class Vector {//向量
public:
double x, y;
Vector() {}
Vector(double x, double y) :x(x), y(y) {}
double Cross(Vector other)
{
return x * other.y - y * other.x;
}
};
class Point {//点
public:
double x, y;
Point(double x, double y) :x(x), y(y) {}
Point() {}
Vector operator -(const Point & other) const
{
return Vector(x - other.x, y - other.y);
}
Point operator + (const Point& other)const
{
return Point(x + other.x, y + other.y);
}
};
class Line {//有向直线,左侧为选中范围
public:
Point P;
Vector v;
double ang;
Line(){}
Line(Point p, Vector v) :P(p), v(v) { ang = atan2(v.y, v.x); }
bool operator < (const Line& L)const {
return ang < L.ang;//按角度排序,从小到大
}
};
class HalfplaneIntersection {//半平面交
private:
Line q[600];//双端队列保存交集用到的线
Point p[600];//保存q[i]和q[i+1]的交点
int first, last;//双端队列指针
public:
vector<Point> poly;//存有效面积的端点
bool OnLeft(Line L, Point p)//判定点p在线L的左边
{
return L.v.Cross(p-L.P) > 0;
}
Point GetInerserction(Line a, Line b)//求两Line的交点
{
Vector u = a.P - b.P;
double t = b.v.Cross(u) / a.v.Cross(b.v);
return Point(a.P.x + a.v.x * t, a.P.y + a.v.y * t);
}
void run(Line* L, int n)
{
sort(L, L + n);//按极角排序
q[first = last = 0] = L[0];
for (int i = 1; i < n; i++)
{
while (first < last && !OnLeft(L[i], p[last - 1]))last--;//删去末尾无效点(p[last]是上一个首尾交点,所以在这不考虑)
while (first < last && !OnLeft(L[i], p[first]))first++;//删去开头无效点(转回来一圈)
q[++last] = L[i];//加入这条线
if (fabs(q[last].v.Cross(q[last - 1].v))<1e-8) {//和上次加入的平行了,取更靠左侧的
last--;
if (OnLeft(q[last], L[i].P))q[last] = L[i];
}
if (first < last)p[last - 1] = GetInerserction(q[last - 1], q[last]);
}
while (first < last && !OnLeft(q[first], p[last - 1]))last--;//删去末尾无效点
if (last - first <= 1)return;//空集
p[last] = GetInerserction(q[last], q[first]);//计算首尾交点,加入
for (int i = first; i <= last; i++)
{
poly.push_back(p[i]);//保存点集
}
}
}HPI;
int pos[2][55];//点数
Line L[600];
int cnt = 0;
int main()
{
ios::sync_with_stdio(false);
int n, m;
cin >> n;
for (int o = 0; o < n; o++)//每次输入一个凸多边形,逆时针顺序
{
cin >> m;
for (int i = 0; i < m; i++)cin >> pos[0][i] >> pos[1][i];
for (int i = 1; i < m; i++)L[cnt++] = Line(Point(pos[0][i],pos[1][i]), Vector(pos[0][i] - pos[0][i - 1], pos[1][i] - pos[1][i - 1]));
L[cnt++]= Line(Point(pos[0][0], pos[1][0]), Vector(pos[0][0] - pos[0][m-1], pos[1][0] - pos[1][m - 1]));
}
HPI.run(L, cnt);
//HPI下的poly现在保存了可行域的顶点
double ans = 0.0;
Point pre = Point(0, 0);//从pre开始,分割三角形用外积求面积
int siz = HPI.poly.size();
for (int i = 1; i <= HPI.poly.size(); i++)
{
ans += (HPI.poly[i - 1] - pre).Cross(HPI.poly[i%siz] - pre);
}
ans /= 2;//外积结果是平行四边形,三角形二倍
if (ans < 0)ans = -ans;
cout << fixed<<setprecision(3)<<ans;
return 0;
}