• 【笔记篇】最良心的计算几何学习笔记(四)


    本文的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]));
    }
    

    注意事项

    1. 按夹角排序的时候建议使用叉积.. 因为好写、不使用除法精度好等原因.

    2. 记得加上最后一条边(还记得多边形的首尾相接吧).

    3. (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]]));
    }
    

    注意事项

    1. 就是枚举的时候不要枚举自己.
    2. 共线的时候要选最远的.(看到里面内个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. 以上的算法都没有涉及到1个点 2个点的情况 如果可能出现的话要记得特判..
    2. 还有一种可能需要进行特判的就是构不成凸多边形(就是都在一条直线上的情况)...
    3. 实数比较的时候记得用cmp... 主要是判断相等的时候= =
    4. 凸包中最后一条边记得考虑.
    5. 有些时候要考虑正负号和绝对值的问题. 比如面积. 比如距离.

    代码已经上传到github

    好的, 完结撒花~

  • 相关阅读:
    apk逆向工具总结
    DAY...
    DAY1-Flask项目
    Django笔记
    DAY8-Python学习笔记
    DAY7-Python学习笔记
    DAY6-Python学习笔记
    DAY5-Python学习笔记
    DAY4-Python学习笔记
    DAY3-Python学习笔记
  • 原文地址:https://www.cnblogs.com/enzymii/p/8413431.html
Copyright © 2020-2023  润新知