初学计算几何（五）——半平面交

前言

学了这么久计算几何，其实就是为了学半平面交。

而学半平面交，其实就是为了做洛谷试炼场里的这道题：【洛谷3222】[HNOI2012] 射箭。

概念

• 半平面：一条直线会把一个平面分成两个半平面，通常我们把直线左边的半平面视为我们所需要的半平面。
• 半平面交：即若干半平面的交集，显然它有以下几种可能性：点/线段/直线/凸多边形/无穷平面。为了避免产生无穷平面，我们一般会加入四个半平面作为限制。

那么半平面交有什么用？

比较常见的作用是求多边形的核（如果多边形中存在一个区域使得在该区域中可以看到多边形中任意位置，则这个区域就是多边形的核）。

大致思路

首先，我们将直线进行极角排序，其中极角相同的直线我们让较左边的直线在前面，这一函数如下：

inline bool operator < (Line A,Line B) {return dcmp(A.Ang-B.Ang)?!~dcmp(A.Ang-B.Ang):!~dcmp(Cro(A.B-A.A,B.B-A.A));}

接下来，我们用一个双向队列\(ql\)来维护边，同时开一个数组\(qp\)，用\(qp_i\)存储\(ql_i\)与\(ql_{i+1}\)的交点。

则每当要加入一条新边，我们就要维护一个凸多边形，大致操作如下：

if(!dcmp(L[i].Ang-L[i-1].Ang)) continue;//如果出现极角相同的情况，将后出现的边跳过
while(H<T&&IsOnRight(qp[T-1],L[i])) --T;while(H<T&&IsOnRight(qp[H],L[i])) ++H;//每次将首尾在当前边右侧的点弹出 
ql[++T]=L[i],qp[T-1]=LineIntersection(ql[T-1],ql[T]);//将边加入队列中，同时更新qp[T-1]

注意，其中的\(IsOnRight()\)函数利用了叉积性质，具体实现如下：

inline bool IsOnRight(Point P,Line L) {return !~dcmp(Cro(L.B-L.A,P-L.A));}//判断一个点是否在一条直线右侧

这样操作一遍之后，前后其实还是有一些多余的点，因此我们需要用首尾两条边再去判断一次，将在这两条边右边的点弹出：

while(H<T&&IsOnRight(qp[T-1],ql[H])) --T;while(H<T&&IsOnRight(qp[H],ql[T])) ++H;

最后判断，如果队列中只剩下小于等于\(2\)个点，就说明面积为\(0\)。

否则，可以把剩下点看成一个凸包，求凸包面积即可。

代码（板子题：【POJ2451】Uyuw's Concert）

#include<cstdio>
#include<algorithm>
#include<cmath>
#define N 20000
using namespace std;
int n;
namespace ComputationGeometry//计算几何
{
	#define eps 1e-10
	inline int dcmp(double x) {return fabs(x)<eps?0:(x>0?1:-1);}
	struct Point
	{
		double x,y;
		Point(double nx=0,double ny=0):x(nx),y(ny){}
	};
	typedef Point Vector;
	inline Vector operator + (Point A,Point B) {return Vector(A.x+B.x,A.y+B.y);}
	inline Vector operator - (Point A,Point B) {return Vector(A.x-B.x,A.y-B.y);}
	inline Vector operator * (Vector A,double x) {return Vector(A.x*x,A.y*x);}
	inline Vector operator / (Vector A,double x) {return Vector(A.x/x,A.y/x);} 
	inline bool operator < (Vector A,Vector B) {return fabs(A.x-B.x)>eps?A.x<B.x:A.y<B.y;}
	inline double Cro(Vector A,Vector B) {return A.x*B.y-A.y*B.x;}
	inline double Angle(Vector A) {return atan2(A.y,A.x);}
	struct Line
	{
		Point A,B;double Ang;
		Line(Point x=Point(0,0),Point y=Point(0,0)):A(x),B(y),Ang(Angle(y-x)){}
	};
	typedef Line Segment;
	inline bool operator < (Line A,Line B) {return dcmp(A.Ang-B.Ang)?!~dcmp(A.Ang-B.Ang):!~dcmp(Cro(A.B-A.A,B.B-A.A));}
	inline Point LineIntersection(Line L1,Line L2) {return L1.A+(L1.B-L1.A)*Cro(L2.B-L2.A,L1.A-L2.A)/Cro(L1.B-L1.A,L2.B-L2.A);}
	inline bool IsOnRight(Point P,Line L) {return !~dcmp(Cro(L.B-L.A,P-L.A));}
	struct Polygon
	{
		int n;Point p[N+5];
		Polygon() {n=0;}
	};
	typedef Polygon ConvexHull;
	inline double PolygonArea(Polygon P)
	{
		register int i;register double res=0;
		for(i=2;i<P.n;++i) res+=Cro(P.p[i]-P.p[1],P.p[i+1]-P.p[1])/2;
		return res;
	}
	inline double HalfPlaneIntersection(int tot,Line *L)//半平面交
	{
		register int i,H=1,T=1;static Line ql[N+5];static Point qp[N+5];sort(L+1,L+tot+1),ql[1]=L[1];
		for(i=2;i<=tot;++i)
		{
			if(!dcmp(L[i].Ang-L[i-1].Ang)) continue;
			while(H<T&&IsOnRight(qp[T-1],L[i])) --T;while(H<T&&IsOnRight(qp[H],L[i])) ++H; 
			ql[++T]=L[i],qp[T-1]=LineIntersection(ql[T-1],ql[T]);
		}
		while(H<T&&IsOnRight(qp[T-1],ql[H])) --T;while(H<T&&IsOnRight(qp[H],ql[T])) ++H;
		if(T-H<=1) return 0;
		register ConvexHull C;
		for(qp[T]=LineIntersection(ql[H],ql[T]),i=H;i<=T;++i) C.p[++C.n]=qp[i];
		return PolygonArea(C);
	}
}; 
using namespace ComputationGeometry;
Line L[N+5];
int main()
{
	register int i;register double X1,Y1,X2,Y2;
	for(scanf("%d",&n),i=1;i<=n;++i) scanf("%lf%lf%lf%lf",&X1,&Y1,&X2,&Y2),L[i]=Line(Point(X1,Y1),Point(X2,Y2));//读入直线
	L[n+1]=Line(Point(0,0),Point(10000,0)),L[n+2]=Line(Point(10000,0),Point(10000,10000)),L[n+3]=Line(Point(10000,10000),Point(0,10000)),L[n+4]=Line(Point(0,10000),Point(0,0));//加入四个半平面作为限制
	return printf("%.1lf",HalfPlaneIntersection(n+4,L)),0;//求解答案
}