计算几何基础–半平面求交

大家好，欢迎来到IT知识分享网。

半平面，平面的一半。半平面求交，就是求得n个半平面的相交的平面。

一开始可能不太好理解。举个例子：

一条线，有两个半平面          两条线，四个半平面，有一个相     三条线，7个半平面，然后有一个相交的区域，这个区域就是半平面

                                                  交，及存在蓝色和绿色的区域         相交的区域。

我是这样理解的。如有错误，希望路过的大佬能够指出。

在程序实现的时候，向量法是一种普遍的表示方法，一般默认向量左侧是所需要的半平面，则可以通过控制向量的方向来选择半平面。

对于求半平面的交，存在一个结论，即：n个半平面的交或者是一个开放的无穷平面，或者是一个封闭的凸多边形（上面例子很明确）。

n个半平面的交有三种基本求法：1.增量法（O(n^2)）。2.分治法（O(nlog(n))）。3.排序增量法（O(nlog(n))）。

简单介绍一下：其中终点介绍一下3.增量排序算法，因为法3可以说是法1的优化，法2不好扩展（个人认为）。

1.增量法：假设已经得到前（n-1）个半平面的交，对于第n个平面，只需用它来“切割”前n-1个半平面交出的多边形。“切割”的过程需要遍历整个多边形，时间复杂度为O(n)，则算法的总时间复杂度是（O(n^2)）。

2.分治法：将n个半平面分成两部分，分别求完交后再将两部分的交合并求交集。整个算法的效率依赖于能够高效地求取两个凸多边形的。利用凸多边形的交可以在O(n)时间内求取的技术，这个算法的时间复杂度为O(nlog(n))。

3.排序增量法：算法类似求取凸包的Graham-Scan算法。先将所有半平面交，按照半平面边界直线与x正半轴的夹角排序，然后扫描所有的半平面。在扫描的过程中使用一个双端队列维护半平面交。排序的复杂度为O(nlog(n))，扫描的过程为O(n)。排序部分为该算法的瓶颈。

排序增量法求半平面交–算法步骤：

Step 1：将所有的半平面按照极角排序，排序过程还要将平行的半平面去重。

Step 2：使用一个双端队列deque，加入极角最小的半平面。

Step 3：扫描过程每次考虑一个新的半平面：

               A：while deque顶端两个半平面的交点在当前半平面外：删除deque顶端的半平面。

               B：while deque底部两个半平面的交点在当前半平面外：删除deque底部的半平面。

               C：将当前半平面加入deque顶端。

Step 4：删除deque两端延伸出的对于半平面：

               A：while deque顶端两个半平面的交点在底部半平面外：删除deque顶端的半平面。

              B：while deque底部两个半平面的交点在顶部半平面外：删除deque底部的半平面。

Step 5：按顺序求取deque中下那个林半平面的交点，得到n个半平面交出的凸多边形。

实现代码：

struct Point{//点的表示 double x,y; Point(double _x=0,double _y=0){ x=_x;y=_y; } friend Point operator + (const Point &a,const Point &b){ return Point(a.x+b.x,a.y+b.y); } friend Point operator - (const Point &a,const Point &b){ return Point(a.x-b.x,a.y-b.y); } friend Point operator == (const Point &a,const Point &b){ return fabs(a.x-b.x)<EPS&&fabs(a.y-b.y)<EPS; } }convex[MAXN]; struct V{//向量的表示 Point start,end; double ang;//角度，[-180,180] V(Point _start=Point(0,0),Point _end=Point(0,0),double _ang=0){ start=_start;end=_end;ang=_ang; } friend V operator + (const V &a,const V &b){ return V(a.start+b.start,a.end+b.end); } friend V operator - (const V &a,const V &b){ return V(a.start-b.start,a.end-b.end); } }l[MAXN],st[MAXN]; struct Triangle{ Point A,B,C; }; int n,ccnt; double DotMul(V a,V b){//点积 a.end=a.end-a.start;b.end=b.end-b.start; return a.end.x*b.end.x+a.end.y*b.end.y; } double CroMul(V a,V b){//叉积 a×b a.end=a.end-a.start;b.end=b.end-b.start; return a.end.x*b.end.y-b.end.x*a.end.y; } int IsLineInter(V l1,V l2){//相交 if(max(l1.start.x,l1.end.x)>=min(l2.start.x,l2.end.x)&& max(l2.start.x,l2.end.x)>=min(l1.start.x,l1.end.x)&& max(l1.start.y,l1.end.y)>=min(l2.start.y,l2.end.y)&& max(l2.start.y,l2.end.y)>=min(l1.start.y,l1.end.y)){ if(CroMul(l2,V(l2.start,l1.start))*CroMul(l2,V(l2.start,l1.end))<=0&& CroMul(l1,V(l1.start,l2.start))*CroMul(l1,V(l1.start,l2.end))<=0){ return 1; } }return 0; } Point LineInterDot(V l1,V l2){//交点 Point p; double S1=CroMul(V(l1.start,l2.end),V(l1.start,l2.start)); double S2=CroMul(V(l1.end,l2.start),V(l1.end,l2.end)); p.x=(l1.start.x*S2+l1.end.x*S1)/(S1+S2); p.y=(l1.start.y*S2+l1.end.y*S1)/(S1+S2); return p; } int JudgeOut(const V &x,const Point &p){//点在线的左侧 return CroMul(V(x.start,p),x)>EPS;//点在左侧返回0,右侧返回1 } int Parellel(const V &x,const V &y){//平行 return fabs(CroMul(x,y))<EPS; } int Cmp(V a,V b){ if(fabs(a.ang-b.ang)<EPS){//角度相同时，不同的边在不同的位置 //此时有两种return方式， //return CorMul(a,V(b.end-a.start))>=0; //左边的边在后面的位置，这样的话，进行计算的时候就可以忽略 相同角度边的影响了 return CroMul(V(b.end-a.start),V(a.end-b.start))>EPS; //左边的边在前面的位置，要进行进行去重判断 。 } return a.ang<b.ang; } double HplaneIntersection(){ int top=1,bot=0; sort(l,l+n,Cmp); int tmp=1; for(int i=1;i<n;++i){ if(l[i].ang-l[i-1].ang>EPS){//去重,如果该边和前面的边平行，则忽略。 l[tmp++]=l[i]; } }n=tmp; st[0]=l[0];st[1]=l[1]; for(int i=2;i<n;++i){ if(Parellel(st[top],st[top-1])||Parellel(st[bot],st[bot+1])) return 0; while(bot<top&&JudgeOut(l[i],LineInterDot(st[top],st[top-1]))) --top; while(bot<top&&JudgeOut(l[i],LineInterDot(st[bot],st[bot+1]))) ++bot; st[++top]=l[i]; } while(bot<top&&JudgeOut(st[bot],LineInterDot(st[top],st[top-1]))) --top; while(bot<top&&JudgeOut(st[top],LineInterDot(st[bot],st[bot+1]))) ++bot; if(top<=bot+1) return 0.00; st[++top]=st[bot]; ccnt=0; for(int i=bot;i<top;++i){ convex[ccnt++]=LineInterDot(st[i],st[i+1]); }double ans=0; convex[ccnt]=convex[0]; for(int i=0;i<ccnt;++i){ ans+=CroMul(V(Point(0,0),convex[i]),V(Point(0,0),convex[i+1])); }return ans/2; }

然后验证一下板子的正确性，进行OJ测试，例题：

这里是关于半平面求交的题目汇总的链接：计算几何基础–半平面求交