zoukankan      html  css  js  c++  java
  • [学习笔记]计算几何

    堪称明明看得见,就是写不出的一类恶心题。

    通常细节颇多。

    一旦方法选择合适,代码量和效率都会提升不少。

    推荐:

    「计算几何」计算几何从入门到入土

     计算几何入门

     

    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);
        }
    };

    (加减为了向量方便赋值)

    两(多)点距离公式

    double dis(po a,po b){
        return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
    }
    距离公式

    中点

    (x=(x1+x2)/2,y=(y1+y2)/2)

    向量

    struct vec{
        double x,y;
        vec(){}
        vec(po a){
            x=a.x,y=a.y;
        }
        double len(){
            return sqrt(x*x+y*y);
        }
    };
    向量

    向量加减和数乘略了。

    len是向量的模长。

    点积

    double dot(vec a,vec b){
        return a.x*b.x+a.y*b.y;
    }
    点积

    点积用处:可以通过两个向量的点积正负,判断向量夹角是锐角、直角、钝角

    叉积

    double cross(vec a,vec b){
        return a.x*b.y-a.y*b.x;
    }
    叉积

    意义,两个向量构成的平行四边形的有向面积

    AxB ,如果有A在B的逆时针方向(不超180度),那么面积是负的。否则是正的。

    取绝对值就是面积。

    向量旋转

    https://zhidao.baidu.com/question/562549647.html

    常用方法

    精度控制

    设置eps 1e-7~1e-10

    int Fabs(double t){
        if(fabs(t)<eps) return 0;
        if(t>0) return 1;
        return -1;
    }
    Fabs

    点到直线距离(给出三个点)

    叉积算平行四边形面积,然后除以底的长度,得到高(距离)

    double hei(po p,po a,po b){
        vec t1=vec(a-b),t2=vec(p-b);
        return fabs(cross(t1,t2))/dis(a,b);
    }
    点到直线距离

    点到线段距离

    点积判断和端点形成的角是否是钝角。是钝角,就只能是和线段两个端点的距离取min

    线段是否相交

     跨立实验

    如果两个线段端点互相都在对方的两侧,那么相交。

    点在多边形内部

    上面第二篇博客有写到。

    就是找一条射线,判断和多边形交点奇偶,

    要枚举每条多边形上的边。(O(n))

    为了排除各种和点相交,和边重合的情况,

    为了统一起见,我们在计算射线L和多边形的交点的时候,1。对于多边形的水平边不作考虑;2。对于多边形的顶点和L相交的情况,如果该顶点是其所属的边上纵坐标较大的顶点,则计数,否则忽略;3。对于P在多边形边上的情形,直接可判断P属于多边行。 其中做射线L的方法是:设P'的纵坐标和P相同,横坐标为正无穷大(很大的一个正数),则P和P'就确定了射线L。

    ——https://blog.csdn.net/qq_35776579/article/details/54836612

    求任意多边形面积

     三角剖分

    任意选择一个点,然后逆时针/顺时针枚举所有边,和这个点构成三角形,叉积计算有向面积,然后做和即可。

    记得取绝对值,然后除以二

    多边形重心

    同样方法将多边形三角剖分

    算出每个三角形的重心,套用质点组的重心公式即可

    公式见上面第一篇博客。

    三角形重心:

    $((x1+x2+x3)/3,(y1+y2+y3)/3)$

    两条直线的交点

    https://blog.csdn.net/dgq8211/article/details/7952825

    可以根据三角形相似。

    得到DP和PC的比值,就可以通过D、C的坐标,算出来P的坐标。

    三角形的面积必须是有向面积,因为情况不止图的这一种。

     这个时候,必须求ABD,ABC的有向面积。这样正负号才是对的。

    式子:

    $P=frac{AB×AC*D+AD×AB*C}{AB×AC+AD×AB}$

    $AB$等都是向量,从A指向B

    $×$是叉乘,$*$是和点乘(横纵坐标分别相乘)

    凸包

    长这个样子的一个包。

     

    Graham 算法

     ①选择左下点,按极角排序,

    ②用一个栈维护凸包上的点,相邻的凹凸情况,用叉积判断。

    由于排序,所有O(nlogn)

     代码见下面模板处

    旋转卡壳

     

     用于求凸包最大直径,或者平面最远点对。

    做法:枚举一条边,然后找到和这条边平行的边在另一端的切点。(距离这个边最远的点)

    具体移动点的话,可以通过叉积判断面积,找最大的位置。(当然,你也可以三分。。。)

    然后这个位置和这个边的距离较大的一个,就是这次可以更新ans的值。

    (证明:

    只需证明,最远点对之一A,一定是和"另一个点B和与B相邻的某个点C"形成的B-C边距离最远的点。

    应该可以证明因为找不出反例233333

    发现,这个切点的位置是具有单调性的,即可O(n)解决这个问题。

    (其实,Graham的排序还是O(nlogn)的)

    半平面交

    另一篇博客

    [学习笔记]半平面交

    其他

    [学习笔记]最小圆覆盖

    [学习笔记]自适应辛普森积分

    [WC2013]平面图——平面图点定位

     [学习笔记]闵可夫斯基和

    例题

     【模板】二维凸包 / [USACO5.1]圈奶牛Fencing the Cows

    #include<bits/stdc++.h>
    #define reg register int
    #define il inline
    #define numb (ch^'0')
    using namespace std;
    typedef long long ll;
    il void rd(int &x){
        char ch;x=0;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=10000+5;
    const double eps=1e-8;
    int n;
    struct po{
        double x,y;
        po(){};
        po(double xx,double yy){
            x=xx;y=yy;
        }
        po operator +(po &b){
            return po(x+b.x,y+b.y);
        }
        po operator -(po &b){
            return po(x-b.x,y-b.y);
        }
        double dis(const po &b){
            return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y));
        }
    }p[N];
    struct vec{
        double x,y;
        void init(const po &b){
            x=b.x,y=b.y;
        }
        double cross(const vec &b){
            return x*b.y-y*b.x;
        }
        double len(){
            return sqrt(x*x+y*y);
        }
    };
    int Fabs(double t){
        if(fabs(t)<eps) return 0;
        if(t>0) return 1;
        return -1;
    }
    bool cmp(po a,po b){
        vec a1,b1;
        a1.init(a-p[1]);b1.init(b-p[1]);
        double t=a1.cross(b1);
        if(Fabs(t)) return t>0;
        return a1.len()<b1.len();
    }
    double ans;
    vector<po>mem;
    int main(){
        scanf("%d",&n);
        if(n==0||n==1){
            printf("0.00");return 0;
        }
        for(reg i=1;i<=n;++i){
            scanf("%lf%lf",&p[i].x,&p[i].y);
        }int id=1;
        for(reg i=2;i<=n;++i){
            if(p[i].x<p[id].x||((p[i].x==p[id].x)&&(p[i].y<p[id].y))) id=i;
        }
        if(id!=1) swap(p[1],p[id]);
        sort(p+2,p+n+1,cmp);
        mem.push_back(p[1]);
        for(reg i=2;i<=n;++i){
            while(mem.size()>2){
                vec t1,t2;
                t1.init(p[i]-mem.back());t2.init(mem.back()-mem[mem.size()-2]);
                double t=t2.cross(t1);
                if(Fabs(t)==-1) mem.pop_back();
                else break;
            }
            mem.push_back(p[i]);
        }
        for(reg i=0;i<mem.size();++i){
            ans+=mem[i].dis(mem[(i+1)%mem.size()]);
        }
        if(n==2){
            ans/=2.0;
        }
        printf("%.2lf",ans);
        return 0;
    }
    
    }
    int main(){
        Miracle::main();
        return 0;
    }
    
    /*
       Author: *Miracle*
       Date: 2018/11/24 15:03:36
    */
    二维凸包

    [SCOI2007]最大土地面积

    直接处理四边形不好处理。

    考虑,拼成两个三角形。

    我们枚举一条四边形对角线A-B(固定A),那么三角形就是在对角线两边各选择一个点C、D,使得和A-B形成的三角形面积最大。

    发现,这个C、D位置,随着B的移动,也是单调移动的。

    本质是维护一个三指针,然后O(n^2)解决。

    // luogu-judger-enable-o2
    #include<bits/stdc++.h>
    #define reg register int
    #define il inline
    #define numb (ch^'0')
    using namespace std;
    typedef long long ll;
    il void rd(int &x){
        char ch;x=0;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=2002;
    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);
        }
    }a[N];
    struct vec{
        double x,y;
        vec(){}
        vec(po a){
            x=a.x,y=a.y;
        }
        double len(){
            return sqrt(x*x+y*y);
        }
    };
    
    double dis(po a,po b){
        return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
    }
    double cross(vec a,vec b){
        return a.x*b.y-a.y*b.x;
    }
    int Fabs(double t){
        if(fabs(t)<eps) return 0;
        if(t>0) return 1;
        return -1;
    }
    double hei(po p,po a,po b){
        vec t1=vec(a-b),t2=vec(p-b);
        return fabs(cross(t1,t2))/dis(a,b);
    }
    bool cmp(po x,po y){
        double t=cross(vec(x-a[1]),vec(y-a[1]));
        if(Fabs(t)) return t>0;
        return vec(x-a[1]).len()<vec(y-a[1]).len();
    }
    po sta[N];
    int top;
    double ans;
    int main(){
        scanf("%d",&n);
        for(reg i=1;i<=n;++i){
            scanf("%lf%lf",&a[i].x,&a[i].y);
        }int id=1;
        for(reg i=2;i<=n;++i){
            if(a[i].x<a[id].x||((a[i].x==a[id].x)&&(a[i].y<a[id].y))) id=i;
        }
        if(id!=1) swap(a[1],a[id]);
        sort(a+2,a+n+1,cmp);
        sta[++top]=a[1];
        for(reg i=2;i<=n;++i){
            while(top>2&&cross(vec(sta[top]-sta[top-1]),vec(a[i]-sta[top]))<=0) --top;
            sta[++top]=a[i];
        }
    //    cout<<" top "<<top<<endl;
    //    for(reg i=1;i<=top;++i){
    //        cout<<i<<" : "<<sta[i].x<<" "<<sta[i].y<<endl;
    //    }
        for(reg i=1;i<=top;++i){
            int p1=i,p2=i+1;
            for(reg j=i+1;j<=top;++j){
                while(p1+1<=j&&(Fabs(hei(sta[p1+1],sta[i],sta[j])-hei(sta[p1],sta[i],sta[j]))>0)) ++p1;
                while(p2%top+1!=i+1&&(Fabs(hei(sta[p2%top+1],sta[i],sta[j])-hei(sta[p2],sta[i],sta[j])>0))) p2=p2%top+1;
                if((hei(sta[p1],sta[i],sta[j])+hei(sta[p2],sta[i],sta[j]))*dis(sta[i],sta[j])/2>ans){
                    //cout<<" new "<<i<<" "<<p1<<" "<<j<<" "<<p2<<endl;
                    ans=max(ans,(hei(sta[p1],sta[i],sta[j])+hei(sta[p2],sta[i],sta[j]))*dis(sta[i],sta[j])/2);
                    //cout<<" ans is "<<ans<<endl;
                }
            }
        }
        printf("%.3lf",ans);
        return 0;
    }
    
    }
    int main(){
        Miracle::main();
        return 0;
    }
    
    /*
       Author: *Miracle*
       Date: 2018/11/24 16:06:58
    */
    最大土地面积

    Beauty Contest

    很裸的旋转卡壳求平面最远点对。

    #include<bits/stdc++.h>
    #define reg register int
    #define il inline
    #define numb (ch^'0')
    using namespace std;
    typedef long long ll;
    il void rd(int &x){
        char ch;x=0;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=50000+5;
    struct po{
        int x,y;
        po(){}
        po(int xx,int 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);
        }
    }a[N];
    struct vec{
        int x,y;
        vec(){}
        vec(po a){
            x=a.x,y=a.y;
        }
        int len(){
            return x*x+y*y;
        }
    };
    int dis(po a,po b){
        return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
    }
    int cross(vec a,vec b){
        return a.x*b.y-a.y*b.x;
    }
    bool cmp(po x,po y){
        vec t1=vec(x-a[1]),t2=vec(y-a[1]);
        int t=cross(t1,t2);
        if(t) return t>0;
        return t1.len()<t2.len();
    }
    po sta[N];
    int top;
    int ans;
    int n;
    int main(){
        scanf("%d",&n);
        for(reg i=1;i<=n;++i){
            rd(a[i].x);rd(a[i].y);
        }int id=1;
        for(reg i=2;i<=n;++i){
            if(a[i].x<a[id].x||((a[i].x==a[id].x)&&a[i].y<a[id].y)) id=i;
        }
        if(id!=1)swap(a[1],a[id]);
        sort(a+2,a+n+1,cmp);
        
        sta[++top]=a[1];
        for(reg i=2;i<=n;++i){
            while(top>2&&cross(vec(sta[top]-sta[top-1]),vec(a[i]-sta[top]))<=0) --top;
            sta[++top]=a[i];
        }
        int T=1;int B;
        for(reg A=1;A<=top;++A){
            B=A-1;if(B==0) B=top;
            while(abs(cross(vec(sta[T%top+1]-sta[A]),vec(sta[B]-sta[A])))>abs(cross(vec(sta[T]-sta[A]),vec(sta[B]-sta[A]))))T=T%top+1;
            ans=max(ans,max(dis(sta[T],sta[A]),dis(sta[T],sta[B])));
        }
        printf("%d",ans);
        return 0;
    }
    
    }
    int main(){
        Miracle::main();
        return 0;
    }
    
    /*
       Author: *Miracle*
       Date: 2018/11/24 21:17:16
    */
    旋转卡壳

    [HNOI2007]最小矩形覆盖

    好题。见另一篇博客

     

  • 相关阅读:
    【java学习笔记】反射基础
    【java学习笔记】线程
    【java学习笔记】Properties
    【java学习笔记】序列化、反序列化
    Java Review (二十六、集合----- Set 集合)
    Java Review (二十五、集合----- Iterator接口)
    Java Review (二十四、集合-----Collection 接口)
    Java Review (二十三、集合-----概述)
    Java Review (二十二、正则表达式)
    Java Review (二十一、基础类库----日期、时间类)
  • 原文地址:https://www.cnblogs.com/Miracevin/p/10013826.html
Copyright © 2011-2022 走看看