【笔记篇】最良心的计算几何学习笔记（七）
动态凸包
本文的github传送门在这里~
======================================================================
不会凸包的赶紧去学一下哦~
======================================================================
好的我们已经会求凸包了. 那我们来看这样一道题.
题目大意(英文题必须要有的翻译过程OvO):
写一个程序支持一下两种操作:
- 1 x y 添加一个点
(x,y) - 2 x y 询问点
(x,y) 是否在当前点集的凸包中
这就不是很好搞了吧?
如果每加一个点就求一遍凸包那么复杂度必然上天.
所以我们就要动态维护.
我们来考虑一下加入一个点的时候会对之前求好的凸包造成什么影响.
很容易从图上看出:
- Case 1. 插入的点在凸包内, 不对凸包造成影响. (蓝)
- Case 2. 凸包中其他点不动, 将新点插入. (红)
- Case 3. 插入新点使得原凸包上的一些点不再在凸包上. (绿)
而考虑到我们的Graham扫描法, 我们可以知道这个就与新加入的点与两侧的边的凹凸有关.
所以我们要找到新点的位置, 然后往两边扫. 因为只有插入, 所以点被删了就不会再回来, 复杂度显然就是
那么如何找新点的位置呢? 我们在Graham扫描法中用到了极角排序. 那么这里我们也可以用极角来搞.
我们只需要用动态的数据结构维护即可. 那我们就想到了平衡树.
但是这里的极角序还是和Graham扫描不太一样… 因为Graham中可以找到
我们要用做半平面交时候用的分
那么基准点就选最开始的凸包上的点就好了. 因为这个点一直都会在凸包里.
由于题目中说开始一定会给一个不退化的三角形, 我们用重心就好啦.
这样我们加入一个点的时候就:
- 插入这个点
- 找到这个点的前驱和后继
- 分别往前和往后判断是否满足凸包性质(就是叉积嘛), 不是就删点.
时间复杂度是查找
而查询的时候就找到前驱和后继然后求一下叉积看看是不是在里面就完了, 时间复杂度也是查找的
但是非常不想码splay….
然后惊奇的发现这个题的操作只有插入、删除、前驱、后继, 那不就可以用set水过去了么= =
于是就愉快地码set…. (其实并不怎么愉快, 因为STL还是有些地方挺反人类的)但是码力太弱set版都写不对还调了好久…. 显然是药丸.
要注意的地方就是我们要手动把set搞成环, 就是把end放在begin的前面, 因为极角排序本身是绕圈圈的.
大约就这样吧.
代码:
#include <set>
#include <cmath>
#include <cstdio>
using namespace std;
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){}
}p[4],sq;
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.x+a.y*b.y;
}
inline double operator*(const vec&a,const vec&b){
return a.x*b.y-a.y*b.x;
}
inline double len(const vec&a){
return sqrt(a^a);
}
inline bool aboveX(const vec&a){
int d=dcmp(a*vec(1,0));
if(!d) return dcmp(a*vec(0,1))>0;
return d<0;
}
inline bool operator<(const vec&A,const vec&B){
vec a=A-sq,b=B-sq;
if(!dcmp(len(a))) return 1;
bool x1=aboveX(a),x2=aboveX(b);
if(x1!=x2) return x1;
int d=dcmp(a*b);
if(!d) return dcmp(len(b)-len(a))>0;
return d>0;
}
inline int gn(int a=0,char c=0,int f=1){
for(;(c<48||c>57)&&c!='-';c=getchar());if(c=='-')f=-1,c=getchar();
for(;c>='0'&&c<='9';c=getchar()) a=a*10+c-'0'; return a*f;
}
typedef set<vec> pset;
typedef pset::iterator pt;
pset s;
pt pre(pt x){
if(x==s.begin()) x=s.end();
return --x;
}
pt suc(pt x){
++x;
if(x==s.end()) return s.begin();
return x;
}
inline bool query(const vec&a){
pt it=s.lower_bound(a);
if(it==s.end()) it=s.begin();
int d=dcmp((a-*it)*(*pre(it)-*it));
if(!d) return dcmp(len(a-*it)-len(*pre(it)-*it))<0;
return d>0;
}
inline void inser(const vec&a){
if(query(a)) return;
pt it=s.lower_bound(a);
if(it==s.end()) it=s.begin();
s.insert(a);
while(s.size()>3&&dcmp((*it-a)*(*suc(it)-*it))<1){
s.erase(it); it=suc(s.find(a));
}
it=pre(s.find(a));
while(s.size()>3&&dcmp((*it-a)*(*pre(it)-*it))>-1){
s.erase(it); it=pre(s.find((a)));
}
}
int main(){
int n=gn(),m;
for(int i=1;i<=3;++i){
m=gn(),p[i].x=gn(),p[i].y=gn();
sq=sq-p[i];
}
sq=vec(-sq.x/3,-sq.y/3);
for(int i=1;i<=3;++i)
s.insert(p[i]);
for(int i=4;i<=n;++i){
int t=gn(),x=gn(),y=gn();
if(t==1) inser(vec(x,y));
else puts(query(vec(x,y))?"YES":"NO");
}
}你看一整道题码完都比一个splay板子短= =
非常感谢写STL的大牛们…
然后我们就可以去水题了.
当然一般的题目是没有这么裸的, 一般要维护一些东西, 比如周长啊 面积啊什么的..
你看这道题就是这样.
咦? 这题不是删除么= =
删除是不好做的, 删掉的点如果在凸壳上就要重新求一遍, 复杂度就上天了, 但是我们可以把操作倒着做, 视为不断的加点和查询凸壳长…
这道题目的出题人非常的良心, 只需要维护一个上凸壳, 而且还做了各种各样的保证, 这样就不用判边界了.. 但是我还是WA了好多次, 原因也很扯淡我们一会再说.
这样的话我们就考虑周长的维护.
很明显地, 加入一个点势必会导致打叉的边被去除, 因为至少有一个左边的点和一个右边的点分别要连向新点. 要删掉的点一定与要删掉的边是一一对应的关系. 这个从图上也能很显然地看出来. 这样我们就可以复杂度不变的维护周长, 最后记得加上两条蓝色边的长度即可.
由于这个题不需要判x轴上下而且特殊情况比较少, 所以写起来会比较舒服.
起码在我交之前是有这么个感觉的.
#include <set>
#include <cmath>
#include <cstdio>
using namespace std;
const int N=101010;
const double eps=1e-9;
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){}
}sq,p[N];
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.x+a.y*b.y;}
inline double operator *(const vec &a,const vec &b){return a.x*b.y-a.y*b.x;}
inline double len(const vec &a){return sqrt(a^a);}
inline bool operator <(const vec &A,const vec &B){
if(A.x==B.x&&A.y==B.y) return 0;
vec a=A-sq,b=B-sq;
int d=dcmp(a*b);
if(!d){
int dd=dcmp(len(a)-len(b));
if(!dd) return dcmp(a*vec(0,1))>0;
return dd>0;
} return d>0;
}
inline int gn(int a=0,char c=0){
for(;c<'0'||c>'9';c=getchar());
for(;c>47&&c<58;c=getchar())a=a*10+c-48;return a;
}
set<vec> s;
typedef set<vec>::iterator pt;
int qu[N<<1]; double ans[N<<1],cur;
bool del[N];
pt pre(pt x){return --x;}
pt suc(pt x){return ++x;}
void inser(vec x){
pt it=s.lower_bound(x);
if(dcmp((x-*it)*(*pre(it)-*it))>-1) return;
cur-=len(*it-*pre(it));
while(suc(it)!=s.end()&&dcmp((x-*it)*(*it-*suc(it)))<1){
cur-=len(*it-*suc(it));
it=suc(it); s.erase(pre(it));
} it=pre(s.lower_bound(x));
while(it!=s.begin()&&dcmp((x-*it)*(*it-*pre(it)))>-1){
cur-=len(*it-*pre(it));
it=pre(it); s.erase(suc(it));
} s.insert(x); it=s.find(x);
cur+=len(*it-*pre(it))+len(*it-*suc(it));
}
int main(){
int n=gn(),x=gn(),y=gn(),m=gn(); sq=vec(n*0.5,0);
vec p1=vec(0,0),p2=vec(n,0),p3=vec(x,y); cur=len(p3-p2)+len(p3-p1);
s.insert(p1); s.insert(p2); s.insert(p3);
for(int i=1;i<=m;++i) p[i].x=gn(),p[i].y=gn(); int q=gn(),tot=0;
for(int i=1;i<=q;++i){
int t=gn();
if(t==1){
int w=gn();
del[w]=1; qu[i]=w;
}
else qu[i]=0;
}
for(int i=1;i<=m;++i) if(!del[i]) inser(p[i]);
for(int i=q;i;--i)
if(!qu[i]) ans[++tot]=cur;
else inser(p[qu[i]]);
for(int i=tot;i;--i) printf("%.2lf\n",ans[i]);
}
然后有一个要注意的地方就是set的弱排序原理, 就是重载<号之后, 出现了
所以把重载<的地方改得那么鬼畜…
写完过了个样例交上去WA了… 然后就调呗.. 但是试了各种边界数据都没事…
甚至还把上面cf70D的代码贴过来但还是WA..
最后发现: 竟然是
for(int i=1;i<=m;++i) if(!del[i]) inser(p[i]);里面的
void inser(vec p);的函数传一个int变量不会报错啊 !!!!!
C++的强制转换太高端了吧…
不过这个题只维护一个上凸壳, 所以还有一种看上去好像更简单一点的做法是用平衡树来维护坐标, 这样坐标的排序就是以
但是最近觉得极角序非常地好用, 所以就写极角序咯…
用极角序这种判个<就要慢好多的 竟然还能跑进第一页也是挺神奇的…
有些时候凸包上的信息满足一些特殊条件的时候还可以用cdq分治来维护什么的…
据说代码会比平衡树好写, 但是平衡树如果写set就不一样了吧2333…..
不过这个应该很快就会去学…
那么学完了再说…