3D Convex Hull
Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 1266 Accepted Submission(s): 658
Problem Description
There are N points in 3D-space which make up a 3D-Convex hull*. How many faces does the 3D-convexhull have? It is guaranteed that all the points are not in the same plane.
In case you don’t know the definition of convex hull, here we give you a clarification from Wikipedia:
*Convex hull: In mathematics, the convex hull, for a set of points X in a real vector space V, is the minimal convex set containing X.
In case you don’t know the definition of convex hull, here we give you a clarification from Wikipedia:
*Convex hull: In mathematics, the convex hull, for a set of points X in a real vector space V, is the minimal convex set containing X.
Input
There are several test cases. In each case the first line contains an integer N indicates the number of 3D-points (3< N <= 300), and then N lines follow, each line contains three numbers x, y, z (between -10000 and 10000) indicate the 3d-position of a point.
Output
Output the number of faces of the 3D-Convex hull.
Sample Input
7 1 1 0 1 -1 0 -1 1 0 -1 -1 0 0 0 1 0 0 0 0 0 -0.1 7 1 1 0 1 -1 0 -1 1 0 -1 -1 0 0 0 1 0 0 0 0 0 0.1
Sample Output
8 5#include"stdio.h" #include"string.h" #include"iostream" #include"map" #include"string" #include"queue" #include"stack" #include"vector" #include"stdlib.h" #include"algorithm" #include"math.h" #define M 309 #define eps 1e-10 #define inf 0x3f3f3f3f #define mod 1070000009 #define PI acos(-1.0) using namespace std; struct node { double x,y,z,dis; node(){} node(double xx,double yy,double zz):x(xx),y(yy),z(zz){} node operator +(const node p)//向量间求和操作 { return node(x+p.x,y+p.y,z+p.z); } node operator -(const node p)//向量间相减操作 { return node(x-p.x,y-p.y,z-p.z); } node operator *(const node p)//向量间叉乘操作 { return node(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x); } node operator *(const double p)//向量乘以一个数 { return node(x*p,y*p,z*p); } node operator /(const double p)//向量除以一个数 { return node(x/p,y/p,z/p); } double operator ^(const node p)//向量间点乘操作 { return x*p.x+y*p.y+z*p.z; } }; struct threeD_convex_hull//三维凸包 { struct face { int a,b,c; int ok; }; int n;//初始点数 int cnt;//凸包三角形数 node p[M];//初始点 face f[M*8];//凸包三角形 int to[M][M];//点i到j是属于哪个面 double len(node p)//向量的长度 { return sqrt(p.x*p.x+p.y*p.y+p.z*p.z); } double area(node a,node b,node c)//三个点的面积*2 { return len((b-a)*(c-a)); } double volume(node a,node b,node c,node d)//四面体面积*6 { return (b-a)*(c-a)^(d-a); } double ptof(node q,face f)//点与面同向 { node m=p[f.b]-p[f.a]; node n=p[f.c]-p[f.a]; node t=q-p[f.a]; return m*n^t; } void dfs(int q,int cur)//维护凸包，若点q在凸包外则更新凸包 { f[cur].ok=0;//删除当前面，因为此时它在更大的凸包内部 deal(q,f[cur].b,f[cur].a); deal(q,f[cur].c,f[cur].b); deal(q,f[cur].a,f[cur].c); } //因为每个三角形的的三边是按照逆时针记录的，所以把边反过来后对应的就是与ab边共线的另一个面 void deal(int q,int a,int b) { int fa=to[a][b];//与当前面cnt共边的另一个面 face add; if(f[fa].ok)//若fa面目前是凸包的表面则继续 { if(ptof(p[q],f[fa])>eps)//若点q能看到fa面继续深搜fa的三条边，更新新的凸包面 dfs(q,fa); else//当q点可以看到cnt面的同时看不到a，b共边的fa面，则p和a，b点组成一个新的表面三角形 { add.a=b; add.b=a; add.c=q; add.ok=1; to[b][a]=to[a][q]=to[q][b]=cnt; f[cnt++]=add; } } } int same(int s,int t)//判断两个三角形是否共面 { node a=p[f[s].a]; node b=p[f[s].b]; node c=p[f[s].c]; if(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) return 1; return 0; } void make()//构建3D凸包 { cnt=0; if(n<4) return; int sb=1; for(int i=1;i<n;i++)//保证前两个点不共点 { if(len(p[0]-p[i])>eps) { swap(p[1],p[i]); sb=0; break; } } if(sb)return; sb=1; for(int i=2;i<n;i++)//保证前三个点不共线 { if(len((p[1]-p[0])*(p[i]-p[0]))>eps) { swap(p[2],p[i]); sb=0; break; } } if(sb)return; sb=1; for(int i=3;i<n;i++)//保证前四个点不共面 { if(fabs(volume(p[0],p[1],p[2],p[i]))>eps) { swap(p[3],p[i]); sb=0; break; } } if(sb)return; face 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(ptof(p[i],add)>eps) swap(add.c,add.b); to[add.a][add.b]=to[add.b][add.c]=to[add.c][add.a]=cnt; f[cnt++]=add; } for(int i=4;i<n;i++)//倍增法更新凸包 { for(int j=0;j<cnt;j++)//判断每个点是在当前凸包的内部或者外部 { if(f[j].ok&&ptof(p[i],f[j])>eps)//若在外部且看到j面继续 { dfs(i,j); break; } } } int tmp=cnt;//把不是凸包上的面删除即ok=0； cnt=0; for(int i=0;i<tmp;i++) if(f[i].ok) f[cnt++]=f[i]; } double Area()//表面积 { double S=0; if(n==3) { S=area(p[0],p[1],p[2])/2.0; return S; } for(int i=0;i<cnt;i++) S+=area(p[f[i].a],p[f[i].b],p[f[i].c]); return S/2.0; } double Volume()//体积 { double V=0; node mid(0,0,0); for(int i=0;i<cnt;i++) V+=volume(p[f[i].a],p[f[i].b],p[f[i].c],mid); V=fabs(V)/6.0; return V; } int tringleCnt()//凸包表面三角形数目 { return cnt; } int faceCnt()//凸包表面多边形数目 { int num=0; for(int i=0;i<cnt;i++) { int flag=1; for(int j=0;j<i;j++) { if(same(i,j)) { flag=0; break; } } num+=flag; } return num; } double pf_dis(face f,node q)//点到面的距离 { double V=volume(p[f.a],p[f.b],p[f.c],q); double S=area(p[f.a],p[f.b],p[f.c]); return fabs(V/S); } double min_dis(node q)//暴力搜索内部的点q到面的最短距离即体积/面积 { double mini=inf; for(int i=0;i<cnt;i++) { double h=pf_dis(f[i],q); if(mini>h) mini=h; } return mini; } node barycenter()//凸包的重心 { node ret(0,0,0),mid(0,0,0); double sum=0; for(int i=0;i<cnt;i++) { double V=volume(p[f[i].a],p[f[i].b],p[f[i].c],mid); ret=ret+(mid+p[f[i].a]+p[f[i].b]+p[f[i].c])/4.0*V; sum+=V; } ret=ret/sum; return ret; } }hull; int main() { while(scanf("%d",&hull.n)!=-1) { 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.make(); printf("%d\n",hull.faceCnt()); } return 0; }