半平面指一个平面的一半。比如在二维空间中,一条直线可以把整个平面分为两部分,左边是一个半平面右边是一个半平面。而半平面交就是一堆半平面的交集,可以想象成线性规划方程组
A0x+B0y+C0>=0,A1x+B1y+C1<=0......
的解集
如果一根直线是有向的,我们可以定义它的左边是它对应的半平面(也可以是右边),那么半平面交就是一堆有向直线的左半平面的交集。例如一个多边形可以想象成以逆时针为顺序的直线的半平面交;用这个方法我们可以得到多个多边形的交集面积。
有向直线的定义:
struct Line
{
poi p; vec v;
double ang;//极角
Line(){}
Line(poi p,vec v) :p(p),v(v) {ang = atan2(v.y,v.x);}//极角用斜率的反正切函数求出
bool operator <(const Line &L) const//按极角排序
{ return ang < L.ang; }
}
求半平面交的O(nlogn)算法和凸包类似(下面有更好理解的O(n^2)算法),用一个双端队列来保存决策。先按照极角排序,那么每次新加入一个半平面后会让队尾一些半平面变得无用。由于极角序是环形的,新加入的半平面也可能让队首的半平面变得无用。
代码:
bool Onleft(Line A, poi p) {return cross(A.v, p-A.p) > 0;}
//用叉积判断点是否在直线左边
int HalfplaneIntersection(Line *L,Line *Q,poi *a,poi *out,int cnt)
{
sort(L,L+cnt);//极角排序
int head = 0,tail = 0;
Q[0] = L[0];//双端队列初始化
for(int i = 1; i < cnt; i++)
{
while(head < tail && !Onleft(L[i],a[tail-1])) --tail;//队尾变得无用
while(head < tail && !Onleft(L[i],a[head])) ++head;//队首变得无用
Q[++tail] = L[i];
if(fabs(cross(Q[tail].v,Q[tail-1].v)) < eps){
--tail;
if(Onleft(Q[tail],L[i].p)) Q[tail] = L[i];//新的直线是否在原来直线的左边
}//两向量平行取内侧一个
if(head < tail) a[tail-1] = GetLineIntersection(Q[tail-1],Q[tail]);//求交点
}
while(head < tail && !Onleft(Q[head],a[tail-1])) --tail;//删除无用平面如下图
if(tail - head <= 1) return 0;//空集
a[tail] = GetLineIntersection(Q[head],Q[tail]);//首尾两个平面的交点
int m = 0;
for(int i = head; i <= tail; i++) out[m++] = a[i];//输出
return m;
}
如果你觉得上面算法太难理解,那么这个O(n^2)的算法就好理解多了,即每次用一根直线去切已知的半平面交集(由于最开始的平面为无限大,要手动添加一个很大很大的虚拟平面,去切割这个平面)。而且求多边形的核也可以用到。
多边形的内核:内核是多边形内的一个点集,点集内任意一点与多边形边上任意一点的连线都在多边形内。你可以把多边形想成一个房间,在内核内任意一个点放上一个全方位360度无死角的摄像机,这个摄像机能够看到房间的任意角落。
内核的求法是用半平面交,按照顺时针或者逆时针顺序储存原来多边形。每次取出多边形两个点,以它为直线割去已知半平面的一部分,剩下的就是内核。
代码(半平面交,多边形的核):
struct polygon{
poi a[MAXN];
int sz;
};//多边形
polygon CutPolygon(polygon poly, poi A, poi B)//用有向直线去切割多边形
{
polygon newpoly;
int m = 0,n = poly.sz;
for(int i = 0; i < n; i++)
{
poi C = poly.a[i], D = poly.a[(i+1)%n];//切割如下图
if(dcmp(cross(B-A, C-A)) <= 0) newpoly.a[m++] = C;//顺时针写<=0表示C在AB右边,逆时针写>=0,表示C在AB左边
if(dcmp(cross(B-A, C-D)) != 0) {//直线不与线段重合
poi ip = GetLineIntersection(A, B-A, C, D-C);
if(OnSegment(ip, C, D)) newpoly.a[m++] = ip;//判断ip是不是在线段CD上
}
}
newpoly.sz = m;
return newpoly;
}
double Get_Kernel(polygon poly)//求多边形的核
{
polygon knl;
knl.a[0] = poi(-1e6,-1e6),knl.a[1] = poi(-1e6,1e6),knl.a[2] = poi(1e6,1e6),knl.a[3] = poi(1e6,-1e6);//初始大平面
knl.sz = 4;
for(int i = 0; i < poly.sz; i++)//poly是要求内核的多边形,不停地用它的边去切割平面
{
poi A = poly.a[i], B = poly.a[(i+1)%poly.sz];
knl = CutPolygon(knl,A,B);
}
double area = 0;//求核的面积,可以自己修改
for(int i = 1; i < knl.sz-1; i++)
area += cross(knl.a[i]-knl.a[0], knl.a[i+1]-knl.a[0]);
return area;
}
例题:
POJ1279 http://poj.org/problem?id=1279
求多边形内核面积(用C++编译器,G++貌似有问题)
POJ1474 http://poj.org/problem?id=1474
判断多边形内核是否存在
题解在此