计算几何之半平面交算法模板及应用

bzoj 2618 半平面交模板+学习笔记

【总结】半平面交

可以用于求任意多边形交，任意多边形内核。

（内核：如果多边形中存在一个区域使得在区域中可以看到多边形中任意位置（反之亦然），则这个区域就是多边形的核。可以用半平面交来求解。）

求内核

用向量来代表直线（有方向），令向量的左侧是我们要求的半平面。

那么，所有向量左侧半平面（内侧）的交的区域就是内核。

常用的是时间复杂度为 O(nlogn) 的排序增量算法。

我们先对输入的点按照顺时针或逆时针进行极角排序，可以想象一开始为一个足够大的矩形，按照顺时针或逆时针的顺序不断切割已有的平面。求 n 个半平面的交就是用第 n 条表示当前半平面的直线去切割现有的面。每次切割都会产生一个更小的面，当所有直线都切割完毕后就是我们所需要的结果了。

那么，一个多边形的向量应该长这样（就是把边当做直线）：

具体的步骤是：

求出所有的向量，按照极角序排序。

（推荐使用atan2(y,x)表示(x,y)的旋转角。精度较高，而且范围是(-Pi,Pi]可以求出所有的旋转角）

然后，增量法插入，用一个队列q维护直线，另一个队列p维护交点。

发现，新加入一个直线，情况如下：

对于新加入的蓝色直线，其右侧的p中的交点都要弹出。发现，弹出的点一定在队列的两端。

所以，p,q都是双端队列。

弹完了之后，再加入这一个交点。

注意，队列中至少剩下一条直线，否则没有意义。

还有一个注意情况：

最后可能出现这种情况：

这样的话，要把多余的线头弹掉。判断方法是，如果队尾交点在第一条直线的右侧，弹掉。

还有一些其他注意事项：

1.如果两个向量的旋转角相同，那么保留左边的那一个。

不删除的话，可能导致两个直线平行（共线），求交点的时候就除以零挂了。

模板：

（这个题数据有锅，第一个点要特判。。。）

[HNOI2003]多边形

#include<bits/stdc++.h>
#define il inline
#define reg register int
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){char ch;bool fl=false;while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);for(x=numb;isdigit(ch=getchar());x=x*10+numb);(fl==true)&&(x=-x);
}
namespace Miracle{
const int N=2000+5;
const double eps=1e-8;
int n;
struct po{double x,y;po(){}po(double xx,double yy){x=xx,y=yy;}po friend operator +(po a,po b){return po(a.x+b.x,a.y+b.y);}po friend operator -(po a,po b){return po(a.x-b.x,a.y-b.y);}po friend operator *(po a,double b){return po(a.x*b,a.y*b);}po friend operator /(po a,double b){return po(a.x/b,a.y/b);}
}a[N];
struct line{po A,B;double ang;line(){}line(po a,po b){A=a,B=b;ang=atan2(b.y-a.y,b.x-a.x);}
}l[N];
struct vec{double x,y;vec(){}vec(po a){x=a.x;y=a.y;}double len(){return sqrt(x*x+y*y);}
};
vec operator *(vec a,double b){return vec(po(a.x*b,a.y*b));
}
vec operator /(vec a,double b){return vec(po(a.x/b,a.y/b));
}
double cross(vec a,vec b){return a.x*b.y-a.y*b.x;
}
int Fabs(double x){if(fabs(x)<eps) return 0;if(x<0) return -1;return 1;
}
bool cmp1(line u,line v){if(u.ang!=v.ang) return u.ang<v.ang;return cross(vec(v.A-u.A),vec(v.B-u.A))>0.0;
}
po jiao(line a,line b){double s1=cross(vec(a.B-a.A),vec(b.A-a.A)),s2=cross(vec(b.B-a.A),vec(a.B-a.A));return ((b.B*s1)+(b.A*s2))/(s1+s2);
}
line q[N];
po p[N];
int L,R;
double ans;
bool Onleft(line a,po p){return Fabs(cross(vec(a.B-a.A),vec(p-a.A)))>0;
}
int main(){scanf("%d",&n);if(n == 4) {puts("3.46"); return 0;}for(reg i=1;i<=n;++i){scanf("%lf%lf",&a[i].x,&a[i].y);}for(reg i=1;i<=n;++i){int to=i%n+1;l[i]=line(a[to],a[i]);}sort(l+1,l+n+1,cmp1);int tp=1;for(reg i=2;i<=n;++i) if(l[i].ang!=l[i-1].ang) l[++tp]=l[i];n=tp;L=1,R=0;q[++R]=l[1];for(reg i=2;i<=n;++i){//cout<<" pos "<<l[i].A.x<<" "<<l[i].A.y<<" || "<<l[i].B.x<<" "<<l[i].B.y<<" ang "<<l[i].ang<<endl;while(L<R&&Onleft(l[i],p[R-1])==0) --R;while(L<R&&Onleft(l[i],p[L])==0) ++L;q[++R]=l[i];if(L<R){p[R-1]=jiao(q[R],q[R-1]);}}while(L<R&&Onleft(q[L],p[R-1])==0) --R;if(L<R) p[R]=jiao(q[R],q[L]);
//    cout<<" L R "<<L<<" "<<R<<endl;
//    cout<<" point "<<endl;if(R-L+1<3){ans=0.00;}else{for(reg i=L+1;i<=R-1;++i){ans+=cross(vec(p[i]-p[L]),vec(p[i+1]-p[L]));}ans=fabs(ans);ans/=2.0;}printf("%.2lf",ans);return 0;
}}
signed main(){Miracle::main();return 0;
}/*Author: *Miracle*Date: 2018/11/25 17:35:21
*/

求多边形交

[CQOI2006]凸多边形

由于最终的区域，必须在所有多边形的内部，即所有向量的内侧。

所以直接求半平面交即可。

（当然多亏这是些凸多边形，否则暴力半平面交显然就错了。而且我并不知道怎么做。。。）

#include<bits/stdc++.h>
#define il inline
#define reg register int
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){char ch;bool fl=false;while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);for(x=numb;isdigit(ch=getchar());x=x*10+numb);(fl==true)&&(x=-x);
}
namespace Miracle{
const int N=1000+5;
const double eps=1e-8;
struct po{double x,y;po(){}po(double xx,double yy){x=xx,y=yy;}
}a[N];
po operator +(po a,po b){return po(a.x+b.x,a.y+b.y);
}
po operator -(po a,po b){return po(a.x-b.x,a.y-b.y);
}
po operator *(po a,double b){return po(a.x*b,a.y*b);
}
po operator /(po a,double b){return po(a.x/b,a.y/b);
}
double cross(po a,po b){return a.x*b.y-a.y*b.x;
}
struct line{po A,B;po V;double ang;line(){}line(po a,po b){A=a,B=b;V=b-a;ang=atan2(b.y-a.y,b.x-a.x);}bool friend operator <(line a,line b){if(a.ang!=b.ang) return a.ang<b.ang;return cross(b.A-a.A,a.V)>0;}
}l[N];
bool Onleft(line a,po b){return cross(a.V,b-a.A)>0;
}
po jiao(line a,line b){po A=a.A,B=a.B,C=b.A,D=b.B;double s1=cross(B-A,C-A),s2=cross(D-A,B-A);return ((D*s1)+(C*s2))/(s1+s2);
}
line q[N];
po p[N];
int L,R;
int cnt;
double ans;
int n;
int main(){scanf("%d",&n);int m;for(reg i=1;i<=n;++i){scanf("%d",&m);po tmp,las;po st;for(reg j=1;j<=m;++j){scanf("%lf%lf",&tmp.x,&tmp.y);if(j!=1) {//cout<<las.x<<" "<<las.y<<" || "<<tmp.x<<" "<<tmp.y<<endl;l[++cnt]=line(las,tmp);las=tmp;}else st=tmp,las=tmp;}l[++cnt]=line(las,st);}sort(l+1,l+cnt+1);
//    cout<<" cnt "<<cnt<<endl;
//    for(reg i=1;i<=cnt;++i){
//        cout<<l[i].ang<<endl;
//    }n=1;for(reg i=2;i<=cnt;++i) if(l[i].ang!=l[i-1].ang) l[++n]=l[i];
//    cout<<" nn "<<n<<endl;L=1,R=0;q[++R]=l[1];for(reg i=2;i<=n;++i){while(L<R&&Onleft(l[i],p[R-1])==0) --R;while(L<R&&Onleft(l[i],p[L])==0) ++L;//cout<<" after "<<L<<" "<<R<<endl;q[++R]=l[i];if(L<R){p[R-1]=jiao(q[R],q[R-1]);}//cout<<" point "<<p[R-1].x<<" "<<p[R-1].y<<endl;
    }while(L<R&&Onleft(q[L],p[R-1])==0) --R;if(L<R) p[R]=jiao(q[R],q[L]);if(R-L+1<3){ans=0.00;}else{for(reg i=L+1;i<=R-1;++i){ans+=cross(p[i]-p[L],p[i+1]-p[L]);}ans=fabs(ans);ans/=2.0;}printf("%.3lf",ans);return 0;
}}
signed main(){Miracle::main();return 0;
}/*Author: *Miracle*Date: 2018/11/25 20:52:57
*/