[CQOI2006]凸多边形 /【模板】半平面交

洛谷

题意：逆时针给出\(n(n<=10)\)个凸多边形的顶点坐标，求它们交的面积。

学长博客，计算几何知识全面

半平面交问题详细讲解

还有一些前置知识。两向量\((x_1,y_1),(x_2,y_2)\)的叉乘为\(x_1y_2-x_2y_1\)，结果为正说明向量\((x_2,y_2)\)在向量\((x_1,y_1)\)逆时针方向，结果为负则在顺时针方向。求两直线交点的公式如下图所示：

代码与原博客稍有不同。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline int read() {
	char ch = getchar(); int x = 0, f = 1;
	while (ch < '0' || ch>'9') { if (ch == '-') f = -1; ch = getchar(); }
	while ('0' <= ch && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
	return x * f;
}
const int N = 1e3 + 10;
const double EPS = 1e-5;
int T, tot;
struct node {//一个点，两个坐标
	double x, y;
};
node p[15][55];
node ppx[N];
node operator - (node a, node b) {//两个点相减，得到向量ab
	node t;
	t.x = a.x - b.x;
	t.y = a.y - b.y;
	return t;
}
double operator ^ (node a, node b) { return a.x * b.y - a.y * b.x; }
//这里的ab应该是向量，然后求叉乘
struct Line {//一个向量，起点s，终点e
	node s, e;
};
Line L[N], que[N];
double getAngle(node a) { return atan2(a.y, a.x); }//这里的a应该是向量
double getAngle(Line a) { return atan2(a.e.y - a.s.y, a.e.x - a.s.x); }//求向量a的arctan值
bool cmp(Line a, Line b) {//对所有向量进行极角排序
	double A = getAngle(a);
	double B = getAngle(b);
	return A < B;
}
node getIntersectPoint(Line a, Line b) {//求两直线交点
	double a1 = a.s.y - a.e.y, b1 = a.s.x - a.e.x, c1 = 1.0 * a.s.x * a.e.y - 1.0 * a.e.x * a.s.y;
	double a2 = b.s.y - b.e.y, b2 = b.s.x - b.e.x, c2 = 1.0 * b.s.x * b.e.y - 1.0 * b.e.x * b.s.y;
	node t;
	t.x = (1.0 * c1 * b2 - 1.0 * c2 * b1) / (1.0 * a2 * b1 - 1.0 * a1 * b2);
	t.y = (1.0 * c1 * a2 - 1.0 * c2 * a1) / (1.0 * a2 * b1 - 1.0 * a1 * b2);
	return t;
}
//判断向量a是否在向量bc交点的右侧
bool onRight(Line a, Line b, Line c) {
	node o = getIntersectPoint(b, c);
	if (((a.e - a.s) ^ (o - a.s)) < 0) return true;//可以自己画图a.s a.e o三个点
	return false;
}
double HalfPlaneIntersection() {
	sort(L + 1, L + tot + 1, cmp);
	int head = 1, tail = 1;
	que[1] = L[1];//构造双端队列
	for (int i = 2; i <= tot; i++) {
		while (head < tail && onRight(L[i], que[tail], que[tail - 1])) tail--;
		while (head < tail && onRight(L[i], que[head], que[head + 1])) head++;
		que[++tail] = L[i];
//极角相同的向量，保留靠左的那一个
		if (fabs(getAngle(que[tail]) - getAngle(que[tail - 1])) < EPS) {
			tail--;
			if (((que[tail].e - que[tail].s) ^ (L[i].e - que[tail].s)) > EPS)que[tail] = L[i];
		}
	}
	while (head < tail && onRight(que[head], que[tail], que[tail - 1])) tail--;
	while (head < tail && onRight(que[tail], que[head], que[head + 1])) head++;
	if (tail - head < 2) return 0;//剩下的直线无法构成多边形
	double ans = 0;
	int tot_jd = 0;
	for (int i = head; i < tail; ++i) {
		ppx[++tot_jd] = getIntersectPoint(que[i], que[i + 1]);
	}
	ppx[++tot_jd] = getIntersectPoint(que[tail], que[head]);
	for (int i = 2; i < tot_jd; ++i) {
		double x1 = ppx[1].x, y1 = ppx[1].y;
		double x2 = ppx[i].x, y2 = ppx[i].y;
		double x3 = ppx[i + 1].x, y3 = ppx[i + 1].y;
		ans = ans + (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1);
	}
	return ans / 2.0;
}
int main() {
	int T; cin >> T;
	for (int t = 1; t <= T; ++t) {
		int n; cin >> n;
		for (int i = 1; i <= n; ++i) {
			cin >> p[t][i].x;
			cin >> p[t][i].y;
		}
		for (int i = 1; i < n; i++) {
			L[++tot].s.x = p[t][i].x;
			L[tot].s.y = p[t][i].y;
			L[tot].e.x = p[t][i + 1].x;
			L[tot].e.y = p[t][i + 1].y;
		}
		L[++tot].s.x = p[t][n].x;
		L[tot].s.y = p[t][n].y;
		L[tot].e.x = p[t][1].x;
		L[tot].e.y = p[t][1].y;
	}
	printf("%.3lf\n", HalfPlaneIntersection());
	return 0;
}