三维凸包模板_HDU 3662 3D Convex Hull

http://acm.hdu.edu.cn/showproblem.php?pid=3662

typedef long long LL;
typedef unsigned long long ULL;
typedef vector <int> VI;
const int INF = 0x3f3f3f3f;
const double eps = 1e-10;
const int MOD = 100000007;
const int MAXN = 400;
const double PI = acos(-1.0);

struct Point
{
    double x, y, z;
    Point(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}

    inline void read()
    {
        scanf("%lf%lf%lf", &x, &y, &z);
    }

    //两向量之差
    inline Point operator- (Point p)
    {
        return Point(x - p.x, y - p.y, z - p.z);
    }

    //两向量之和
    inline Point operator+ (Point p)
    {
        return Point(x + p.x, y + p.y, z + p.z);
    }

    //叉乘
    inline Point operator* (Point p)
    {
        return Point(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);
    }

    inline Point operator* (double d)
    {
        return Point(x * d, y * d, z * d);
    }

    inline Point operator/ (double d)
    {
        return Point(x / d, y / d, z / d);
    }

    //点乘
    inline double operator^ (Point p)
    {
        return (x * p.x + y * p.y + z * p.z);
    }
};

struct CH3D
{
    struct face
    {
        //表示凸包一个面上的三个点的编号
        int a,b,c;
        //表示该面是否属于最终凸包上的面
        bool ok;
    };
    //初始顶点数
    int n;
    //初始顶点
    Point P[MAXN];
    //凸包表面的三角形数
    int num;
    //凸包表面的三角形
    face F[8*MAXN];
    //凸包表面的三角形
    int g[MAXN][MAXN];
    //向量长度
    inline double Length(Point a)
    {
        return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
    }
    //叉乘
    inline Point cross(Point a, Point b, Point c)
    {
        return Point((b.y - a.y) * (c.z - a.z) - (b.z - a.z) * (c.y - a.y),
                     (b.z - a.z) * (c.x - a.x) - (b.x - a.x) * (c.z - a.z),
                     (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)
                    );
    }
    //三角形面积*2
    inline double area(Point a, Point b, Point c)
    {
        return Length((b - a) * (c - a));
    }
    //四面体有向体积*6
    inline double volume(Point a, Point b, Point c, Point d)
    {
        return (b - a) * (c - a) ^ (d - a);
    }
    //正：点在面同向
    inline 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];//搜索与该边相邻的另一个平面
        face add;
        if(F[f].ok)
        {
            if(dblcmp(P[p],F[f])>eps)
                dfs(p,f);
            else
            {
                add.a = b;
                add.b = a;
                add.c = p;//这里注意顺序，要成右手系
                add.ok = true;
                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;
    }
    //构建三维凸包
    void create()
    {
        face add;

        num = 0;
        if(n < 4) return;
        //**********************************************
        //此段是为了保证前四个点不共面
        bool flag = true;
        FF(i, 1, n)
        {
            if(Length(P[0] - P[i]) > eps)
            {
                swap(P[1], P[i]);
                flag=false;
                break;
            }
        }
        if(flag) return;
        flag = true;
        //使前三个点不共线
        FF(i, 2, n)
        {
            if(Length((P[0] - P[1]) * (P[1] - P[i])) > eps)
            {
                swap(P[2], P[i]);
                flag = false;
                break;
            }
        }
        if(flag) return;
        flag = true;
        //使前四个点不共面
        FF(i, 3, n)
        {
            if(fabs((P[0] - P[1]) * (P[1] - P[2]) ^ (P[0] - P[i])) > eps)
            {
                swap(P[3], P[i]);
                flag = false;
                break;
            }
        }
        if(flag) return;
        //*****************************************
        REP(i, 4)
        {
            add.a = (i + 1) % 4;
            add.b = (i + 2) % 4;
            add.c = (i + 3) % 4;
            add.ok = true;
            if(dblcmp(P[i], add) > 0)
                swap(add.b, add.c);
            g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
            F[num++] = add;
        }
        FF(i, 4, n)
        {
            REP(j, num)
            {
                if(F[j].ok && dblcmp(P[i],F[j]) > eps)
                {
                    dfs(i, j);
                    break;
                }
            }
        }
        int tmp = num;
        num = 0;
        REP(i, tmp)
            if(F[i].ok)
                F[num++] = F[i];

    }
    //表面积
    double area()
    {
        double res = 0;
        if(n == 3)
        {
            Point p = cross(P[0], P[1], P[2]);
            res = Length(p) / 2.0;
            return res;
        }
        REP(i, num)
            res += area(P[F[i].a], P[F[i].b], P[F[i].c]);
        return res / 2.0;
    }
    double volume()
    {
        double res = 0;
        Point tmp(0, 0, 0);
        REP(i, num)
            res += volume(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);
        return fabs(res / 6.0);
    }
    //表面三角形个数
    inline int triangle()
    {
        return num;
    }
    //表面多边形个数
    int polygon()
    {
        int res = 0, flag;
        REP(i, num)
        {
            flag = 1;
            REP(j, i)
                if(same(i, j))
                {
                    flag = 0;
                    break;
                }
            res += flag;
        }
        return res;
    }
    //三维凸包重心
    Point barycenter()
    {
        Point ans(0,0,0), t(0,0,0);
        double all = 0, vol;
        REP(i, num)
        {
            vol = volume(t, P[F[i].a], P[F[i].b], P[F[i].c]);
            ans = ans + (t + P[F[i].a] + P[F[i].b] + P[F[i].c]) / 4.0 * vol;
            all += vol;
        }
        return ans / all;
    }
    //点到面的距离
    inline double ptoface(Point p, int i)
    {
        double Len = Length((P[F[i].b] - P[F[i].a]) * (P[F[i].c] - P[F[i].a]));
        return fabs(volume(P[F[i].a], P[F[i].b], P[F[i].c],p) / Len);
    }
} hull;

int main()
{
    // freopen("in.txt","r",stdin);
    // freopen("out.txt","w",stdout);
    while(~RI(hull.n))
    {
        REP(i, hull.n) hull.P[i].read();
        hull.create();
        printf("%d\n", hull.polygon());
    }
    return 0;
}