hdu 3662 3D Convex Hull (求三维凸包面数)模版

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. 

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<cstdio>  
#include<vector>  
#include<cstring>  
#include<iostream>  
#include<algorithm>  
#define PR 1e-8  
#define N 510  
using 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;  
}