有向直线的定义:
struct Line
{
Point P;//直线上任意一点
Vector v;//方向向量,左边是对应的半平面
double ang;//极角。从x正半轴旋转到向量v所需要的角度(弧度)
Line() {}
Line(Point P, Vector v):P(P),v(v) { ang=atan2(v.y,v.x); }
bool operator < (const Line& L) const { return ang<L.ang; }//极角排序
};
半平面交算法:
bool OnLeft(Line L, Point p) { return Cross(L.v,p-L.P)>0; }
//和凸包类似的算法,半平面可能“绕了一圈”所以队首可能变得无用
//如果可能出现“无界”情况,需要在外面加一个大框(4个半平面)
int HalfplaneIntersection(Line* L, int n, Point* poly)
{
sort(L,L+n);//极角排序
int first,last;
Point *p=new Point[n];//p[i]为q[i]和q[i+1]的交点
Line *q=new Line[n];//双端队列
q[first=last=0]=L[0];//初始时只有一个平面
for (int 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 (!fcmp(Cross(q[last].v,q[last-1].v)))
{
last--;
if (OnLeft(q[last],L[i].P)) q[last]=L[i];
}
if (first<last) p[last-1]=GetIntersection(q[last-1],q[last]);
}
while (first<last && !OnLeft(q[first],p[last-1])) last--; //删除无用平面
if (last-first<=1) return 0;//空集
p[last]=GetIntersection(q[last],q[first]);//计算首尾两个半平面交点
int m=0;//从deque复制到输出
for (int i=first; i<=last; i++) poly[m++]=p[i];
return m;
}