【学习小记】半平面交——排序增量法

Preface

之前的半平面交的算法是基于分治和凸包合并的，分治两边，计算出半平面交，再合并凸包。
而这种排序增量法好写简洁常数小，适合在比赛中使用。

Text

为了避免半平面交区域*的情况，我们在无穷远处四个方向加上四个半平面的限制。

可以看出，有限的半平面交是一个凸包

方便起见，我们用点+向量的形式来表示一个半平面，向量的左手向就是半平面的方向。
定义半平面的极角为向量的极角，我们将半平面按照极角排序，对于极角相同的，我们只保留最靠里，也就是限制最紧的那一条。

我们维护一个队列，像做凸壳一样的从队尾加入半平面，如果当前队尾半平面边界和队尾前一个的边界交点在新的半平面之外，说明队尾就没用了，将队尾弹掉。

但是半平面交既然是一个凸包，它会绕一圈回来，如果只在队尾加肯定是不对的。
弹完队尾之后，我们再判断队头的半平面是否也没用了，判定方法类似，如果能的话还要将队头的半平面给弹掉。

这样做完以后，我们发现最后队尾会剩下一些半平面，虽然这些半平面是凸的，但是它们有可能会被队头给弹掉，此时还要再判一下，将一些队尾弹掉。

分析复杂度，每一个半平面最多只会进队出队各一次，时间复杂度仅在于排序，为$O(n\log n)$O(nlogn)

PS：尽量避免在排序重载compare函数的时候用atan2去比较，每次比较都要计算atan2是很慢的，我们可以提前计算好再比，或者用叉积+象限来比。

Code

丑到极致了…

namespace geometry
{
	const double eps=1e-10;
	const double pi=acos(-1);
	struct node
	{	
		double x,y;
		node(double _x=0,double _y=0){x=_x,y=_y;}
	};
	node operator +(node x,node y) {return node(x.x+y.x,x.y+y.y);}
	node operator -(node x,node y) {return node(x.x-y.x,x.y-y.y);}
	node operator *(double x,node y) {return node(x*y.x,x*y.y);}
	double crs(node x,node y) {return x.x*y.y-x.y*y.x;}
	double ang(node x) 
	{
		double v=atan2(x.y,x.x);
		if(v<0) v+=2*pi;
		return v;
	}

	struct line
	{
		node p,v;
		line(){};
		line(node _p,node _v) {p=_p,v=_v;}
	};
	node dot(line x,line y) {return y.p+(crs(x.v,x.p-y.p)/crs(x.v,y.v))*y.v;}
	bool inside(line x,node y)
	{
		return (crs(x.v,y-x.p)>eps);
	}
	bool in(line x,node y)
	{
		return abs(crs(x.v,y-x.p)<=eps);
	}
}
using namespace geometry;

int t,n,n1,n2;
node a[N],ab[N];
line d[4*N],db[4*N];
struct px
{
	double v;
	int p;
	friend bool operator <(px x,px y) {return x.v<y.v;}
}u2[N];


int pre(int x) {return (x==2)?n:x-1;}
int suf(int x) {return (x==n)?2:x+1;}

bool cmp2(px x,px y)
{
	return(x.v+eps<y.v||(abs(x.v-y.v)<=eps&&inside(d[y.p],d[x.p].p)));
}
double get()
{
	fo(i,1,n1) 
	{
		if(in(d[i],node(0,0))) return 0;
		u2[i]=(px){ang(d[i].v),i};
		db[i]=d[i];
	}
	static int dq[4*N],d2[4*N];
	n2=0;
	sort(u2+1,u2+n1+1,cmp2);

	fo(i,1,n1) d[i]=db[u2[i].p];
	
	fo(i,1,n1) if(i==1||fabs(ang(d[i].v)-ang(d[i-1].v))>eps) d2[++n2]=i;
	int l=1,r=2;dq[1]=d2[1],dq[2]=d2[2];
	fo(i1,3,n2)
	{
		int i=d2[i1];
		while(l<r&&!inside(d[i],dot(d[dq[r]],d[dq[r-1]]))) r--;
		while(l<r&&!inside(d[i],dot(d[dq[l]],d[dq[l+1]]))) l++;
		dq[++r]=i;
	}
	while(l+2<r&&!inside(d[dq[l]],dot(d[dq[r]],d[dq[r-1]]))) r--;
	double s=0;//计算面积
	fo(i,l,r-2) s+=crs(dot(d[dq[i]],d[dq[i+1]]),dot(d[dq[i+1]],d[dq[i+2]]));
	s+=crs(dot(d[dq[r-1]],d[dq[r]]),dot(d[dq[r]],d[dq[l]]));
	s+=crs(dot(d[dq[r]],d[dq[l]]),dot(d[dq[l]],d[dq[l+1]]));
	return s/2;
}