Hdu 3662 3D Convex Hull（三维凸包）

来源：互联网发布：纳网科技域名续缴费编辑：程序博客网时间：2024/04/30 12:37
题目地址：http://acm.hdu.edu.cn/showproblem.php?pid=3662
思路：三维凸包模板。
#include<cstdio>#include<vector>#include<cstring>#include<iostream>#include<algorithm>#define PR 1e-8#define N 510using namespace std;struct TPoint{    double x,y,z;    TPoint() {}    TPoint(double _x,double _y,double _z):x(_x),y(_y),z(_z) {}    TPoint operator - (const TPoint p)    {        return TPoint(x-p.x,y-p.y,z-p.z);    }    TPoint operator * (const TPoint p)    {        return TPoint(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);    }    double operator ^ (const TPoint p)    {        return x*p.x+y*p.y+z*p.z;    }//点积};struct fac{    int a,b,c;//凸包一个面上的三个点的编号    bool ok;//该面是否是最终凸包中的面};struct T3dull{    int n;//初始点数    TPoint ply[N];//初始点    int trianglecnt;//凸包上三角形数    fac tri[N];//凸包三角形    int vis[N][N];//点i到点j是属于哪个面    double dist(TPoint a)    {        return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);    }//两点长度    double area(TPoint a,TPoint b,TPoint c)    {        return dist((b-a)*(c-a));    }//三角形面积*2    double volume(TPoint a,TPoint b,TPoint c,TPoint d)    {        return (b-a)*(c-a)^(d-a);    }//四面体有向体积*6    double ptoplane(TPoint &p,fac &f)//正：点在面同向    {        TPoint m=ply[f.b]-ply[f.a],n=ply[f.c]-ply[f.a],t=p-ply[f.a];        return (m*n)^t;    }    void deal(int p,int a,int b)    {        int f=vis[a][b];        fac add;        if(tri[f].ok)        {            if((ptoplane(ply[p],tri[f]))>PR) dfs(p,f);            else            {                add.a=b,add.b=a,add.c=p,add.ok=1;                vis[p][b]=vis[a][p]=vis[b][a]=trianglecnt;                tri[trianglecnt++]=add;            }        }    }    void dfs(int p,int cnt)//维护凸包，如果点p在凸包外更新凸包    {        tri[cnt].ok=0;        deal(p,tri[cnt].b,tri[cnt].a);        deal(p,tri[cnt].c,tri[cnt].b);        deal(p,tri[cnt].a,tri[cnt].c);    }    bool same(int s,int e)//判断两个面是否为同一面    {        TPoint a=ply[tri[s].a],b=ply[tri[s].b],c=ply[tri[s].c];        return fabs(volume(a,b,c,ply[tri[e].a]))<PR&&fabs(volume(a,b,c,ply[tri[e].b]))<PR&&fabs(volume(a,b,c,ply[tri[e].c]))<PR;    }    void construct()//构建凸包    {        int i,j;        trianglecnt=0;        if(n<4) return ;        bool tmp=1;        for(int i=1; i<n; i++) //前两点不共点        {            if((dist(ply[0]-ply[i]))>PR)            {                swap(ply[1],ply[i]);                tmp=0;                break;            }        }        if(tmp) return ;        tmp=1;        for(i=2; i<n; i++) //前三点不共线        {            if((dist((ply[0]-ply[1])*(ply[1]-ply[i])))>PR)            {                swap(ply[2],ply[i]);                tmp=0;                break;            }        }        if(tmp) return ;        tmp=1;        for(i=3; i<n; i++) //前四点不共面        {            if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>PR)            {                swap(ply[3],ply[i]);                tmp=0;                break;            }        }        if(tmp) return ;        fac add;        for(int i=0; i<4; i++) //构建初始四面体        {            add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.ok=1;            if((ptoplane(ply[i],add))>0)                swap(add.b,add.c);            vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=trianglecnt;            tri[trianglecnt++]=add;        }        for(int i=4; i<n; i++) //构建更新凸包        {            for(j=0; j<trianglecnt; j++)            {                if(tri[j].ok&&(ptoplane(ply[i],tri[j]))>PR)                {                    dfs(i,j);                    break;                }            }        }        int cnt=trianglecnt;        trianglecnt=0;        for(i=0; i<cnt; i++)        {            if(tri[i].ok)                tri[trianglecnt++]=tri[i];        }    }    double area()//表面积    {        double ret=0;        for(int i=0; i<trianglecnt; i++)            ret+=area(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);        return ret/2.0;    }    double volume()//体积    {        TPoint p(0,0,0);        double ret=0;        for(int i=0; i<trianglecnt; i++)            ret+=volume(p,ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);        return fabs(ret/6);    }    int facetri()//表面三角形数    {        return trianglecnt;    }    int facepolygon()//表面多边形数    {        int ans=0,i,j,k;        for(i=0; i<trianglecnt; i++)        {            for(j=0,k=1; j<i; j++)            {                if(same(i,j))                {                    k=0;                    break;                }            }            ans+=k;        }        return ans;    }} hull;int main(){    while(~scanf("%d",&hull.n))    {        int i;        for(i=0; i<hull.n; i++)            scanf("%lf%lf%lf",&hull.ply[i].x,&hull.ply[i].y,&hull.ply[i].z);        hull.construct();        printf("%d\n",hull.facepolygon());    }    return 0;}
0 0