计算几何之半平面交代码模板及过程证明

野性酷女 2022-11-07 00:57 300阅读 0赞

### 半平面交定义 ###

首先要了解半平面交是什么，简单来讲，就是一堆直线，我们只取他的一边，所构成的区域就是半平面交。比如，有4条**有向直线**，我们都只取直线左边的那一部分，所有直线左边平面的交集就是半平面交。如下图所示，蓝色区域就是四条直线左边平面的交集：  
![eg][]

### 半平面交求法 ###

了解了什么是半平面交后，我们看一下他的求法。和[求凸包][Link 1]一样，分为两步。

**第一步：** 将所有直线按倾角从小到大排序  
**第二步：** 从小到大枚举每条边来**维护**一个双端队列。

看着很简单，但是如何**维护每条边**呢？  
如下图，将队尾的边记作b，b的上一条边记作c，新边记作a，如果b与c的交点在直线a的右侧，那就删去队尾的边b，然后再将新边a入队；如果在左侧，就说明队尾的边也是我们要留的边，那就不用出队，直接将新边入队；如果交点就在直线a上，说明交点重合了，这就具体看题目要求删不删重复的边。注意这里是while操作，如果队尾出队了，我们要判断新的队尾满不满足要求，直到满足为止。（和求凸包类似）  
注：以下只是假设的情况  
![eg][eg 1]

我们这里维护的是双端队列，所以每次都要对队尾和队首都判断交点在其左侧还是右侧，因为当我们转一圈过来后，就要判断是否保留队首的直线了，如下图所示：  
![eg][eg 2]  
假设我们是从直线b开始维护的，b就是队首，c是队首的下一条边，a是新的边，新的边也要和队首的边判断，同样的判断：判断b和c的交点是否在a的右侧，是的话就删去队首的边（图中的情况）；不是的话就不删，保留队首。（我们对队尾和队首都要维护，所以用的是双端队列）

记得循环完后，一定要再判断一下队首的边和队尾的两条边的交点，看看边界是否满足半平面交的要求（类比[求凸包][Link 1]），否则可能会出现如下情况：  
![wr][]  
如果只经过刚才队首队尾的维护，队尾的这条红边最后也会入队，就成了如上情况，但明显红边是不对的，所以最后我们还需要判断一下队首和队尾的两条边，来避免这种边界情况。

当我们维护完所有边后，队列中的边就是我们需要的边。但是可能有些人还有一些疑问，如何求直线的交点：[基础知识][Link 2]中有，这里直线都用点向式存。如何判断点在直线的右侧还是左侧？向量的叉积法，也在基础知识中有。如何求直线的倾角？用[atan2()][atan2]这个函数，用atan()求倾角时，若直线是竖直的，则tan不存在，而atan2函数是arctan的一个封装，解决了竖直直线的情况，但要注意他的参数。  
具体可以看看代码来理解每一步。

### 代码模板 ###

模板题：[Acwing 2803：凸多边形][Acwing 2803]  
只需要算所有多边形的所有边的半平面交即可

#include<iostream>
    #include<cmath>
    #include<algorithm>
    using namespace std;
    const int N = 510;
    const double eps = 1e-8;
    typedef pair<double,double> PDD;	//用pair来存点
    struct Line{ 		//定义Line结构体存直线
    	PDD st,ed;	//存直线上的两个点-起点和终点，虽然存的是两个点，但是是点想法表示的直线，起点就是点，终点-起点就是直线的向量
    }line[N];		//line数组存直线
    
    PDD p[N],ans[N];	//p存所有点，ans存最终半平面交的所有交点
    int n,m,cnt,q[N];	//q是用来维护的双端队列，手动建双端队列
    
    int judge(double a,double b)	//判断两个浮点数的大小，相等返回0，小于返回-1，大于返回1
    { 
    	if(fabs(a-b) < eps)
    		return 0;
    	if(a < b)
    		return -1;
    	return 1;
    }
    
    PDD operator-(PDD a,PDD b)	//重载运算符减号，可以直接用pair相减
    { 
    	return { a.first-b.first , a.second-b.second};
    }
    
    double angle(const Line& a)	//求一个直线角度的函数
    { 
    	return atan2(a.ed.second-a.st.second , a.ed.first-a.st.first);	//注意atan2函数的参数分别为纵坐标y和横坐标x，这里xy是直线从原点出发的向量横纵坐标
    }
    
    double cross(PDD a,PDD b,PDD c)		//求向量b-a和c-a叉积的函数，参数是点，方便调用
    { 
    	PDD u,v;
    	u = b - a;	//两向量
    	v = c - a;
    	return u.first*v.second - v.first*u.second;	//叉积 = x1*y2 - x2*y1
    }
    
    bool cmp(const Line& a,const Line& b)		//比较函数，用来对直线按倾角排序
    { 
    	double A = angle(a),B = angle(b);	//求出直线a、b的倾角
    	if(judge(A,B) == 0)		//如果A和B相等的话
    		return cross(a.st,a.ed,b.ed) < 0;	//将直线a、b按从左到右的顺序排，相同角度的直线只需要最左边的就行
    	return A < B;	//按角度从小到大排
    }
    
    PDD intersection(const Line& a,const Line& b)	//求两直线交点的函数，基础知识里都有
    { 
    	PDD p,q,w,v,u;
    	p = a.st,q = b.st,v = a.ed - a.st,w = b.ed - b.st;	//正常的板子求法
    	u = p - q;
    	double t = cross({ 0,0},w,u) / cross({ 0,0},v,w);
    	return { p.first+t*v.first , p.second+t*v.second};	//先将两条直线转换到板子需要的东西，然后按板子求就行，不懂得可以看一下基础知识的板子
    }
    
    bool onright(const Line& a,const Line& b,const Line& c)		//判断b、c的交点是不是在直线a的右边
    { 
    	PDD f = intersection(b,c);	//先求出直线b、c的交点f
    	return cross(a.st,a.ed,f) < 0;	//然后判断交点f是不是在a的右边，叉积小于0说明在右边
    }
    
    double half_plane_intersection()	//半平面交的函数
    { 
    	sort(line,line+cnt,cmp);	//第一步先排序
    	int hh = 0,tt = -1;		//手建头尾
    	for(int i = 0;i < cnt;i++)	//循环每条边
    	{ 
    		if(i && judge(angle(line[i]),angle(line[i-1])) == 0)	//如果当前直线和上一条直线的角度相同就跳过，因为我们已经在cmp中将相同角度的直线按从左到右的顺序排过了，相同角度的直线只取最左边的就行，其余的都没用
    			continue;
    		while(hh+1 <= tt && onright(line[i],line[q[tt-1]],line[q[tt]]))	//维护队尾，当队列元素大于等于2，且队尾两直线的交点在新直线的右边时
    			tt--;	//队尾出队
    		while(hh+1 <= tt && onright(line[i],line[q[hh]],line[q[hh+1]]))	//同理维护队首，注意是while
    			hh++;	//队首出队
    		q[++tt] = i;	//新直线入队
    	}
    	while(hh+1 <= tt && onright(line[q[hh]],line[q[tt-1]],line[q[tt]]))
    		tt--;		//最后判断边界问题的，就是上述的最后一种情况
    	q[++tt] = q[hh];	//将队首入队尾，方便下一行求最开始的第一个交点
    	int k = 0;
    	for(int i = hh;i < tt;i++)	//循环队列中的每一条直线
    		ans[k++] = intersection(line[q[i]],line[q[i+1]]);	//将相邻两条直线的交点求出并记入ans数组，队首和队尾的交点也在其中
    	double area = 0;	//面积
    	for(int i = 1;i < k-1;i++)	//循环ans数组中的所有点
    		area += cross(ans[0],ans[i],ans[i+1]);	//叉积计算面积
    	return area / 2;	//面积除以2
    }
    
    int main()
    { 
    	cin >> n;	//n个凸多边形
    	while(n--)
    	{ 
    		cin >> m;	//m条边
    		for(int i = 0;i < m;i++)
    			cin >> p[i].first >> p[i].second;	//将所有点存进p数组
    		for(int i = 0;i < m;i++)
    			line[cnt++] = { p[i],p[(i+1)%m]};	//将所有边存进line数组
    	}
    	double ans = half_plane_intersection();		//求半平面交
    	printf("%.3lf\n",ans);		//输出面积
    
    	return 0;
    }

--------------------

求半平面交虽然看着比较复杂，函数比较多，但其实理解了每一步后其实就是个套路，要想理解好还得多看。

[eg]: /images/20221023/71fac7bf829942abbc315d5ea53e8c93.png
[Link 1]: https://blog.csdn.net/qq_45735851/article/details/114553142
[eg 1]: /images/20221023/2f1dffd578074fcda0cbb3762e32c824.png
[eg 2]: /images/20221023/0ef167529e414abe98728e0b33623df7.png
[wr]: /images/20221023/c5f1ad2f64d04b7088034a1da4ed3ebb.png
[Link 2]: https://blog.csdn.net/qq_45735851/article/details/114382634
[atan2]: https://baike.baidu.com/item/atan2/10931300?fr=aladdin
[Acwing 2803]: https://www.acwing.com/problem/content/2805/