三维凸包

三维凸包模板 hdu4266

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>

using namespace std;
const int MAXN = 1100;
const double eps = 1e-9;
struct Point
{
     double x,y,z;
     Point(){}
     Point(double _x,double _y,double _z):x(_x),y(_y),z(_z){}
     Point operator -(const Point &P)
     {
           return Point(x-P.x, y-P.y, z-P.z);
     }  
     Point operator ^ (const Point &P)//叉积 
     {
           return Point(y*P.z - z * P.y, z * P.x - x * P.z, x*P.y - y * P.x);
     }
     double operator * (const Point &P)//点积 
     {
           return (x * P.x + y * P.y + z * P.z); 
     } 
     
};
struct CH3D
{
       struct face
       {
            int a,b,c;
            bool ok ;// 判断是否是凸包的面 
            face(){}
            face(int _a,int _b,int _c, bool _ok):a(_a),b(_b),c(_c),ok(_ok){}
            void init(int _a,int _b,int _c, bool _ok){
                 a = _a; b = _b; c =_c; ok = _ok;
            }
       };
       int n,num;
       Point p[MAXN]; 
       int g[MAXN][MAXN];//凸包表面的三角形
       face F[MAXN * 8];
       double vlen(const Point &a) //向量长度
       {
              return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
       }
       double area(Point &a,Point &b,Point &c)// 三角形面积*2
       {
              return vlen((b-a) ^ (c-a));
       } 
       double volume(Point &a, Point &b , Point &c ,Point &d)//四面体体积
       {
              return ((b-a)^(c-a)) * (d-a);
       }   
       double dblcmp(Point &P,face &f)                                    
       {
              Point m=p[f.b]-p[f.a];
              Point n=p[f.c]-p[f.a];
              Point t=P-p[f.a];
              return (m ^ n) * t;
       }  
       void deal(int P,int a,int b)
       {
            int f=g[a][b];   
            if(F[f].ok)
            {
                 if(dblcmp(p[P],F[f])>eps) dfs(P,f);
                 else
                 {
                     face add(b,a,P,1); 
                     g[P][b]=g[a][P]=g[b][a]=num;
                     F[num++]=add;
                 }
            }
       }
       void dfs(int p,int now)
       {
            F[now].ok=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);
       }
       bool same(int s,int t)
       {
            Point &a=p[F[s].a];
            Point &b=p[F[s].b];
            Point &c=p[F[s].c];
            return 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;
       }