【笔记篇】最良心的计算几何学习笔记（六）
半平面交
github传送门
简介
Emmmm学完旋转卡壳感觉自己已经是个废人了..
修整了一个周末, 回来接着跟计算几何势力硬干... (这个周末是不是有点长?)
今天就讲讲半平面交吧.
请自己回顾必修五 线性规划相关知识...
什么是半平面?
就是一条直线一侧的点构成的集合..
用解析几何的观点来看就是满足\(Ax+By+C<0\)这个不等式的点的集合.
那么半平面交自然就是许多这样的集合的交集咯~
最后就很像线性规划做出来的可行域一样...
半平面交可以长这样
这样
甚至这样
也就是说, 半平面交的结果可能是无界的一大片, 可能是一个有界的凸多边形, 也可能是点、线段、直线(没图╮(╯_╰)╭), 甚至空集..
引例
题目
先来道不那么标准的板子题:戳这里..
分析
这个很显然就是要求朝上的无界的一片的那种.. md就是个下凸壳啊...
所以我们按照凸壳的方式试着思考一下.
我们先将直线按斜率排序(都按斜截式给你了显然不会不存在....) 如果斜率一样截距小的可以直接不要了(肯定会被挡住..)
维护一个栈, 先将第一条直线压入, 这条直线肯定能露出来(谁让斜率最大了呢)
然后我们考虑新加入的直线对前面的直线的影响.
我们可以看到,当新加入一条斜率比较小的直线(绿色或黄色, 按理说应该是平行的, 如果看着不大平行就凑合看吧)
- 如果比前面两条直线的交点要高(比如绿色) 那么前一条直线就会被遮住(红色左边被绿遮, 右边被黑遮...);
- 而如果比前面两条直线的交点低(比如黄色) 那么前一条直线就会有一部分露出来.
所以出现了绿色的情况我们就把红色退栈. 然后继续判断.
是不跟Graham扫描法一个样纸?
时间复杂度\(O(nlogn)\) 复杂度仍然卡在排序上. 这是整数而且斜率极差不超过1e6所以可以很方便的用基数排序降到\(O(n)\)虐场.. 然而窝还是用sort 因为省事...
现在唯一的问题就是判断直线在点上面了(根据样例可以得知交点一个点是不算子线段的, 也要算在上面).. 但这就是解析几何的问题了.. (其实可以用类似于方向向量的方法解, 只是我不想再写叉积了(想起被叉积支配的恐惧..))
我们令前两条直线分别为\(y=k_1x+b_1, y=k_2x+b_2\), 新加入的直线为\(=k_3x+b_3\)
联立得到前两条直线交点横坐标\(x=\frac{b_2-b_1}{k_1-k_2}\)
那么根据代入法我们可以得到
其实这样已经可以做了, 而且题目坐标都是整数精度问题也不大, 但是我还是想找一个优美的形式...
由于坐标排过序, 所以\(k_1-k_2>0\), 我们都移到左边, 然后两边同时乘\(k_1-k_2\)(乘正数不等号不变号), 就能得到
这样就没有分式, 显得更为有优美一些(其实是闲的(卡精度卡怕了))
代码
下面贴上1A代码:
#include <cmath>
#include <cstdio>
#include <algorithm>
const double eps=1e-9;
inline int dcmp(const double &a){
if(fabs(a)<eps) return 0;
return a<0?-1:1;
}
struct line{
double k,b;
int id;
}p[50050],stk[50050]; int tp,n;
inline bool cmp(const line &a,const line &b){
if(!dcmp(a.k-b.k)) return a.b>b.b;
return a.k>b.k;
}
inline bool cmp2(const line &a,const line &b){
return a.id<b.id;
}
void halfPlanesIntersection(){
stk[++tp]=p[0];
for(int i=1;i<n;++i){
if(!dcmp(p[i].k-stk[tp].k)) continue;
while(tp>1&&dcmp((stk[tp].b-stk[tp-1].b)*(stk[tp-1].k-p[i].k)+(stk[tp-1].k-stk[tp].k)*(stk[tp-1].b-p[i].b))<1) --tp;
stk[++tp]=p[i];
}
}
int main(){
scanf("%d",&n);
for(int i=0;i<n;++i){
scanf("%lf%lf",&p[i].k,&p[i].b);
p[i].id=i+1;
}
std::sort(p, p+n, cmp);
halfPlanesIntersection();
std::sort(stk+1, stk+tp+1, cmp2);
for(int i=1;i<=tp;++i)
printf("%d ",stk[i].id);
}
Emmmm这个题还有一种SAO操作就是Trinkle大神提出的半平面交对偶转凸包做法... orz
大体上就是把直线\(y=kx+b\)视为一个点\((k,b)\), 然后求一波上凸壳.
等下个题看下能不能用上吧.
正题
刚才那个题太解析几何了, 拉个平时数学课不睡觉的人就能做..
我们来看更"计算几何"的东西..
比如半平面除了能用线性规划里面的\(Ax+By+C<0\)以外, 还可以用向量的左边来表示...
向量的左边嘛, 当然就是叉积<0咯(被叉积支配的恐惧...)
方法一
有一种非常简单的\(O(n^2)\)做法, 我们假设这个东西是有界的, 于是就从无穷远处套个框..
然后用每一条边去切割这个多边形就行了.
我们枚举一条边的时候, 分别判断其他的点是否在左边.. 如果在左边显然就不在半平面交了.
然后有一个端点在左边的边会与正在枚举的边有交点(8), 我们就把这个新的交点加入我们的程序即可.
一共要枚举\(n\)条边, 每条边要判断\(n\)次, 所以时间复杂度显然是\(O(n^2)\)的.
这种做法有个非常大的好处就是好理解....
代码(无)
等用到了再写..
方法二
还有一种方法叫分治.. 算法的流程大约是这样的:
-
在最外面套个框
-
递归搞出前一半的半平面交
-
递归高出后一半的半平面交
-
利用凸多边形\(O(n)\)求交(好像是CGI什么的) 来处理整个的半平面交
时间复杂度也好分析: \(T(n)=2T(\frac n2)+O(n)=>O(nlogn)\)...
但是好像比较麻烦(而且凸多边形求交真的不是搞一个半平面交???)
有更好用的\(O(nlogn)\)做法所以这个就不常见了...
代码(无)
代码应该也不太好写 等用到了再写...
方法三
隆重请出我们的方法三...
好像是某ZZY大神自己yy出来的?? orz
然后被某大神优化了一下就出现了现在的版本...
真的是高端大气上档次...
不过这样对于我这种蒟蒻来讲就不大好理解了..
求半平面交的方式和凸包很像很像..
先将直线的方向向量起点都平移到原点, 然后按极角排序..
但是直接用叉积就会出现一些问题(就是不喜欢atan2..), 比如红色不一定能排在紫色的左边了...(为了方便向量的先后关系应该如彩虹的颜色顺序..)
所以我们要分类讨论(几何常见套路), 将向量分成x轴上方(包括x轴正半轴)和x轴下方(包括x轴负半轴)两类. 很显然第一类<第二类(因为要写sort里的cmp函数就用大于小于了).
然后对于每一类分别求叉积就好了.. 比如橙色在黄色的顺时针方向所以橙色<黄色, 比如蓝色在紫色的顺时针方向所以蓝色<紫色, 以此类推.
那么极角相同呢?
由于我们要求的是交集, 所以要把右边的排除. 我们可以通过\(\vec a\)的终点-\(\vec b\)的起点与\(\vec b\)做叉积, 来判断直线的左右关系.
排完序之后就进入了正题.
我们还是看添加一条直线之后会造成什么影响.
(P.S. 图里没有画箭头, 方向都是↗)
根据上面不那么标准的题目可以分析得到, 当我们新加入一条边的时候, 便检测一下前面两条边的交点在不在这条新加入的边所对应的半平面上. 所以对于红点来说, 它在粉色对应的半平面但不在紫色对应的半平面, 所以如果加入的是粉色那就加入了, 而如果加入的是紫色那么浅绿是要退栈的..
而我们的半平面交不只是求一个凸壳(所以说引例不标准嘛), 所以图形可能是封闭的, 那么新加入的边也就可能影响到最开始加入的边...
所以我们要维护一个双端队列.. 每次分别判断队首队尾两条边的交点在不在对应的半平面内. 如果不在就弹出去.
然后就出现了下一个问题, 就是求完以后, 可能会有很多冗余的直线.
比如我们求完黑色的半平面交以后, 红色的也能成功的入队, 但是对答案并没有什么贡献. 队首也是一样, 所以我们要再扫, 把这些冗余的出队.
由于每条直也最多只会进队一次, 出队一次, 所以复杂度是\(O(n)\)的,
总复杂度就自然是排序的\(O(nlogn)\)(所以也可以用基数排序做到\(O(n)\), 但我还是用sort..)
然后细节上来说还是比较麻烦的... 用到了好多前面的基础知识, 然后发现自己全tm不记得了..
代码(这次有了XD)
贴一道算是模板题的题目吧 bzoj2618
把多边形的每一条边拆成直线, 最后求一下凸多边形的面积就行了(其实有点综合的意思~)
WA竟然是因为求面积的时候\(pcnt\)打成了\(n\)... 是不该去打星际了..
代码:
#include <cmath>
#include <cstdio>
#include <algorithm>
const int N=101010;
const double eps=1e-9;
inline int dcmp(const double &a){
if(fabs(a)<eps) return 0;
return a<0?-1:1;
}
struct vec{
double x,y;
vec(double X=0,double Y=0):x(X),y(Y){}
}pt[N],yM; int pcnt;
inline vec operator-(const vec &a,const vec &b){
return vec(a.x-b.x,a.y-b.y);
}
inline double operator*(const vec &a,const vec &b){
return a.x*b.y-a.y*b.x;
}
struct line{
vec a,b;
line(){}
line(vec X,vec Y):a(X),b(Y){}
}p[N],dq[N]; int tp,bt,n; //双端队列
inline bool aboveX(const line &a) //是否在x轴上方
{
if(!dcmp(a.b.y-a.a.y)) return dcmp(a.b.x-a.a.x)>0;
return dcmp(a.b.y-a.a.y)>0;
}
inline bool cmpa(const line &a,const line &b){
if(aboveX(a)!=aboveX(b)) return aboveX(a);
if(!dcmp((a.b-a.a)*(b.b-b.a)))
return dcmp((a.b-a.a)*(b.b-a.a))<0;
return dcmp((a.b-a.a)*(b.b-b.a))>0;
}
inline vec getLinesIntersect(const line &a,const line &b){
double a1,a2,b1,b2,c1,c2,d;
a1=a.b.y-a.a.y; b1=a.a.x-a.b.x; c1=a.a*a.b;
a2=b.b.y-b.a.y; b2=b.a.x-b.b.x; c2=b.a*b.b;
d=a1*b2-a2*b1;
return vec((b2*c1-b1*c2)/d,(a1*c2-a2*c1)/d);
}
inline bool pd(const line &a,const line &b,const line &c){
//判断直线a和直线b的交点是否在c左边(在半平面内)
vec p=getLinesIntersect(a,b);
return dcmp((p-c.a)*(c.b-c.a))>-1; //true表示不在半平面内
}
void halfPlaneIntersection(){
int i,j; vec t;
for(i=0,j=0;i<n;++i)
if(dcmp((p[i].b-p[i].a)*(p[j].b-p[j].a)))
p[++j]=p[i];
n=j+1; //去掉极角相同
dq[bt]=p[0];
dq[++tp]=p[1]; //开始两条直线入队
for(int i=2;i<n;++i){
while(tp>bt&&pd(dq[tp],dq[tp-1],p[i])) --tp;
while(bt<tp&&pd(dq[bt],dq[bt+1],p[i])) ++bt;
dq[++tp]=p[i];
}
while(tp>bt&&pd(dq[tp],dq[tp-1],dq[bt])) --tp;
while(bt<tp&&pd(dq[bt],dq[bt+1],dq[tp])) ++bt;
dq[++tp]=dq[bt]; //已经保存好绕一圈的直线了
//确定半平面交部分的顶点
for(i=bt;i<tp;++i,++pcnt)
pt[pcnt]=getLinesIntersect(dq[i],dq[i+1]);
}
double polyArea(double s=0){
if(pcnt<3) return 0;
pt[pcnt]=pt[0];
for(int i=0;i<pcnt;++i)
s=s+pt[i]*pt[i+1];
return 0.5*fabs(s);
}
int main(){
int T,w=0;
scanf("%d",&T);
while(T--){
int ps; scanf("%d",&ps);
for(int i=1;i<=ps;++i) scanf("%lf%lf",&pt[i].x,&pt[i].y);
pt[0]=pt[ps];
for(int i=0;i<ps;++i) p[n].a=pt[i],p[n++].b=pt[i+1];
}
std::sort(p,p+n,cmpa);
halfPlaneIntersection();
double area=polyArea();
printf("%.3lf",area);
}
其实想清楚的话也不是很难, 但是真的细啊OvO
大约就这样吧?
