关于凸包位置关系的判断

近日恰好和同学谈到多边形之间怎么判断相交关系，便写下这篇博文。

由于非凸多边形的不确定性，这里就只谈论凸多边形间位置关系判断的优化。对于分别有 $n$ 和 $m$ 条边的非凸多边形可以枚举两个多边形的边判断线段是否相交，时间复杂度为 $O (m n)$ 。

凸多边形（以下简称凸包）也可以通过枚举边判断相交关系来判断两凸包是否相交，但 $O (m n)$ 的复杂度有时属实过大。下面对如何优化进行讨论。

对于上图两个凸包，我们将凸包B绕着凸包A移动，并绘制出凸包B最左下角点的轨迹：

不难发现，最后凸包B最左下角点的轨迹为两凸包边重新拼接后形成的大凸包，即闵可夫斯基和。对于为什么绘制两凸包的闵可夫斯基和，可以从图中看出，当凸包B最左下角点位于该大凸包内时，两个凸包相交。这里判断点在凸包内则可以用叉乘解决，最终复杂度为 $O ((m + n) \log (m + n) + (m + n))$ （包括求闵可夫斯基和对边进行极角排序和遍历进行叉积判断单点是否在凸包内），即 $O ((m + n) \log (m + n))$ 。

#include<bits/stdc++.h>
#define pii pair<int, int>
#define double long double
#define eps 1e-9

using namespace std;

//使用结构体表示向量方便进行极角排序
struct vec
{
	pii v;
	double ang;

	vec(pii v)
	{
		this->v = v;
		ang = atan2(v.second, v.first);
	}

	bool operator < (vec& r)
	{
		return this->ang < r.ang;
	}
};

pii operator + (pii& l, pii& r)
{
	return make_pair(l.first + r.first, l.second + r.second);
}

pii operator - (pii& l, pii& r)
{
	return make_pair(l.first - r.first, l.second - r.second);
}

bool cmp(vec l, vec r)
{
	return l.ang < r.ang;
}

//向量叉乘
double cross(pii l, pii r)
{
	return l.first * r.second - l.second * r.first;
}

int posOfCh(vector<pii> pochA, int size1, vector<pii> pochB, int size2) //pochA, pochB分别按顺时针顺序存储两凸包上的点
{
	vector<vec> vecs;
	pii maxp = pochA[0], minp = pochB[0];
	for (int i = 0; i < size1; i++)
	{
		vecs.push_back(vec(pochA[(i + 1) % size1] - pochA[i])); //将凸包的边按向量的形式存入，凸包A按顺时针，凸包B按逆时针
		maxp = max(maxp, pochA[i]); //同时查找凸包A中的最右上角点
	}
	for (int i = size2; i > 0; i--)
	{
		vecs.push_back(vec(pochB[i - 1] - pochB[i % size2]));
		minp = min(minp, pochB[i - 1]); //同时查找凸包B中的最左下角点
	}
	sort(vecs.begin(), vecs.end()); //极角排序
	
	//找到从极角-π/2开始（包括-π/2）顺时针方向第一条向量
	int pos = upper_bound(vecs.begin(), vecs.end(), vec({ 0, -1 }), cmp) - vecs.begin() - 1;
	
	//计算闵可夫斯基和的凸包顶点
	vector<pii> pomsch;
	pomsch.push_back(maxp);
	for (int i = pos, j = 0; j < vecs.size(); i = (i + vecs.size() - 1) % vecs.size(), j++)
	{
		pii p = pomsch.back() + vecs[i].v;
		while (pomsch.size() > 1 && cross(vecs[i].v, pomsch.back() - pomsch[pomsch.size() - 2]) < eps) //去除三点共线的中间点
			pomsch.pop_back();

		pomsch.push_back(p);
	}
	pomsch.pop_back();

	//叉积判断凸包B最左下角点是否在闵可夫斯基和中
	int flag = 1;
	for (int i = 0; i < pomsch.size(); i++)
		if (cross(minp - pomsch[i], pomsch[(i + 1) % pomsch.size()] - pomsch[i]) < eps)
		{
			flag = 0;
			break;
		}

	return flag;
}

上面讲述了如何判断凸包间的相交关系，但如果两凸包之间是包含关系（一个凸包的边与另一个凸包的边无任何交点但两凸包并非相离关系），以上代码依然会返回整数1。因此，我们再次以相似的方式画出凸包B最左下角点的轨迹来分析。

在这张图中我们较难看出什么，因此我在图中将一些向量标出来，如下图：

由瞪眼法可以看出，向量的排列依然与闵可夫斯基和相似，但这次凸包B的边在排列后相加时会先取反向；同时，若凸包A无法包含凸包B，所画新凸包C会与凸包A外切，仅首尾两点会在凸包A上，其余点均在凸包A外。因此，我们将上面代码稍加修改，得到可以判断两个凸包包含关系的新代码，如下：

/*
* 此处将所有pair<int, int>替换为pair<double, double>，
* 以便计算新凸包。
*/

#include<bits/stdc++.h>
#define double long double
#define pdd pair<double, double>
#define eps 1e-9

using namespace std;

//使用结构体表示向量方便进行极角排序
//增加整型变量num用于记录边向量所属凸包
struct vec
{
	pdd v;
	int num;
	double ang;

	vec(pdd v)
	{
		this->v = v;
		this->num = 0;
		ang = atan2(v.second, v.first);
	}
	vec(pdd v, int num)
	{
		this->v = v;
		this->num = num;
		ang = atan2(v.second, v.first);
	}

	bool operator < (vec& r)
	{
		return this->ang < r.ang;
	}
};

//直线结构体
struct line
{
	pdd p, v;

	line() {}
	line(pdd p, pdd v)
	{
		this->p = p;
		this->v = v;
	}
};

pdd operator + (pdd& l, pdd& r)
{
	return make_pair(l.first + r.first, l.second + r.second);
}

pdd operator - (pdd& l, pdd& r)
{
	return make_pair(l.first - r.first, l.second - r.second);
}

bool cmp(vec l, vec r)
{
	return l.ang < r.ang;
}

//向量叉乘
double cross(pdd l, pdd r)
{
	return l.first * r.second - l.second * r.first;
}

//求直线交点
pdd getNode(line l, line r)
{
	pdd p2 = r.p + r.v;
	double s1 = cross(l.v, r.p - l.p), s2 = cross(l.v, p2 - l.p);

	return make_pair((s2 * r.p.first - s1 * p2.first) / (s2 - s1), (s2 * r.p.second - s1 * p2.second) / (s2 - s1));
}

int posOfCh(vector<pdd> pochA, int size1, vector<pdd> pochB, int size2) //pochA, pochB分别按顺时针顺序存储两凸包上的点
{
	/*
	* 当两凸包相离时返回0，
	* 相交返回1，
	* 包含返回2
	*/

	/*
	* 此处为上方代码
	*/
	//return flag
	if (!flag)
		return 0;

	vecs.clear();
	pdd minp1 = pochA[0];
	for (int i = 0; i < size1; i++)
	{
		vecs.push_back(vec(pochA[(i + 1) % size1] - pochA[i], 1)); //重新将凸包的边按向量的形式存入，凸包A、B均按逆时针
		minp1 = min(minp1, pochA[i]);	//同时查找凸包A中的最左下角点
	}
	for (int i = 0; i < size2; i++)
		vecs.push_back(vec(pochB[(i + 1) % size2] - pochB[i], 0));
	sort(vecs.begin(), vecs.end()); //极角排序

	//找到第一条属于凸包A的边
	pdd p = minp1;
	int idx = upper_bound(vecs.begin(), vecs.end(), vec({ 0, 1 }), cmp) - vecs.begin() - 1;
	while (!vecs[idx].num)
	{
		p = p - vecs[idx].v;
		idx = (idx + vecs.size() - 1) % vecs.size();
	}

	//计算新凸包上的点
	vector<pdd> ps;
	line las = line(p, vecs[idx].v);
	p = p + vecs[idx].v;
	idx = (idx + vecs.size() - 1) % vecs.size();
	for (int j = 0; j < size1; idx = (idx + vecs.size() - 1) % vecs.size())
	{
		if (vecs[idx].num)
		{
			ps.push_back(getNode(las, line(p, vecs[idx].v)));
			las = line(p, vecs[idx].v);
			p = p + vecs[idx].v;
			j++;
		}
		else
			p = p - vecs[idx].v;
	}

	//当凸包A无法包含凸包B时，新凸包C会与凸包A相切，除首尾两点在凸包A上其余点均在凸包A外，用Pick定理进行判断
	flag = 1;
	for (int i = 0; i < size1; i++)
		if (cross(ps[1] - pochA[i], pochA[(i + 1) % size1] - pochA[i]) < eps)
		{
			flag = 0;
			break;
		}

	if (!flag)
		return 1;

	//去除新凸包中的三点共线情况
	vector<pdd> ponch;
	for (int i = 0; i < ps.size(); i++)
	{
		while (ponch.size() > 2 && cross(ps[i] - ponch.back(), ponch.back() - ponch[ponch.size() - 2]) < eps)
			ponch.pop_back();

		ponch.push_back(ps[i]);
	}
	if (ponch.size() > 2 && cross(ps[0] - ponch.back(), ponch.back() - ponch[ponch.size() - 2]) < eps)
		ponch.pop_back();

	//叉积判断
	flag = 1;
	for (int i = 0; i < ponch.size(); i++)
		if (cross(minp - ponch[i], ponch[(i + 1) % ponch.size()] - ponch[i]) < eps)
		{
			flag = 0;
			break;
		}

	if (flag)
		return 2;
	else
		return 1;
}