计算几何_三维凸包

1.hdoj3662 3D Convex Hull

　　传送：http://acm.hdu.edu.cn/showproblem.php?pid=3662

题意：给出空间n个点，问凸包表面的多边形个数。

分析：rt。

#include<bits/stdc++.h>
using namespace std;
const int maxn=305;
const double eps=1e-8;
int sgn(double x){
    if (fabs(x)<eps) return 0;
    if (x<0) return -1;
    return 1;
}
struct point3{
    double x,y,z;
    point3(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;}
    bool operator ==(const point3 &b)const{
        return sgn(x-b.x)==0 && sgn(y-b.y)==0 && sgn(z-b.z)==0;
    }
    point3 operator -(const point3 &b)const{
        return point3(x-b.x,y-b.y,z-b.z);
    }
    //点乘 
    double operator *(const point3 &b)const{
        return x*b.x+y*b.y+z*b.z; 
    }
    //叉乘
    point3 operator ^(const point3 &b)const{
        return point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
    } 
    double len(){
        return sqrt(x*x+y*y+z*z);
    }
};

struct CH3D{
    struct face{
        int a,b,c; //表示凸包一个面的三个点的编号 
        bool ok; //表示该面是否属于最终的凸包上的面
    };
    int n;  //初始顶点数 
    point3 p[maxn];
    int num;  //凸包表面三角形数
    face f[8*maxn];  //凸包表面的三角形
    int g[maxn][maxn];
    point3 cross(const point3 &a,const point3 &b,const point3 &c){
        return (b-a)^(c-a);
    } 
    //四面体有向面积 *6
    double volume(point3 a,point3 b,point3 c,point3 d){
        return ((b-a)^(c-a))*(d-a);
    } 
    double dblcmp(point3 &P,face &F){
        point3 p1=p[F.b]-p[F.a],p2=p[F.c]-p[F.a],p3=P-p[F.a];
        return (p1^p2)*p3;
    }
    void deal(int x,int a,int b){
        int k=g[a][b];
        face add;
        if (f[k].ok){
            if (dblcmp(p[x],f[k])>eps) dfs(x,k);
            else{
                add.a=b; add.b=a; add.c=x; add.ok=true;
                g[x][b]=g[a][x]=g[b][a]=num;
                f[num++]=add;
            }
        }
    }
    //递归搜索应该从凸包中删除的面 
    void dfs(int p,int now){
        f[now].ok=false;
        deal(p,f[now].b,f[now].a);
        deal(p,f[now].c,f[now].b);
        deal(p,f[now].a,f[now].c);
    }
    //构建三维凸包
    void create(){
        num=0;
        face add;
        
        //保证前四个点不共面
        bool flag=true;
        for (int i=1;i<n;i++)
            if (!(p[0]==p[i])){
                swap(p[1],p[i]);
                flag=false; break;
            }
        if (flag) return ;
        flag=true;
        for (int i=2;i<n;i++)
            if (((p[1]-p[0])^(p[i]-p[0])).len()>eps){
                swap(p[2],p[i]);
                flag=false; break;
            }
        if (flag) return ;
        flag=true;
        for (int i=3;i<n;i++){
            if (fabs(((p[1]-p[0])^(p[2]-p[0]))*(p[i]-p[0]))>eps){
                swap(p[3],p[i]);
                flag=false; 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.ok=true;
            if (dblcmp(p[i],add)>0) swap(add.b,add.c);
            g[add.a][add.b]=g[add.b][add.c]=g[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].ok && dblcmp(p[i],f[j])>eps){
                    dfs(i,j);
                    break;
                }
        }
        int tmp=num; num=0;
        for (int i=0;i<tmp;i++) if (f[i].ok) f[num++]=f[i];
    } 
    bool same(int s,int t){
        point3 &a=p[f[s].a],&b=p[f[s].b],&c=p[f[s].c];
        return fabs(volume(a,b,c,p[f[t].a]))<eps && fabs(volume(a,b,c,p[f[t].b]))<eps
                && fabs(volume(a,b,c,p[f[t].c]))<eps;
    }
    //凸包上多边形面积
    int polygon(){
        int res=0;
        for (int i=0;i<num;i++){
            int flag=1;
            for (int j=0;j<i;j++)
                if (same(i,j)){flag=0;break;}
            res+=flag;
        } 
        return res;
    }
};
CH3D C;
int main(){
    int n; double x,y,z;
    while (cin >> n){
        C.n=n;
        for (int i=0;i<n;i++) cin >> C.p[i].x >> C.p[i].y >> C.p[i].z;
        C.create();
        cout << C.polygon() << endl;
    } 
    return 0;
} 

 2.hdoj4266 The Worm in the Apple

传送：http://acm.hdu.edu.cn/showproblem.php?pid=4266

题意：对于每个询问，问到凸包表面的最小距离。

分析：rt。

#include<bits/stdc++.h>
using namespace std;
const int maxn=1005;
const double eps=1e-8;
const double inf=1e20;
int sgn(double x){
    if (fabs(x)<eps) return 0;
    if (x<0) return -1;
    return 1;
}
struct point3{
    double x,y,z;
    point3(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;}
    bool operator ==(const point3 &b)const{
        return sgn(x-b.x)==0 && sgn(y-b.y)==0 && sgn(z-b.z)==0;
    }
    point3 operator -(const point3 &b)const{
        return point3(x-b.x,y-b.y,z-b.z);
    }
    //点乘 
    double operator *(const point3 &b)const{
        return x*b.x+y*b.y+z*b.z; 
    }
    //叉乘
    point3 operator ^(const point3 &b)const{
        return point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
    } 
    double len(){
        return sqrt(x*x+y*y+z*z);
    }
};

struct CH3D{
    struct face{
        int a,b,c; //表示凸包一个面的三个点的编号 
        bool ok; //表示该面是否属于最终的凸包上的面
    };
    int n;  //初始顶点数 
    point3 p[maxn];
    int num;  //凸包表面三角形数
    face f[8*maxn];  //凸包表面的三角形
    int g[maxn][maxn];
    point3 cross(const point3 &a,const point3 &b,const point3 &c){
        return (b-a)^(c-a);
    } 
    //四面体有向面积 *6
    double volume(point3 a,point3 b,point3 c,point3 d){
        return ((b-a)^(c-a))*(d-a);
    } 
    double dblcmp(point3 &P,face &F){
        point3 p1=p[F.b]-p[F.a],p2=p[F.c]-p[F.a],p3=P-p[F.a];
        return (p1^p2)*p3;
    }
    void deal(int x,int a,int b){
        int k=g[a][b];
        face add;
        if (f[k].ok){
            if (dblcmp(p[x],f[k])>eps) dfs(x,k);
            else{
                add.a=b; add.b=a; add.c=x; add.ok=true;
                g[x][b]=g[a][x]=g[b][a]=num;
                f[num++]=add;
            }
        }
    }
    //递归搜索应该从凸包中删除的面 
    void dfs(int p,int now){
        f[now].ok=false;
        deal(p,f[now].b,f[now].a);
        deal(p,f[now].c,f[now].b);
        deal(p,f[now].a,f[now].c);
    }
    //构建三维凸包
    void create(){
        num=0;
        face add;
        
        //保证前四个点不共面
        bool flag=true;
        for (int i=1;i<n;i++)
            if (!(p[0]==p[i])){
                swap(p[1],p[i]);
                flag=false; break;
            }
        if (flag) return ;
        flag=true;
        for (int i=2;i<n;i++)
            if (((p[1]-p[0])^(p[i]-p[0])).len()>eps){
                swap(p[2],p[i]);
                flag=false; break;
            }
        if (flag) return ;
        flag=true;
        for (int i=3;i<n;i++){
            if (fabs(((p[1]-p[0])^(p[2]-p[0]))*(p[i]-p[0]))>eps){
                swap(p[3],p[i]);
                flag=false; 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.ok=true;
            if (dblcmp(p[i],add)>0) swap(add.b,add.c);
            g[add.a][add.b]=g[add.b][add.c]=g[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].ok && dblcmp(p[i],f[j])>eps){
                    dfs(i,j);
                    break;
                }
        }
        int tmp=num; num=0;
        for (int i=0;i<tmp;i++) if (f[i].ok) f[num++]=f[i];
    } 
    double ptoface(point3 pp,int i){
        double tmp1=fabs(volume(p[f[i].a],p[f[i].b],p[f[i].c],pp));
        double tmp2=((p[f[i].b]-p[f[i].a])^(p[f[i].c]-p[f[i].a])).len();
        return tmp1/tmp2;
    }
    double slove(point3 pp){
        double mi=inf;
        for (int i=0;i<num;i++) mi=min(mi,ptoface(pp,i));
        return mi;
    }    
};
CH3D C;
int main(){
    int n,q; point3 pp; 
    while (cin >> n && n){
        C.n=n;
        for (int i=0;i<n;i++) cin >> C.p[i].x >> C.p[i].y >> C.p[i].z;
        C.create();
        cin >> q; 
        while (q--){
            cin >> pp.x >> pp.y >> pp.z;
            double ans=inf;
            double tmp=C.slove(pp);
            ans=min(ans,tmp);
            printf("%.4f
",(ans));
        }
    }
    return 0;
}

 3.hdoj4273 Rescue

传送：http://acm.hdu.edu.cn/showproblem.php?pid=4273

题意：给出一个三维凸包，求凸包重心到凸包表面的最短距离 。

分析：rt。

#include<bits/stdc++.h>
using namespace std;
const int maxn=305;
const double eps=1e-8;
const double inf=1e20;
int sgn(double x){
    if (fabs(x)<eps) return 0;
    if (x<0) return -1;
    return 1;
}
struct point3{
    double x,y,z;
    point3(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;}
    bool operator ==(const point3 &b)const{
        return sgn(x-b.x)==0 && sgn(y-b.y)==0 && sgn(z-b.z)==0;
    }
    point3 operator +(const point3 &b)const{
        return point3(x+b.x,y+b.y,z+b.z);
    }
    point3 operator -(const point3 &b)const{
        return point3(x-b.x,y-b.y,z-b.z);
    }
    point3 operator *(const double k)const{
        return point3(x*k,y*k,z*k);
    }
    point3 operator /(const double k)const{
        return point3(x/k,y/k,z/k);
    } 
    //点乘 
    double operator *(const point3 &b)const{
        return x*b.x+y*b.y+z*b.z; 
    }
    //叉乘
    point3 operator ^(const point3 &b)const{
        return point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);
    } 
    double len(){
        return sqrt(x*x+y*y+z*z);
    }
};

struct CH3D{
    struct face{
        int a,b,c; //表示凸包一个面的三个点的编号 
        bool ok; //表示该面是否属于最终的凸包上的面
    };
    int n;  //初始顶点数 
    point3 p[maxn];
    int num;  //凸包表面三角形数
    face f[8*maxn];  //凸包表面的三角形
    int g[maxn][maxn];
    point3 cross(const point3 &a,const point3 &b,const point3 &c){
        return (b-a)^(c-a);
    } 
    //四面体有向面积 *6
    double volume(point3 a,point3 b,point3 c,point3 d){
        return ((b-a)^(c-a))*(d-a);
    } 
    double dblcmp(point3 &P,face &F){
        point3 p1=p[F.b]-p[F.a],p2=p[F.c]-p[F.a],p3=P-p[F.a];
        return (p1^p2)*p3;
    }
    void deal(int x,int a,int b){
        int k=g[a][b];
        face add;
        if (f[k].ok){
            if (dblcmp(p[x],f[k])>eps) dfs(x,k);
            else{
                add.a=b; add.b=a; add.c=x; add.ok=true;
                g[x][b]=g[a][x]=g[b][a]=num;
                f[num++]=add;
            }
        }
    }
    //递归搜索应该从凸包中删除的面 
    void dfs(int p,int now){
        f[now].ok=false;
        deal(p,f[now].b,f[now].a);
        deal(p,f[now].c,f[now].b);
        deal(p,f[now].a,f[now].c);
    }
    //构建三维凸包
    void create(){
        num=0;
        face add;
        
        //保证前四个点不共面
        bool flag=true;
        for (int i=1;i<n;i++)
            if (!(p[0]==p[i])){
                swap(p[1],p[i]);
                flag=false; break;
            }
        if (flag) return ;
        flag=true;
        for (int i=2;i<n;i++)
            if (((p[1]-p[0])^(p[i]-p[0])).len()>eps){
                swap(p[2],p[i]);
                flag=false; break;
            }
        if (flag) return ;
        flag=true;
        for (int i=3;i<n;i++){
            if (fabs(((p[1]-p[0])^(p[2]-p[0]))*(p[i]-p[0]))>eps){
                swap(p[3],p[i]);
                flag=false; 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.ok=true;
            if (dblcmp(p[i],add)>0) swap(add.b,add.c);
            g[add.a][add.b]=g[add.b][add.c]=g[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].ok && dblcmp(p[i],f[j])>eps){
                    dfs(i,j);
                    break;
                }
        }
        int tmp=num; num=0;
        for (int i=0;i<tmp;i++) if (f[i].ok) f[num++]=f[i];
    } 
    //求凸包重心 
    point3 barycenter(){
        point3 ans=point3(0,0,0),o=point3(0,0,0);
        double all=0;
        for (int i=0;i<num;i++){
            double vol=volume(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;
        }
        ans=ans/all;
        return ans;
    } 
    //点到平面距离 
    double ptoface(point3 pp,int i){
        double tmp1=fabs(volume(p[f[i].a],p[f[i].b],p[f[i].c],pp));
        double tmp2=((p[f[i].b]-p[f[i].a])^(p[f[i].c]-p[f[i].a])).len();
        return tmp1/tmp2;
    }
};
CH3D C;
int main(){
    int n; double x,y,z;
    while (cin >> n){
        C.n=n;
        for (int i=0;i<n;i++) cin >> C.p[i].x >> C.p[i].y >> C.p[i].z;
        C.create();
        point3 pp=C.barycenter();
        double ans=inf;
        for (int i=0;i<C.num;i++)
            ans=min(ans,C.ptoface(pp,i));
        printf("%.3f
",ans);
    }
    return 0;
}