P4196 半平面交

题面

逆时针顺序给出一些凸多边形,求它们交的面积
多边形个数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;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值