[计算几何]【学习笔记】

深碍√TFBOYSˉ_ 2022-03-30 08:35 262阅读 0赞

学计算几何专题已经两年了(还不是因为刚刚过年了.....)

**哼 都写完了网页崩溃全没了 生气了 哼 本来已经12点了超想睡觉的现在又要重写一遍 哼 我不睡觉了**

先扔一个完整版方便复制

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

//By Candy?
    #include <iostream>
    #include <cstdio>
    #include <cstring>
    #include <algorithm>
    #include <cmath>
    using namespace std;
    typedef long long ll;
    const int N=1e4,INF=1e9;
    inline int read(){
        char c=getchar();int x=0,f=1;
        while(c<'0'||c>'9'){
         if(c=='-')f=-1; c=getchar();}
        while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}
        return x*f;
    }
    
    const double eps=1e-8;
    const double pi=acos(-1);
    
    inline int sgn(double x){
        if(abs(x)<eps) return 0;
        else return x<0?-1:1;
    }
    
    struct Vector{
        double x,y;
        Vector(double a=0,double b=0):x(a),y(b){}
        bool operator <(const Vector &a)const{
            return sgn(x-a.x)<0||(sgn(x-a.x)==0&&sgn(y-a.y)<0);
        }
        void print(){printf("%lf %lf\n",x,y);}
    };
    typedef Vector Point;
    Vector operator +(Vector a,Vector b){
         return Vector(a.x+b.x,a.y+b.y);}
    Vector operator -(Vector a,Vector b){
         return Vector(a.x-b.x,a.y-b.y);}
    Vector operator *(Vector a,double b){
         return Vector(a.x*b,a.y*b);}
    Vector operator /(Vector a,double b){
         return Vector(a.x/b,a.y/b);}
    bool operator ==(Vector a,Vector b){
         return sgn(a.x-b.x)==0&&sgn(a.y-b.y)==0;}
    
    double Dot(Vector a,Vector b){
         return a.x*b.x+a.y*b.y;}
    double Cross(Vector a,Vector b){
         return a.x*b.y-a.y*b.x;}
    
    double Len(Vector a){
         return sqrt(Dot(a,a));}
    double Len2(Vector a){
         return Dot(a,a);}
    double Angle(Vector a,Vector b){
        return acos(Dot(a,b)/Len(a)/Len(b));
    }
    Vector Normal(Vector a){
        return Vector(-a.y,a.x);//counterClockwise
    }
    Vector Rotate(Vector a,double rad){
        return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
    }
    
    //Line:use p and v
    //struct Line{
    //    Point p;
    //    Vector v;
    //    Line(){}
    //    Line(Point p,Vector v):p(p),v(v){}
    //    bool operator <(const Line a)const{
    //        return sgn(Cross(v,a.v))>=0;
    //    }
    //};
    //bool OnLeft(Line l,Point p){
    //    return sgn(Cross(l.v,p-l.p))>=0;
    //}
    //Point LI(Line a,Line b){
    //    Vector v=a.p-b.p;
    //    double t=Cross(b.v,v)/Cross(a.v,b.v);
    //    return a.p+a.v*t;
    //}
    
    //Line:use s and t
    struct Line{
        Point s,t;
        Line(){}
        Line(Point a,Point b):s(a),t(b){}
        bool operator <(Line a)const{
            return sgn(Cross(t-s,a.t-a.s))>=0;
        }
    };
    
    double DisTL(Point p,Point a,Point b){
        Vector v1=b-a,v2=p-a;
        return abs(Cross(v1,v2)/Len(v1));
    }
    double DisTS(Point p,Point a,Point b){
        if(a==b) return Len(p-a);
        Vector v1=b-a,v2=p-a,v3=p-b;
        if(sgn(Dot(v1,v2))<0) return Len(v2);
        else if(sgn(Dot(v1,v3))>0) return Len(v3);
        else return abs(Cross(v1,v2)/Len(v1));
    }
    bool OnLeft(Line l,Point p){
        return sgn(Cross(l.t-l.s,p-l.s))>=0;
    }
    bool OnSeg(Point p,Point a,Point b){
        return DisTL(p,a,b)==0&&sgn(Dot(p-a,p-b)<0&&!(p==a))&&!(p==b);
    }
    
    Point LI(Line a,Line b){
        Vector v=a.s-b.s,v1=a.t-a.s,v2=b.t-b.s;
        double t=Cross(v2,v)/Cross(v1,v2);
        return a.s+v1*t;
    }
    bool isLSI(Line l1,Line l2){
        Vector v=l1.t-l1.s,u=l2.s-l1.s,w=l2.t-l1.s;
        return sgn(Cross(v,u))!=sgn(Cross(v,w));
    }
    bool isSSI(Line l1,Line l2){
        return isLSI(l1,l2)&&isLSI(l1,l2);
    }
    //---chong he
    //bool isSSI(Line l1,Line l2){
    //    Vector v1=l1.t-l1.s,v2=l2.t-l2.s;
    //    if(sgn(Cross(v1,v2))==0){
    //        int flag=0;
    //        Vector u=l2.s-l1.s,w=l2.t-l1.s;
    //        if(sgn(Dot(u,w))<0) flag=1;
    //        u=l2.s-l1.t,w=l2.t-l1.t;
    //        if(sgn(Dot(u,w))<0) flag=1;
    //        return flag;
    //    }
    //    else return isLSI(l1,l2)&&isLSI(l2,l1);
    //}
    
    Point Circumcenter(Point a,Point b,Point c){
        Point p=(a+b)/2,q=(a+c)/2;
        Vector v=Normal(b-a),u=Normal(c-a);
        if(sgn(Cross(v,u))==0){
            if(sgn(Len(a-b)+Len(b-c)-Len(a-c))==0) return (a+c)/2;
            if(sgn(Len(a-b)+Len(a-c)-Len(b-c))==0) return (b+c)/2;
            if(sgn(Len(a-c)+Len(b-c)-Len(a-b))==0) return (a+b)/2;
        }
        return LI(Line(p,p+v),Line(q,q+u));
    }
    
    Point Barycenter(Point a,Point b,Point c){
        return (a+b+c)/3;
    }
    
    bool cmpPolar(Point a,Point b){
        return sgn(Cross(a,b))>0;
    }
    
    int PointInPolygon(Point p,Point poly[],int n){
        int wn=0;
        for(int i=1;i<=n;i++){
            if(sgn(DisTS(p,poly[i],poly[i%n+1]))==0) return -1;
            int k=sgn(Cross(poly[i%n+1]-poly[i],p-poly[i])),
            d1=sgn(poly[i].y-p.y),d2=sgn(poly[i%n+1].y-p.y);
            if(k>0&&d1<=0&&d2>0) wn++;
            if(k<0&&d2<=0&&d1>0) wn--;
        }
        return (bool)wn;
    }
    
    double PolygonArea(Point p[],int n){
        double s=0;
        for(int i=2;i<n;i++) s+=Cross(p[i]-p[1],p[i+1]-p[1]);
        return abs(s/2);
    }
    
    bool isConvex(Point poly[],int n){
        int last=0,now=0;
        for(int i=1;i<=n;i++){
            now=sgn(Cross(poly[i%n+1]-poly[i],poly[(i+1)%n+1]-poly[i%n+1]));
            if(last==0||now==0||now*last>0) last=now;
            else return false;
        }
        return true;
    }
    
    int ConvexHull(Point p[],int n,Point ch[]){
         //cannot handle repeat point
        sort(p+1,p+1+n);
        int m=0;
        for(int i=1;i<=n;i++){
            while(m>1&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--;
            ch[++m]=p[i];
        }
        int k=m;
        for(int i=n-1;i>=1;i--){
         //n-1
            while(m>k&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--;
            ch[++m]=p[i];
        }
        if(n>1) m--;//the first  rpoint
        return m;
    }
    
    double RotatingCalipers(Point p[],int n){
        if(n==1) return 0;
        if(n==2) return Len(p[1]-p[2]);
        int now=1;
        double ans=0;
        p[n+1]=p[1];
        for(int i=1;i<=n;i++){
            while(sgn(DisTL(p[now],p[i],p[i+1])-DisTL(p[now+1],p[i],p[i+1]))<=0) now=now%n+1;
            ans=max(ans,Len(p[now]-p[i]));
            ans=max(ans,Len(p[now]-p[i+1]));
        }
        return ans;
    }
    
    void iniPolygon(Point p[],int &n,double inf){
        n=0;
        p[++n]=Point(-inf,-inf);
        p[++n]=Point(inf,-inf);
        p[++n]=Point(inf,inf);
        p[++n]=Point(-inf,inf);
    }
    Point t[N];int tn;
    void CutPolygon(Point p[],int &n,Point a,Point b){
         //get the left of a->b
        tn=0;
        Point c,d,e;
        for(int i=1;i<=n;i++){
            c=p[i],d=p[i%n+1];
            if(sgn(Cross(b-a,c-a))>=0) t[++tn]=c;
            if(isLSI(Line(a,b),Line(c,d))){
                e=LI(Line(a,b),Line(c,d));
                t[++tn]=e;
            }
        }
        n=tn;for(int i=1;i<=n;i++)p[i]=t[i];
    }
    //    iniPolygon(q,m,INF);
    //    for(int i=1;i<=n;i++) CutPolygon(q,m,p[i%n+1],p[i]);
    
    double minCircleCover(Point p[],int n,Point &c){
        random_shuffle(p+1,p+1+n);
        c=p[1];
        double r=0;
        for(int i=2;i<=n;i++)
            if(sgn(Len(c-p[i])-r)>0){
                c=p[i],r=0;
                for(int j=1;j<i;j++)
                    if(sgn(Len(c-p[j])-r)>0){
                        c=(p[i]+p[j])/2,r=Len(c-p[i]);
                        for(int k=1;k<j;k++)
                            if(sgn(Len(c-p[k])-r)>0){
                                c=Circumcenter(p[i],p[j],p[k]);
                                r=Len(c-p[i]);
                            }
                    }
            }
        return r;
    }
    
    
    //jie xi ji he
    double a,b,c;
    double f(Point p){
         return a*p.x+b*p.y+c;}
    Point abcLI(Line l){
        double u=abs(f(l.s)),v=abs(f(l.t));
        return Point(l.s.x*v+l.t.x*u,l.s.y*v+l.t.y*u)/(u+v);
    }
    
    double F(double x){
         return x;} //function
    inline double cal(double l,double r,double fl,double fr,double fm){
        return (fl+fr+4*fm)*(r-l)/6;
    }
    double Simpson(double l,double r,double now,double fl,double fr,double fm){
        double mid=(l+r)/2,flm=F((l+mid)/2),frm=F((mid+r)/2);
        double p=cal(l,mid,fl,fm,flm),q=cal(mid,r,fm,fr,frm);
        if(sgn(now-p-q)==0) return now;
        else return Simpson(l,mid,p,fl,fm,flm)+Simpson(mid,r,q,fm,fr,frm);
    }
    
    //sao miao xian

计算几何模板

--------------------

并不想多说什么 ljw的课件已经很好了

**1.精度问题**

const double eps=1e-8;
    inline int sgn(double x){
        if(abs(x)<eps) return 0;
        else return x<0?-1:1;
    }

**2.向量 点**

![ContractedBlock.gif][] 向量+-\*／ 比较

点积 cosΘ

叉积 sinΘ aXb>0 b在a的左面(逆时针方向)

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double Dot(Vector a,Vector b){
         return a.x*b.x+a.y*b.y;}
    double Cross(Vector a,Vector b){
         return a.x*b.y-a.y*b.x;}
    
    double Len(Vector a){
         return sqrt(Dot(a,a));}
    double Len2(Vector a){
         return Dot(a,a);}
    double Angle(Vector a,Vector b){
        return acos(Dot(a,b)/Len(a)/Len(b));
    }

点积叉积 长度 夹角

法向量 (-y,x) 逆时针90

旋转 想一个与坐标轴平行的直角然后一推就可以了

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

Vector Normal(Vector a){
        return Vector(-a.y,a.x);//counterClockwise
    }
    Vector Rotate(Vector a,double rad){
        return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
    }

法向量 旋转

**3.直线 线段**

感觉用两点比较方便

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

struct Line{
        Point s,t;
        Line(){}
        Line(Point a,Point b):s(a),t(b){}
        bool operator <(Line a)const{
            return sgn(Cross(t-s,a.t-a.s))>=0;
        }
    };

线

点到直线距离：面积法

线段最短距离：用点积特判是点到端点到距离的情况

点线位置关系 点在线段上

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double DisTL(Point p,Point a,Point b){
        Vector v1=b-a,v2=p-a;
        return abs(Cross(v1,v2)/Len(v1));
    }
    double DisTS(Point p,Point a,Point b){
        if(a==b) return Len(p-a);
        Vector v1=b-a,v2=p-a,v3=p-b;
        if(sgn(Dot(v1,v2))<0) return Len(v2);
        else if(sgn(Dot(v1,v3))>0) return Len(v3);
        else return abs(Cross(v1,v2)/Len(v1));
    }
    bool OnLeft(Line l,Point p){
        return sgn(Cross(l.t-l.s,p-l.s))>=0;
    }
    bool OnSeg(Point p,Point a,Point b){
        return DisTL(p,a,b)==0&&sgn(Dot(p-a,p-b)<0&&!(p==a))&&!(p==b);
    }

点线

直线相交交点：面积求比例，正负画个图看一下是否符合就行了

直线与线段相交：端点叉积不同

线段与线段相交：两遍直线与线段(有时候特判重合)

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

Point LI(Line a,Line b){
        Vector v=a.s-b.s,v1=a.t-a.s,v2=b.t-b.s;
        double t=Cross(v2,v)/Cross(v1,v2);
        return a.s+v1*t;
    }
    bool isLSI(Line l1,Line l2){
        Vector v=l1.t-l1.s,u=l2.s-l1.s,w=l2.t-l1.s;
        return sgn(Cross(v,u))!=sgn(Cross(v,w));
    }
    bool isSSI(Line l1,Line l2){
        return isLSI(l1,l2)&&isLSI(l1,l2);
    }

相交 交点

**4.三角形**

外心：两个中垂线交点，特判三点共线

重心：三个点平均

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

Point Circumcenter(Point a,Point b,Point c){
        Point p=(a+b)/2,q=(a+c)/2;
        Vector v=Normal(b-a),u=Normal(c-a);
        if(sgn(Cross(v,u))==0){
            if(sgn(Len(a-b)+Len(b-c)-Len(a-c))==0) return (a+c)/2;
            if(sgn(Len(a-b)+Len(a-c)-Len(b-c))==0) return (b+c)/2;
            if(sgn(Len(a-c)+Len(b-c)-Len(a-b))==0) return (a+b)/2;
        }
        return LI(Line(p,p+v),Line(q,q+u));
    }
    
    Point Barycenter(Point a,Point b,Point c){
        return (a+b+c)/3;
    }

心

## 5.多边形 ##

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

bool cmpPolar(Point a,Point b){
        return sgn(Cross(a,b))>0;
    }

极角排序

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double PolygonArea(Point p[],int n){
        double s=0;
        for(int i=2;i<n;i++) s+=Cross(p[i]-p[1],p[i+1]-p[1]);
        return abs(s/2);
    }

多边形面积

选一个基准点，不停叉积

判断点在多边形内：**射线法**

>  *  引一条向右的水平射线
>  *  注意判断点在多边形上的情况
>  *  逆时针穿过边界，转数+1，顺时针则-1
>  *  注意射线穿过顶点的情况 下算上不算

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

int PointInPolygon(Point p,Point poly[],int n){
        int wn=0;
        for(int i=1;i<=n;i++){
            if(sgn(DisTS(p,poly[i],poly[i%n+1]))==0) return -1;
            int k=sgn(Cross(poly[i%n+1]-poly[i],p-poly[i])),
            d1=sgn(poly[i].y-p.y),d2=sgn(poly[i%n+1].y-p.y);
            if(k>0&&d1<=0&&d2>0) wn++;
            if(k<0&&d2<=0&&d1>0) wn--;
        }
        return (bool)wn;
    }

射线法

**凸包：**

判断凸包：相邻边叉积同号

求凸包：Graham扫描法 我不重写了！细节见代码

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

bool isConvex(Point poly[],int n){
        int last=0,now=0;
        for(int i=1;i<=n;i++){
            now=sgn(Cross(poly[i%n+1]-poly[i],poly[(i+1)%n+1]-poly[i%n+1]));
            if(last==0||now==0||now*last>0) last=now;
            else return false;
        }
        return true;
    }
    
    int ConvexHull(Point p[],int n,Point ch[]){
         //cannot handle repeat point
        sort(p+1,p+1+n);
        int m=0;
        for(int i=1;i<=n;i++){
            while(m>1&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--;
            ch[++m]=p[i];
        }
        int k=m;
        for(int i=n-1;i>=1;i--){
         //n-1
            while(m>k&&sgn(Cross(ch[m]-ch[m-1],p[i]-ch[m-1]))<=0) m--;
            ch[++m]=p[i];
        }
        if(n>1) m--;//the first  rpoint
        return m;
    }

凸包

**旋转卡壳**

>  *  被凸包上被一对平行直线卡住的点对叫**对踵点**
>  *  答案一定出在一对对踵点上
>  *  尝试用通过旋转一对平行直线，枚举所有对踵点
> 
>  *  按逆时针顺序枚举凸包上一条边，则凸包上到该边所在直线最远的点也单调逆时针旋转，相当于用一条直线卡对面一个顶点

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double RotatingCalipers(Point p[],int n){
        if(n==1) return 0;
        if(n==2) return Len(p[1]-p[2]);
        int now=1;
        double ans=0;
        p[n+1]=p[1];
        for(int i=1;i<=n;i++){
            while(sgn(DisTL(p[now],p[i],p[i+1])-DisTL(p[now+1],p[i],p[i+1]))<=0) now=now%n+1;
            ans=max(ans,Len(p[now]-p[i]));
            ans=max(ans,Len(p[now]-p[i+1]));
        }
        return ans;
    }

旋转卡壳

注意特判；变种很多 见题目

**半平面交**

>  *  初始化时加上一个范围巨大的“框”
>  *  每次拿一个新的半平面切割原先的凸集
>  *  保留在新加直线左边的点，删除右边的，有向直线与多边形相交产生的新的点加入到新多边形内

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

void iniPolygon(Point p[],int &n,double inf){
        n=0;
        p[++n]=Point(-inf,-inf);
        p[++n]=Point(inf,-inf);
        p[++n]=Point(inf,inf);
        p[++n]=Point(-inf,inf);
    }
    Point t[N];int tn;
    void CutPolygon(Point p[],int &n,Point a,Point b){
         //get the left of a->b
        tn=0;
        Point c,d,e;
        for(int i=1;i<=n;i++){
            c=p[i],d=p[i%n+1];
            if(sgn(Cross(b-a,c-a))>=0) t[++tn]=c;
            if(isLSI(Line(a,b),Line(c,d))){
                e=LI(Line(a,b),Line(c,d));
                t[++tn]=e;
            }
        }
        n=tn;for(int i=1;i<=n;i++)p[i]=t[i];
    }

半平面交

**最小圆覆盖：**三点确定圆 随即增量法 O(n)

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double minCircleCover(Point p[],int n,Point &c){
        random_shuffle(p+1,p+1+n);
        c=p[1];
        double r=0;
        for(int i=2;i<=n;i++)
            if(sgn(Len(c-p[i])-r)>0){
                c=p[i],r=0;
                for(int j=1;j<i;j++)
                    if(sgn(Len(c-p[j])-r)>0){
                        c=(p[i]+p[j])/2,r=Len(c-p[i]);
                        for(int k=1;k<j;k++)
                            if(sgn(Len(c-p[k])-r)>0){
                                c=Circumcenter(p[i],p[j],p[k]);
                                r=Len(c-p[i]);
                            }
                    }
            }
        return r;
    }

最小圆覆盖

## 6.辛普森积分 ##

用二次函数拟合

我也是看过托马斯微积分的....

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double F(double x){
         return x;} //function
    inline double cal(double l,double r,double fl,double fr,double fm){
        return (fl+fr+4*fm)*(r-l)/6;
    }
    double Simpson(double l,double r,double now,double fl,double fr,double fm){
        double mid=(l+r)/2,flm=F((l+mid)/2),frm=F((mid+r)/2);
        double p=cal(l,mid,fl,fm,flm),q=cal(mid,r,fm,fr,frm);
        if(sgn(now-p-q)==0) return now;
        else return Simpson(l,mid,p,fl,fm,flm)+Simpson(mid,r,q,fm,fr,frm);
    }

自适应辛普森公式

## 7.解析几何 ##

求直线与直线方程交点的神奇方法

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

double a,b,c;
    double f(Point p){
         return a*p.x+b*p.y+c;}
    Point abcLI(Line l){
        double u=abs(f(l.s)),v=abs(f(l.t));
        return Point(l.s.x*v+l.t.x*u,l.s.y*v+l.t.y*u)/(u+v);
    }

ABC

8.**扫描线：**见题目

[ContractedBlock.gif]: /images/20220330/24846dd74e3d4408922d1779c51022f4.png
[ExpandedBlockStart.gif]: /images/20220330/5e238f1f3411435aa750f3ee70d68d72.png