首页 > 代码库 > HDU 3662

HDU 3662

求的是凸包有多少个面。

注意,求的是面。这就需要把同一个面的三角形合并。只需判断两个三角形的法向量是否同向平行。

/*增量法求凸包。选取一个四面体,同时把它各面的方向向量向外,增加一个点时,若该点与凸包上的某些面的方向向量在同一侧,则去掉那些面,并使某些边与新增点一起连成新的凸包上的面。 */ #include <iostream>#include <cstring>#include <cstdio>#include <algorithm>#include <cmath> using namespace std;const int MAXN=350;const double eps=1e-8;struct point {    double x,y,z;};struct face {    int a,b,c;    bool ok;};int n;  //初始点数 point p[MAXN]; //空间点int trianglecnt; //凸包上三角形数face tri[6*MAXN]; //凸包上被创建的三角形int vis[MAXN][MAXN]; //点i到点j是属于哪一个三角形。此处是有方向bool vit[MAXN*6];point operator -(const point &x, const point &y){    point ret;    ret.x=x.x-y.x; ret.y=x.y-y.y; ret.z=x.z-y.z;    return ret;} point operator * (const point &u,const point &v){  //叉积     point ret;    ret.x=u.y*v.z-u.z*v.y;    ret.y=u.z*v.x-u.x*v.z;    ret.z=u.x*v.y-u.y*v.x;    return ret;} double  operator ^(const point &u,const point &v){    return (u.x*v.x+u.y*v.y+u.z*v.z);} double dist(point t){    return sqrt(t.x*t.x+t.y*t.y+t.z*t.z);} double ptoplane(point &tmp,face &f){    //若结果大于0,证明点面的同向,即法向量方向     point m=p[f.b]-p[f.a]; point n=p[f.c]-p[f.a];    point t=tmp-p[f.a];    return (m*n)^t;} double farea(point a,point b,point c ){    point t1=a-c; point t2=b-c;    return fabs(dist(t1*t2));}void dfs(int pt, int ct);void deal(int pt,int a,int b){    int f=vis[a][b];   //所属三角形,即原来的ab。     face add;    if(tri[f].ok){        if((ptoplane(p[pt],tri[f]))>eps) dfs(pt,f);   //若点同样在该f三角形方向一侧,继续调整         else {            add.a=b; add.b=a; add.c=pt; add.ok=1;            vis[pt][b]=vis[a][pt]=vis[b][a]=trianglecnt;            tri[trianglecnt++]=add;        }    }} void dfs(int pt, int ct){    tri[ct].ok=0;   //去掉该面     deal(pt,tri[ct].b,tri[ct].a);   //因为有向边ab所属三角形去掉,则反方向边必定属于另一个三角形.     deal(pt,tri[ct].c,tri[ct].b);    deal(pt,tri[ct].a,tri[ct].c);} void construct (){    int i,j;    trianglecnt=0;    if(n<4) return ; //不可能构成一个多面体    bool tmp=true;     for(i=1;i<n;i++){    //不共点两点         if(dist(p[0]-p[i])>eps){            swap(p[1],p[i]); tmp=false; break;        }    }    if(tmp) return ;    tmp=true;    for(i=2;i<n;i++){   //不共线         if(dist((p[0]-p[1])*(p[1]-p[i]))>eps){            swap(p[2],p[i]); tmp=false; break;        }    }    if(tmp) return ;    tmp=true;    for(i=3;i<n;i++){   //四点不共面K         if(fabs((p[0]-p[1])*(p[1]-p[2])^(p[0]-p[i]))>eps){            swap(p[3],p[i]); tmp=false; break;        }    }    if(tmp) return ;    face add;    for(i=0;i<4;i++){   //使各三角形的方向向量向外,同时记录下三角形的序号         add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=1;  //等于1表示在凸包上         if(ptoplane(p[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(i=4;i<n;i++){   //构建凸包         for(j=0;j<trianglecnt;j++){            if(tri[j].ok&&(ptoplane(p[i],tri[j]))>eps){  //增加点可见该平,即在面方向一侧                 dfs(i,j); break;            }        }    }    int cnt=trianglecnt;    trianglecnt=0;    for(i=0;i<cnt;i++){    //只有ok为1的才属于凸包上的三角形         if(tri[i].ok){            tri[trianglecnt++]=tri[i];        }    }}void clean(){	memset(vit,false,sizeof(vit));	int counted=0;	for(int i=0;i<trianglecnt;i++){		if(vit[i]) continue;		vit[i]=true;		point vect1=(p[tri[i].b]-p[tri[i].a])*(p[tri[i].c]-p[tri[i].a]);		counted++;		for(int j=i+1;j<trianglecnt;j++){			point vect2=(p[tri[j].b]-p[tri[j].a])*(p[tri[j].c]-p[tri[j].a]);			if(dist(vect1*vect2)<eps&&(vect1^vect2)>0)			vit[j]=true;		}	}	printf("%d\n",counted);}int main(){    while(scanf("%d",&n)!=EOF){        memset(vis,-1,sizeof(vis));        for(int i=0;i<n;i++)        scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);        construct();        clean();    }}