计算几何全纪录

最新推荐文章于 2023-03-09 18:35:30 发布

Coco_T_

最新推荐文章于 2023-03-09 18:35:30 发布

阅读量294

点赞数

分类专栏：全纪录系列文章标签：计算几何

本文链接：https://blog.csdn.net/wu_tongtong/article/details/79737447

版权

全纪录系列专栏收录该内容

19 篇文章 0 订阅

订阅专栏

计算几何，分为二维和三维

二维计算几何

二维计算几何基础入门

二维计算几何中，应用比较广泛的就是凸包和半平面交

凸包

几何意义上的凸包比较好理解
但是很多dp，分治题目中状态转移的最优点满足某种单调性，也就是说优秀转移点在凸包上
这种性质可能不好判断，不过斜率优化dp就是凸包和dp的完美结合

经典例题：凸包
经典例题：线段树分治+凸包
经典例题：分块+凸包+三分

注意栈的边界条件
实际上，Cross是不满足交换律的： $Cross(u,w)=-Cross(w,u)$
$Cross(p[i]-p[S[top-1]],p[S[top]]-p[S[top-1]])$
但是不管顺序是什么，上下两个while里的Cross一样，不等号方向一样就可以了

const int N=100010;
const double eps=1e-8;
struct node{
    double x,y;
    node(double xx=0,double yy=0) {
        x=xx;y=yy;
    }
}p[N];

node operator +(const node &A,const node &B) {return node(A.x+B.x,A.y+B.y);}
node operator -(const node &A,const node &B) {return node(A.x-B.x,A.y-B.y);}
node operator *(const node &A,const double &B) {return node(A.x*B,A.y*B);}
node operator /(const node &A,const double &B) {return node(A.x/B,A.y/B);}

int S[N],top=0,n;

int dcmp(double x) {
    if (fabs(x)<eps) return 0;
    else if (x>0) return 1;
    else return -1;
}

int cmp(const node &A,const node &B) {
    if (fabs(A.x-B.x)<eps) return A.y<B.y;
    else return A.x<B.x;
}

double Cross(node A,node B) {return A.x*B.y-A.y*B.x;}
double Dot(node A,node B) {return A.x*B.x+A.y*B.y;}

void TuB() {
    sort(p+1,p+1+n,cmp);
    for (int i=1;i<=n;i++) {
        while (top>1 && dcmp(Cross(p[i]-p[S[top-1]],p[S[top]]-p[S[top-1]]))<=0) top--;
        //上凸壳 top>1 
        S[++top]=i;
    }
    int k=top;
    for (int i=n-1;i>=1;i--) {
        while (top>k && dcmp(Cross(p[i]-p[S[top-1]],p[S[top]]-p[S[top-1]]))<=0) top--;
        //下凸壳 top>k 
        S[++top]=i;
    }
    if (n>1) top--;   //
}

旋xuán,xuàn转zhuǎn,zhuàn卡qiǎ,kǎ壳ké,qiào

凸包的扩展操作：求最远点对
求出凸包之后维护两条平行线卡住凸包，维护平行线上点对之间的最远距离即可

旋转卡壳讲解
 经典例题：旋转卡壳（难）

感性认知，我们需要求点到直线的距离： $dis(a,p_1,p_2)=fabs({Cross(a-p_1,p_2-p_1) \over len(p_2-p_1)})$
但是我们不需要这么麻烦
因为叉积的定义： $Cross(a,b)=absinθ$ ，直接求叉积即可
（不过我觉得还是保险一点求dis好）

double len(node a) {return sqrt(Dot(a,a));}
double dis(node a,node p1,node p2) {             //点到直线的距离 
    return fabs(Cross(a-p1,p2-p1))/len(p2-p1);
}

node q[N];

void solve() {
    TuB();
    for (int i=1;i<=top;i++) q[i-1]=p[S[i]];
    int now=1;
    double ans=0;                     //最远点对距离 
    for (int i=0;i<top;i++) {         //枚举每一条边 
        int j=(i+1)%top;
        while ( dis(q[now],q[i],q[j]) < dis(q[(now+1)%top],q[i],q[j]) ) now=(now+1)%top;
        ans=max(ans,len(q[now]-q[i]));
    }
}

半平面交

半平面交的重要程度和凸包可以说是不分上下，而且变化更多

二维平面上的半平面交，没什么可说的
一些问题会出现一次函数的形式，可以把一次函数视为平面中的直线，最优解一般就是在半平面交的顶点上
同样的，一些双元不等式也可以用半平面交完成（判断是否有解）
（不等式不是用单纯形吗？话是这么说没错啦，不过单纯形复杂度高啊）

经典例题：二分+半平面+直线平移
经典例题：半平面交（一）
经典例题：半平面交（二）铁人双项
 经典例题：问题转化=>半平面交+单调栈
经典例题：问题转化=>半平面交（难）

在二维平面上，只有点按照顺序（顺时针，逆时针）给出我们才能安心的半平面

有时候，半平面交的题目中会遇到直线平移的问题

已知平移距离为 $d$ ，则前后两条直线 $c$ 值之差为

| c 1 - c 2 | = d * a 2 + b 2 - - - - - - \sqrt

$|c_1-c_2|={d*\sqrt{a^2+b^2}}$

直线向ta的左边平移，则 $c(l')=c(l)+{d*\sqrt{a^2+b^2}}$
直线向ta的右边平移，则 $c(l')=c(l)-{d*\sqrt{a^2+b^2}}$

const int N=100010;

struct node{
    double x,y;
    node (double xx=0,double yy=0) {
        x=xx;y=yy;
    }
}po[N],p[N],q[N];

int n,mm;
double a,b,c;

void getline(node A,node B) {
    a=B.y-A.y;
    b=A.x-B.x;
    c=A.y*B.x-A.x*B.y;
}

node insert(node A,node B) {
    double u=fabs(a*A.x+b*A.y+c);    //fabs
    double v=fabs(a*B.x+b*B.y+c);
    node ans;
    ans.x=(u*B.x+v*A.x)/(u+v);
    ans.y=(u*B.y+v*A.y)/(u+v);
    return ans;
}

void cut() {
    int cnt=0;
    for (int i=1;i<=m;i++) {
        if (a*p[i].x+b*p[i].y+c<=0) q[++cnt]=p[i];   //逆时针给点,符合要求的点在直线的左边 <0 
        else {
            if (a*p[i-1].x+b*p[i-1].y+c<0) 
                q[++cnt]=insert(p[i-1],p[i]);
            if (a*p[i+1].x+b*p[i+1].y+c<0) 
                q[++cnt]=insert(p[i],p[i+1]);
        }
    }
    for (int i=1;i<=cnt;i++) p[i]=q[i];
    m=cnt;
    p[0]=p[m]; p[m+1]=p[1];
}

void solve() {
    for (int i=1;i<=n;i++) p[i]=po[i];
    po[0]=po[n]; po[n+1]=po[1];
    p[0]=p[n]; p[n+1]=p[1];
    m=n;
    for (int i=1;i<=n;i++) {
        getline(po[i],po[i+1]);
        cut();
        if (!m) break;             //没有凸核 
    }
}

最小圆覆盖

最小圆覆盖讲解
最小圆覆盖：用最小面积的圆覆盖平面上n个点
应用可能比较窄，但是其中渗透的二维计算几何操作（向量旋转，直线交点）还是不错的

不要忘了random_shuffle
有一个小技巧，要取fabs我们可以直接套一个lenth

需要注意的是，我们在求直线交点的时候，任何一个变量的顺序都不要错（mmp记住就好了）：

node jiao(node P,node v,node Q,node w) {
    node u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}

const double Pi=acos(-1.0);

double lenth(node a) {return sqrt(Dot(a,a));}

node jiao(node P,node v,node Q,node w) {
    node u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}

node rotate(node P,double a) {
    return node(P.x*cos(a)-P.y*sin(a),P.x*sin(a)+P.y*cos(a));
}

node get(node a,node b,node c) {
    node P=(a+b)/2;
    node Q=(a+c)/2;
    node v=rotate(a-b,Pi/2.0);    //b-a
    node w=rotate(a-c,Pi/2.0);    //c-a
    if (dcmp(Cross(v,w))==0) {
        if (dcmp(lenth(a-b)-lenth(a-c)-lenth(c-b))==0) 
            return (a+b)/2;
        if (dcmp(lenth(a-c)-lenth(a-b)-lenth(b-c))==0) 
            return (a+c)/2;
        if (dcmp(lenth(b-c)-lenth(b-a)-lenth(a-c))==0) 
            return (b+c)/2;
    }
    return jiao(P,v,Q,w);
}

void min_circular() {
    random_shuffle(p+1,p+1+n);    //随机化 保证复杂度 
    node o=p[1];
    double r=0;
    for (int i=2;i<=n;i++) 
        if (dcmp(lenth(p[i]-o)-r)>0) 
        {
            o=p[i],r=0.0;
            for (int j=1;j<i;j++)
                if (dcmp(lenth(p[j]-o)-r)>0)
                {
                    o=(p[i]+p[j])/2,r=lenth(p[i]-o);
                    for (int k=1;k<j;k++) 
                        if (dcmp(lenth(p[k]-o)-r)>0) 
                        {
                            o=get(p[i],p[j],p[k]);
                            r=lenth(p[i]-o);
                        }
                }
        }
}

三维计算几何

三维计算几何基础入门

三维凸包讲解
三维计算几何中比wei较yi重xue要guo的就是三维凸包了
三维凸包的基本操作可以说就囊括了三维计算几何的重要操作了

const int N=500;
const double eps=1e-7;

struct node{
    double x,y,z;
    node (double xx=0,double yy=0,double zz=0) {
        x=xx;y=yy;z=zz;
    }
}; 

node operator +(const node &A,const node &B) {return node(A.x+B.x,A.y+B.y,A.z+B.z);}
node operator -(const node &A,const node &B) {return node(A.x-B.x,A.y-B.y,A.z-B.z);}
node operator *(const node &A,const double &B) {return node(A.x*B,A.y*B,A.z*B);}
node operator /(const node &A,const double &B) {return node(A.x/B,A.y/B,A.z/B);}

node Cross(node A,node B) {return node(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);}
double Dot(node A,node B) {return A.x*B.x+A.y*B.y+A.z*B.z;}

struct CH3D{
    struct face{
        int a,b,c;
        //凸包一个面上的三个点编号 
        bool ff;
        //该面是否属于最终凸包上的面
    };

    int n,num;
    node P[N];
    face F[N<<3];
    int vis[N][N];

    //向量长度
    double lenth(node a) {
        return sqrt(Dot(a,a));
    }
    //三角形面积*2 
    double S(node a,node b,node c) {
        return lenth(Cross(b-a,c-a));         //注意要取lenth
    }
    //四面体有向体积*6  混合积 
    double V(node a,node b,node c,node d) {
        return Dot(d-a,Cross(b-a,c-a));
    }
    //面是否可见
    int cansee(node p,face f) {
        double t=Dot(p-P[f.a],Cross(P[f.b]-P[f.a],P[f.c]-P[f.a]));
        if (t>eps) return 1;
        else return 0;
    } 

    //搜索与该边相邻的另一个平面
    void deal(int p,int a,int b) {
        int f=vis[a][b];
        face add;
        if (F[f].ff) {
            if (cansee(P[p],F[f])) dfs(p,f);
            else {
                add.a=b;
                add.b=a;
                add.c=p;
                add.ff=1;
                vis[b][a]=vis[a][p]=vis[p][b]=num;
                F[num++]=add;
            }
        }
    }
    //递归搜索所有应该从凸包内删除的面
    void dfs(int i,int j) {
        F[j].ff=0;                //把面的方向改变一下 
        deal(i,F[j].b,F[j].a);
        deal(i,F[j].c,F[j].b);
        deal(i,F[j].a,F[j].c);
    }
    //构建三维凸包
    void create() {
        face add;
        num=0;

        if (n<4) return;

        int flag=0;
        for (int i=1;i<n;i++) 
            if (lenth(P[0]-P[i])>eps) {
                swap(P[1],P[i]);
                flag=1;
                break;
            }
        if (!flag) return;
        flag=0;
        for (int i=2;i<n;i++)
            if (S(P[0],P[1],P[i])>eps) {
                swap(P[2],P[i]);
                flag=1;
                break;
            }
        if (!flag) return;
        flag=0;
        for (int i=3;i<n;i++) 
            if (fabs(V(P[0],P[1],P[2],P[i]))>eps) {
                swap(P[3],P[i]);
                flag=1;
                break;
            }
        if (!flag) return;

        for (int i=0;i<4;i++) {
            add.a=(i+1)%4;
            add.b=(i+2)%4;
            add.c=(i+3)%4;
            add.ff=1;
            if (cansee(P[i],add)) swap(add.b,add.c);
            vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=num;
            F[num++]=add;
        }

        for (int i=4;i<n;i++)
            for (int j=0;j<num;j++)
                if (F[j].ff&&cansee(P[i],F[j])) {   //F[j].ff
                    dfs(i,j);
                    break;
                }
        int tmp=num; num=0;
        for (int i=0;i<tmp;i++)
            if (F[i].ff)
                F[num++]=F[i];
    }

    //表面积 
    double Ssum() {
        if (n==3) {
            return S(P[0],P[1],P[2])/2.0;
        }
        double ans=0;
        for (int i=0;i<num;i++) 
            ans+=S(P[F[i].a],P[F[i].b],P[F[i].c]);
        return ans/2.0;
    }
    //体积
    double Vsum() {
        double ans;
        node o(0,0,0);
        for (int i=0;i<num;i++)
            ans+=V(o,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(ans)/6.0;
    }
    bool same(int i,int j) {
        node a=P[F[i].a];
        node b=P[F[i].b];
        node c=P[F[i].c];
        return fabs(V(a,b,c,P[F[j].a]))<eps &&
               fabs(V(a,b,c,P[F[j].b]))<eps &&
               fabs(V(a,b,c,P[F[j].c]))<eps;
    }
    //表面多边形个数
    int face_cnt() {
        int ans=0,flag;
        for (int i=0;i<num;i++) {
            flag=1;
            for (int j=0;j<i;j++)
                if (same(i,j)) {
                    flag=0;
                    break;
                }
            ans+=flag;
        }
        return ans;
    }
    //三维凸包重心
    node centre() {
        node ans(0,0,0),o(0,0,0);
        double all=0;
        for (int i=0;i<num;i++) {
            double vol=V(o,P[F[i].a],P[F[i].b],P[F[i].c]);
            ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
            all+=vol;
        }
        return ans/all;
    }
    //点到面的距离
    double dis(node p,face f) {
        return fabs(V(p,P[f.a],P[f.b],P[f.c])/S(P[f.a],P[f.b],P[f.c]));
    }
};