半平面交

前言

由于笔者对半平面交算法的很多数学细节不是非常清楚，所以这篇博客可能会有描述不到位或者错误的地方。欢迎读者向笔者指出文章中的不足，也希望初学半平面交的同学不要用这篇博客来入门，以免产生了一些错误的理解，对计算几何之后的学习产生了影响。

前置知识

我们知道一条直线可以将一个平面分成两部分，我们称其中的一部分为 半平面 。半平面是一个 点集 ，也就是由直线和直线的一侧构成的点集。当半平面包含这条直线时，我们称其为 闭半平面，反之称其为 开半平面。半平面在计算几何中用 向量 来表示，且统一以向量的某一侧为半平面。

多个半平面的交集称之为 半平面交。半平面交仍然是一个点集，它在平面直角坐标系中围成了一个区域。如果把多边形的边看作是首尾相连的向量，那么这些向量在多边形内侧方向的半平面交称之为 多边形的核。

接下来引入和 极角排序 有关的概念。在平面上任取一点 \(O\)，取其为 极点。再引出一条射线 \(Ox\) ，称为 极轴。设平面上有任意一点 \(M\)，令线段 \(OM\) 的长度为 \(p\) ，\(\theta\) 表示 \(Ox\) 到 \(OM\) 的夹角。我们称 \(p\) 是点 \(M\) 的极径，\(\theta\) 是点 \(M\) 的极角。我们在极角排序中取极轴为 \(x\) 轴，这样我们就可以把平面上的若干个点按照与 \(x\) 轴的夹角排序。具体实现可以使用 C 语言中 math.h 库里的 atan2 函数，调用格式为 atan2(double y, double x) ，其中 y 表示已知点的 \(y\) 坐标，x 同理。这个函数的返回值是这个点与远点连边后，这条连边与 \(x\) 轴正方向的夹角。

\(S\&I\) 算法

假设给定了多边形的顶点坐标，需要求出这个多边形的核。假设给定了下图中的多边形。

我们可以先处理出这个多边形的所有边，我们定义 逆时针方向 为正方向，然后把这些边极角排序。这样做可以方便后面双端队列的处理。

我们可以先假定这个平面上所有的面积都属于半平面交，然后考虑多边形的边对于半平面交的影响。我们可以把多边形的边转化成若干向量，然后对于这些向量的左侧取交集。我们用一个双端队列来维护对当前的半平面交有影响的边，同时维护队列中相邻边的交点。

设双端队列中的第 \(k\) 条边为 \(q_k\) ，且当前共有 \(m\) 条边。对于当前需要加入的边 \(p_i\) 和上一次的交点 \(t\) ，如果 \(p\) 在 \(t\) 的逆时针方向，说明当前的边 \(p\) 对于半平面交有影响，而构成上一次的交点的连边，也就是 \(q_{m - 1}\) 对半平面交没有影响，所以需要弹出。

当半平面交处理到最后几条边的时候，这些边可能会覆盖掉双端队列开头的影响。这时我们还需要判断是否需要弹出开头的边，例如下图中的情况：

当编号为 \(8\) 的边加入双端队列以后，编号为 \(1\) 的边对半平面交就没有影响了。因为边已经极角排序过了，所以此时编号为 \(1\) 的边应该在双端队列开头。因此我们需要比较队尾和开头，决定是否需要将开头弹出。另外，我们应该先处理队尾，再处理开头。我们可以感性理解一下：因为先处理开头可能会导致合法的范围扩大，因而原本对半平面交没有影响的边也会加入双端队列。因为队头被其后面的元素约束，队尾则没有，所以队尾应该优先考虑。

如果出现两条平行线，我们应该取在另一条平行线逆时针方向的平行线。因为我们取的是向量的左侧，所以这样做显然不会破坏算法的正确性。

最后，我们应该还要用队首来过滤掉队尾的不合法元素。假如队尾的元素在队首的右侧，那么队尾的元素就是没有影响的。因为我们对边进行过极角排序，所以可能会出现环。

当多边形的核的边数 \(\geq 3\) 时，我们就可以求出这个多边形的核的面积了。具体实现我们可以利用叉积的几何意义，因为笔者也不太熟悉叉积，因此无法对原理进行解释。具体请见代码注释。

参考代码

例题链接

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const int maxn = 1e3 + 5;
const double eps = 1e-10;

int n, m, cnt;
int st, en;
double ans;

struct vec // 定义向量和基本运算 
{
	double x, y;
	vec(): x(), y() {}
	vec(double _x, double _y): x(_x), y(_y) {}
	vec operator + (const vec& rhs) const
	{
		return vec(x + rhs.x, y + rhs.y);
	}
	vec operator - (const vec& rhs) const
	{
		return vec(x - rhs.x, y - rhs.y);
	}
	vec operator += (const vec& rhs)
	{
		this -> x += rhs.x;
		this -> y += rhs.y;
	}
	vec operator -= (const vec& rhs)
	{
		this -> x -= rhs.x;
		this -> y -= rhs.y;
	}
	vec operator * (const double& rhs) const
	{
		return vec(x * rhs, y * rhs);
	}
	vec operator / (const double& rhs) const
	{
		return vec(x / rhs, y / rhs);
	}
	vec operator *= (const double& rhs)
	{
		this -> x *= rhs;
		this -> y *= rhs;
	}
	vec operator /= (const double& rhs)
	{
		this -> x /= rhs;
		this -> y /= rhs; 
	}
} a[maxn], t[maxn]; // a ->读入的若干多边形的顶点，t ->半平面交中相邻边的交点 

struct line // 定义多边形的边
{
	vec a, b;
	double k;
	line(): a(), b() {}
	line(vec _a, vec _b): a(_a), b(_b) {k = atan2(b.y, b.x);} // k ->极角 
	bool operator < (const line& rhs) const // 极角排序 
	{
		return k < rhs.k;
	}
} p[maxn], q[maxn]; // p ->若干多边形的边，q ->双端队列 

double cross(vec a, vec b) // 计算叉积 
{
	return a.x * b.y - b.x * a.y;
}

vec get(line a, line b) // 计算边 a 和边 b 的交点 
{
	vec c = a.a - b.a;
	double t = cross(b.b, c) / cross(a.b, b.b);
	return a.a + a.b * t; 
}

void work() // 求半平面交 
{
	st = en = 1;
	q[st] = p[1];
	for (int i = 2; i <= cnt; i++)
	{
		while (st < en && cross(p[i].b, t[en - 1] - p[i].a) <= eps) // 当前边是否在队尾交点的逆时针方向 
			en--;
		while (st < en && cross(p[i].b, t[st] - p[i].a) <= eps) // 当前边是否在队首交点的逆时针方向 
			st++;
		q[++en] = p[i];
		if (fabs(cross(q[en].b, q[en - 1].b)) <= eps) // 判断平行线的情况 
		{
			en--;
			if (cross(q[en].b, p[i].a - q[en].a) > eps) // 当前边在平行线的逆时针方向（合法） 
				q[en] = p[i];
		}
		if (st < en) // 如果边数 > 1 则加入交点 
			t[en - 1] = get(q[en - 1], q[en]);
	}
	while (st < en && cross(q[st].b, t[en - 1] - q[st].a) <= eps) // 用队首来过滤队尾 
		--en;
	if (en - st <= 1) // 半平面交是否存在面积 
		return;
	t[en] = get(q[st], q[en]);
}

int main()
{
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
	{
		scanf("%d", &m);
		for (int j = 1; j <= m; j++)
			scanf("%lf%lf", &a[j].x, &a[j].y);
		for (int j = 1; j <= m; j++)
		{
			cnt++;
			int to = (j == m ? 1 : j + 1);
			p[cnt] = line(a[j], a[to] - a[j]);
		}
	}
	sort(p + 1, p + cnt + 1);
	work();
	for (int i = st; i <= en; i++)
	{
		int to = (i == en ? st : i + 1);
		ans += cross(t[i], t[to]); // 叉积求半平面交面积 
	}
	ans /= 2;
	printf("%.3lf\n", ans);
	return 0;
}