hdu3662 3D Convex Hull(三维凸包【三维计算几何基本操作)

题目连接

分析:
三维凸包模板

瞧好了基本操作
三维和二维简直不是一个级别的。。。orz

#include<bits/stdc++.h>

using namespace std;

const double eps=1e-8;
const int N=503;
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);}

double Dot(node A,node B) {return A.x*B.x+A.y*B.y+A.z*B.z;}
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);}

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

    int n;            //初始顶点数
    node P[N]; 
    int num;          //凸包表面的三角形数
    face F[8*N];
    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 a,face f) {
        double x=Dot(a-P[f.a],Cross(P[f.b]-P[f.a],P[f.c]-P[f.a]));
        if (x>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[p][b]=vis[a][p]=vis[b][a]=num;
                F[num++]=add;
            }
        }
    }
    //递归搜索所有应该从凸包内删除的面
    void dfs(int p,int now) {
        F[now].ff=0;
        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() {
        face add;

        num=0;
        if (n<4) return;

        //保证前四个点不共面 
        bool flag=1;
        for (int i=1;i<n;i++) {
            if (lenth(P[0]-P[i])>eps) {
                swap(P[1],P[i]);
                flag=0; break;
            }
        }
        if (flag) return;
        flag=1;
        for (int i=2;i<n;i++) {
            if (S(P[0],P[1],P[i])>eps) {
                swap(P[2],P[i]);
                flag=0; break;
            }
        }
        if (flag) return;
        flag=1;
        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=0; 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;   //加面就要维护vis 
            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])){
                    dfs(i,j);
                    break;
                }
        int tmp=num;
        for (int i=num=0;i<tmp;i++)
            if (F[i].ff)
                F[num++]=F[i];
    }

    //表面积 
    double Ssum() {
        double ans=0;
        if (n==3) {
            ans=S(P[0],P[1],P[2]);
            return ans/2.0;
        }
        for (int i=0;i<num;i++)  
            ans=ans+S(P[F[i].a],P[F[i].b],P[F[i].c]);
        return ans/2.0;
    } 
    //体积
    double Vsum() {
        double ans=0;
        node tmp(0,0,0);
        for (int i=0;i<num;i++)
            ans+=V(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(ans/6.0);
    }
    bool same(int s,int t) {
        node &a=P[F[s].a];
        node &b=P[F[s].b];
        node &c=P[F[s].c];
        return fabs(V(a,b,c,P[F[t].a]))<eps &&
               fabs(V(a,b,c,P[F[t].b]))<eps &&
               fabs(V(a,b,c,P[F[t].c]))<eps;
    }
    //表面多边形个数
    int facecnt() {
        int flag,ans=0;
        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;
        }
        ans=ans/all;
        return ans;
    }
    //点到面的距离
    double dis_face(node p,int i) {
        return fabs(V(P[F[i].a],P[F[i].b],P[F[i].c],p)/lenth(Cross(P[F[i].b]-P[F[i].a],P[F[i].c]-P[F[i].a])));
    }
};
CH3D hull;

int main()
{
    while (scanf("%d",&hull.n)!=EOF)
    {
        for (int i=0;i<hull.n;i++)
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
        hull.create();
        printf("%d\n",hull.facecnt());
    }
    return 0;
}
  • 3
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值