简单的说,半平面交问题就是给出若干个半平面,求它们的公共部分。其中每个半平面用一条有向直线表示。它的左侧就是它所代表的半平面。
在很多情况下,这个半平面交都是一个凸多边形,但也有时候会得到一个无界多边形,甚至是一条直线、线段或者是点。不管怎么样,结果移动是凸的(因为凸集的交是凸的)。当然,半平面交也可以为空。
计算半平面交的一个方法是增量法,即初始答案为整个平面,然后逐一加入各个半平面,维护当前的半平面交,为了编程方便,我们一般用一个很大的矩形(4个半平面的交)代替整个平面,计算出结果以后删除这4个人工半平面,这样,每加入一个平面都相当于用一条有向直线去切割多边形。
切割的方法很简单:按照逆时针顺序考虑多边形的所有顶点,保留在直线左侧和直线上的点,而删除直线右边的点。如果有向直线和多边形相交产生了新的点,这些点应加在新多边形中。具体来说,没考虑完一个顶点Pi,在考虑Pi+1之前先判断PiPi+1是否和有向直线在PiPpi+1的内部(端点不算)相交。如果是,则还要把交点加入新多边形中。
用有向直线A->B切割多边形poly,返回左侧。如果退化,可能会返回单点或者线段。
Polygon CutPolygon(const Polygon &poly, Point A, Point B)
{
int i, n;
Polygon newpoly;
Point C, D, ip;
n = poly.size();
for(i = 0; i < n; i++)
{
C = poly[i];
D = poly[(i+1)%n];
if(dcmp(Cross(B-A, C-A)) >= 0)
newpoly.push_back(C);
if(dcmp(Cross(B-A, C-D)) != 0)
{
ip = GetLineIntersection(A, B-A, C, D-C);
if(OnSegment(ip, C, D))
newpoly.push_back(ip);
}
}
return newpoly;
}
可惜每次切割需要O(n)时间,因此上述算法的时间复杂度为O(n^2),不是很优秀。
和凸包类似,半平面交也可以通过排序、扫描的方法在O(nlogn)时间内解决,不同的是凸包用的栈,而半平面用的是双端队列。注意,按照极角排序后,每次新加入的半平面可能会让队尾的半平面变得无用,从而需要删除。
注意,由于极角序是环形的,新加的半平面也可能饶了一圈以后让队首的半平面变得无用。
点p在直线L的左边,线上不算。
bool OnLeft(Line L, Point p)
{
return Cross(L.v, p-L.p) > 0;
}
二直线交点,假定交点唯一存在
Point GetIntersection(Line a, Line b)
{
Vector u;
double t;
u = a.p - b.p;
t = Cross(b.v, u) / Cross(a.v, b.v);
return a.p + a.v*t;
}
半平面交的主过程
int HalfplaneIntersection(Line *L, int n, Point *poly)
{
int i, m, first, last;
Point p[maxn]; //p[i]为q[i]和q[i+1]的交点
Line q[maxn]; //双端队列
sort(L, L+n); //按极角排序
first = 0; //双端队列的第一个元素的下标
last = 0; //双端队列的最后一个元素的下标
q[first] = L[0]; //双端队列初始化为只有一个半平面L[0]
for(i = 1; i < n; i++)
{
while(first < last && !OnLeft(L[i], p[last-1])) //新的半平面让队尾变得无用
last--;
while(first < last && !OnLeft(L[i], p[first])) //新的半平面让队首变得无用
first++;
q[++last] = L[i]; //将新的半平面加入队列中
if(dcmp(Cross(q[last].v, q[last-1].v)) == 0)
{
//两向量平行且同向,取内侧的一个
last--;
if(OnLeft(q[last], L[i].p))
q[last] = L[i];
}
if(first < last) //获得交点
p[last-1] = GetIntersection(q[last], q[last-1]);
}
//删除无用平面 首部半平面可能让队尾元素无用
while(first < last && !OnLeft(q[first], p[last-1]))
last--;
if(last - first <= 1) //空集
return 0;
p[last] = GetIntersection(q[first], q[last]); // 计算首尾两个半平面的交点
//从deque复制到输出中
m = 0;
for(i = first; i <= last; i++)
poly[m++] = p[i];
return m;
}
在大多数情况下,若干半平面的交是一个凸多边形,但也有例外,比如可能得到的是一个无界的区域,解决方法和前面一样,在外面加一个坐标很大(注意不要让运算溢出)的包围框(由4个半平面组成),最后再把框去掉。但对于保证不会产生无界区域的情况下,就不需要加这4个特殊的半平面了。