本文的github地址在这里~
Emmmm 今天的主题是:
凸包(Convex Hull)
But what is "凸包"?
You can find it here.
WTF???
没事, 我知道wiki你萌看不下去.(因为我也看不下去←_←
但是我们不是有baidu嘛..
说的通俗但不严谨一些, 就是能包含平面上所有给定点的最小的凸多边形.
利用之前学过的知识, 我们已经能解决不少计算几何方面的问题, 甚至都能计算任意多边形的面积了. 但是如何高效地求凸包还是不那么显然的, 这就需要学习一些dalao们的算法了.
目前来看求凸包的常见方法有: 暴力((O(n^3))) Graham扫描法((O(nlogn))) Jarvis步进法((O(nH))) 快包算法((O(nlogn))) 等
但是
暴力
复杂度显然不怎么科学, 而且做法也特别显然.
所以我们就不说了.
Graham扫描法
wiki上讲得非常好.. 还有图.. 不过生肉啃起来有点累OvO
试图google翻译 但是机翻翻出来的真的不是人话...
但是还是要认真地学一学的(还可以顺便提高英语阅读水平不是→_→)
算法流程
我们一步一步地剖析这个算法.
首先先给出给定的点
step1
找到(y)坐标最小的点(P), 如果(y)相同则找(x)更小的.
这个不需要排序啦~一般读入的时候处理一下就好了, 这一步的复杂度是(O(n)).
step2
然后我们将所有点(P_i)按照(vec{PP_i})与(x)轴的夹角( heta)排序..
当然这个排序比较的时候并不需要真的算出夹角(三角函数是很昂贵的...)
我们知道(cos<vec a,vec b>=frac{vec acdotvec b}{|vec a||vec b|})
而且夹角的范围是([0,pi)), 这个区间内(cos)是单调递减的, 这样我们取(x)轴正方向的单位向量为(vec b)再算即可..
或者考虑斜率似乎也可以.. 或者用叉积..
排序算法只要是(O(nlogn))(或以内)的就可以了. (但是实数的基数排序不靠谱吧)
不过作为C++选手直接sort就okay了~
特别注意: 共线的时候要按照距离排序...
step3
建一个栈.
我们可以很容易的看出和证出(P_0), (P_1)和(P_{n-1})是凸包上的点.
于是将(P_0P_1)入栈.
然后按顺序枚举其他的点, 判断一下线段的拐向是顺时针还是逆时针.
如果是逆时针的话, 这个图形就还是凸的, 将其入栈.
而如果是顺时针的话, 这个图形就不再是凸的了, 我们需要退栈.
比如4-5这条边就在3-4的顺时针方向, 如果连上就变凹了, 所以4这个点是在凸包内的点, 退栈.
退到3后发现3-5就在2-3的逆时针了, 5进栈.
假如3-5仍然在2-3的顺时针方向, 则3也要退栈, 依此类推, 直到有逆时针方向或栈中只有(P)为止.(其实因为排过序了所以退到最后也会剩下(PP_1)这条边的...)
时间复杂度(O(n)), 虽然看上去好像最坏会一直退到最前面, 复杂度像是(O(n^2))的, 但实际上每条边最多被考虑两次(入栈的时候一次, 退栈的时候一次...)
按这个方法一路推到(P_{n-1})(前面说过肯定在凸包上), 就得到了凸包.
其实wiki上的内个动图就画的非常好, 我觉得大家都应该去看一下.
最后惊奇地发现复杂度主要取决于sort的复杂度...
总时间复杂度:(O(nlogn))
凸包的题目(包括但不限于模板题)是很多的, 这里用了luogu的一道例题.. [秘技:传送~]
实现代码
#include <cmath>
#include <cstdio>
#include <algorithm>
const int N=50101;
const double eps=1e-9;
int dcmp(const double &a){
if(fabs(a)<eps)return 0;
return a<0?-1:1;
}
inline double max(const double &a,const double &b){return dcmp(a-b)>0?a:b;}
struct point{
double x,y;
point(const double &X=0,const double &Y=0):x(X),y(Y){}
}p[N],stk[N];int tp,mi;
point operator -(const point &a,const point &b){
return point(a.x-b.x,a.y-b.y);
}
double operator ^(const point &a,const point &b){
return a.x*b.x+a.y*b.y;
}
double operator *(const point &a,const point &b){
return a.x*b.y-a.y*b.x;
}
double len(const point &a){
return sqrt(a^a);
}
//bool cmpa(const point &a,const point &b){
// point X=point(1,0),A=a-p[0],B=b-p[0];
// double coa=(A^X)/len(A),cob=(B^X)/len(B);
// 共线的时候选远的那个.保证前面的在凸包上.
// if(!dcmp(coa-cob)) return dcmp(len(A)-len(B))>0;
// return dcmp(coa-cob)>0;
//} //按夹角排序(点积版)
bool cmpa(const point &a,const point &b){
point A=a-p[0],B=b-p[0];
if(!dcmp(A*B)) return dcmp(len(A)-len(B))>0;
return dcmp(A*B)>0;
} //叉积版
void grahamScan(point* p,int n){
std::sort(p+1,p+n,cmpa);
stk[++tp]=p[0]; stk[++tp]=p[1];
for(int i=2;i<n;++i){
while(dcmp((p[i]-stk[tp])*(stk[tp]-stk[tp-1]))>=0&&tp>2) --tp;
//由于排过序共线的时候一定不优, 所以不逆时针而且能退栈(栈里大于两个点)就退栈
stk[++tp]=p[i]; //进栈
}
}
int main(){
int n; scanf("%d",&n); double my=1e9,ans=0;
for(int i=0;i<n;++i){
scanf("%lf%lf",&p[i].x,&p[i].y);
if(dcmp(p[i].y-my)<0||
(!dcmp(p[i].y-my)&&dcmp(p[i].x-p[mi].x)<0))
my=p[i].y,mi=i;
} std::swap(p[0],p[mi]);
grahamScan(p,n);
for(int i=1;i<tp;++i)
ans+=len(stk[i]-stk[i+1]);
printf("%.2lf",ans+len(stk[tp]-p[0]));
}
注意事项
-
按夹角排序的时候建议使用叉积.. 因为好写、不使用除法精度好等原因.
-
记得加上最后一条边(还记得多边形的首尾相接吧).
-
(P)点就不要参与排序了吧.
其实还是比较好写的..(毕竟只是个板子)
Jarvis步进法
其实wiki上面这个条目叫Gift wrapping algorithm(礼品包装算法)来着..
这个算法是输出敏感的, 复杂度(O(nH)), (H)是凸包上点的个数.
也就是说极限情况下(所有点都在凸包上)是(O(n^2))的, 但很多情况下(H)是小于(logn)的, 那么反而会比Graham扫描要快.(比如人口比较集中分布在城市, 郊区的人口较少之类的). 但是算法竞赛中还是要慎用...
算法流程
给定的点还是上面那些, 这里就不画了.
step1
找到一个一定在凸包上的点(P), 然后令(H=P).
这里我们找最左边的点(这次(y)应该是无所谓的)
方法同上, 时间复杂度(O(n)).
step2
枚举每一个点(H'), 找出所有的(vec{HH'})中最逆时针方向的一个(利用叉积即可)
可以肯定的是, 这个(H')也在凸包上.
令(H=H')
时间复杂度(O(n))
step3
重复step2, 直到(H)重新等于(P)为止.
这样我们就找到了所有在凸包上的点, 也就是说找到了这个凸包.
每次我们会找到一个凸包上的点, 然后进行一次step2, 所以step2总共会进行(H)次, 所以时间复杂度(O(H)*O(n)=O(nH))
还是比较简单的. wiki上说"在二维中, 礼品包装算法类似于在一组点上缠绕线(或包装纸)的过程".
而且这种算法似乎是可以扩展到更高维度的. 等以后再学吧..
实现代码
还是那道题吧.. 本来以为可能会被卡于是开了个O2但似乎并不会卡...
// luogu-judger-enable-o2
#include <cmath>
#include <cstdio>
const int N=10101;
const double eps=1e-9;
int dcmp(const double &a){
if(fabs(a)<eps)return 0;
return a<0?-1:1;
}
struct point{
double x,y;
point(const double &X=0,const double &Y=0):x(X),y(Y){}
}p[N];
point operator -(const point &a,const point &b){
return point(a.x-b.x,a.y-b.y);
}
double operator *(const point &a,const point &b){
return a.x*b.y-a.y*b.x;
}
double len(const point &a){
return sqrt(a.x*a.x+a.y*a.y);
}
int po[N],tot,mi;
void jarvisMarch(point *p,int n){
int h=mi;
do{
int h2=-1;
for(int i=0;i<n;++i)
if(h!=i&&(h2<0||dcmp((p[i]-p[h])*(p[h2]-p[h]))>0
||(!dcmp((p[i]-p[h])*(p[h2]-p[h]))
&&dcmp(len(p[h2]-p[h])-len(p[i]-p[h]))<0)))
h2=i;
po[++tot]=h=h2;
}while(h!=mi);
}
int main(){
int n;scanf("%d",&n); double ans=0,mx=1e9;
for(int i=0;i<n;++i){
scanf("%lf%lf",&p[i].x,&p[i].y);
if(p[i].x<mx) mx=p[i].x,mi=i;
}
jarvisMarch(p,n);
for(int i=1;i<tot;++i)
ans+=len(p[po[i]]-p[po[i+1]]);
printf("%.2lf",ans+len(p[po[tot]]-p[po[1]]));
}
注意事项
- 就是枚举的时候不要枚举自己.
- 共线的时候要选最远的.(看到里面内个if语句写了多么一坨了么←_←
快包算法
wiki传送门
跟快排十分相似, 平均复杂度(O(nlogn)), 最好的情况能达到(O(n)), 最坏的情况会变成(O(n^2)). 但还是一种比较常见的做法.
思路也和快排比较相似, 用的是分治.
算法流程
step1
找到最左边和最右边的点(A,B), 它们一定在凸包上.
然后把他们连接起来.
step2
这条线上面找到这条线最远的点(C), 这个点一定在凸包上.
然后递归处理(AC,BC)之间的凸壳.
step3
找从(B)到(A)的上凸壳, 也就是下凸壳..
这个复杂度的证明方式我是真的不会了..
快包可以比较方便地只找上/下凸壳.
实现代码
这个方法看上去很简单但是写起来很鬼畜啊= =为了方便理解和写, 我们放一下伪代码.
quickHull(点集S,有向线段P[A->B]){
选取距离P最远的点C
有向线段M[A->C] N[C->B]
点集SL{X|X∈S且在M左侧}
点集SR{X|X∈S且在N左侧}
quickHull(SL,M)
C是凸包上的点
//这里务必采用中序的方式,保证凸包上的点是按顺序输出的
quickHull(SR,N)
}
根据伪代码我们可以整合出代码来
#include <cmath>
#include <cstdio>
#include <vector>
using namespace std;
const int N=10101;
const double eps=1e-9;
int dcmp(const double &a){
if(fabs(a)<eps)return 0;
return a<0?-1:1;
}
struct point{
double x,y;
point(const double &X=0,const double &Y=0):x(X),y(Y){}
}hull[N];
typedef vector<point> pset;
typedef pset::iterator pit;
int mi,tot;
inline bool cmp(const point &a,const point &b){
if(!dcmp(a.x-b.x)) return dcmp(a.y-b.y)<0;
return dcmp(a.x-b.x)<0;
}
point operator -(const point &a,const point &b){
return point(a.x-b.x,a.y-b.y);
}
double operator *(const point &a,const point &b){
return a.x*b.y-a.y*b.x;
}
double operator ^(const point &a,const point &b){
return a.x*b.x+a.y*b.y;
}
double len(const point &a){
return sqrt(a.x*a.x+a.y*a.y);
}
double ptDisSeg(const point &p,const point &q0,const point &q1){
if(!dcmp((p-q0)^(q1-q0))) return len(p-q0);
if(!dcmp((p-q1)^(q0-q1))) return len(p-q1);
return fabs((p-q0)*(q1-q0)/len(q1-q0));
}
void quickHull(pset s,const point &a,const point &b){
pset sl,sr; double dis=-1e9; point c;
for(pit i=s.begin();i!=s.end();++i){
double d=ptDisSeg(*i,a,b);
if(d>dis) dis=d,c=*i;
}
for(pit i=s.begin();i!=s.end();++i){
if(dcmp((*i-a)*(c-a))<0)
sl.push_back(*i);
if(dcmp((*i-c)*(b-c))<0)
sr.push_back(*i);
}
if(sl.size())quickHull(sl,a,c);
hull[++tot]=c;
if(sr.size())quickHull(sr,c,b);
}
int main(){
int n; scanf("%d",&n); pset p;
double my=1e9,ans=0;
for(int i=0;i<n;++i){
double x,y; scanf("%lf%lf",&x,&y);
if(x<my) mi=i,my=x; p.push_back(point(x,y));
} swap(p[0],p[mi]); hull[++tot]=p[0]; p.erase(p.begin());
quickHull(p,hull[1],hull[1]);
for(int i=1;i<tot;++i)
ans+=len(hull[i]-hull[i+1]);
printf("%.2lf",ans+len(hull[tot]-hull[1]));
}
总结
其他的还有不少高端大气上档次的算法, 比如似乎有(O(nlogH))的Chan's Algorithm什么的. 有时间再研究吧= = 这几种方法应该够用了.
来对注意事项做一波汇总
- 以上的算法都没有涉及到1个点 2个点的情况 如果可能出现的话要记得特判..
- 还有一种可能需要进行特判的就是构不成凸多边形(就是都在一条直线上的情况)...
- 实数比较的时候记得用cmp... 主要是判断相等的时候= =
- 凸包中最后一条边记得考虑.
- 有些时候要考虑正负号和绝对值的问题. 比如面积. 比如距离.
代码已经上传到github
好的, 完结撒花~