POJ 3528--Ultimate Weapon(三维凸包)

#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
using namespace std;
const int N = 510;
const double eps = 1e-8;
typedef struct point3 {
    double x, y, z;
    point3() {  }
    point3(double a, double b, double c) :x(a), y(b), z(c) {  }
    point3 operator -(const point3 &b)const {       //返回减去后的新点
        return point3(x - b.x, y - b.y, z - b.z);
    }
    point3 operator +(const point3 &b)const {       //返回加上后的新点
        return point3(x + b.x, y + b.y, z + b.z);
    }
    //数乘计算
    point3 operator *(const double &k)const {   //返回相乘后的新点
        return point3(x * k, y * k, z*k);
    }
    point3 operator /(const double &k)const {   //返回相除后的新点
        return point3(x / k, y / k, z / k);
    }
    double operator *(const point3 &b)const {   //点乘
        return (x*b.x + y*b.y + z*b.z);
    }
    point3 operator ^(const point3 &p) const {  //叉积
        return point3(y*p.z - p.y*z, z*p.x - x*p.z, x*p.y - y*p.x);
    }
    double vlen()const {    //向量的模
        return sqrt(x*x + y*y + z*z);
    }
}point3;
struct fac {
    int a, b, c;//凸包一个面上的三个点的编号
    bool ok;    //该面是否是最终凸包中的面
};
struct T3dhull {
    int n;                  //初始点数
    point3 ply[N];          //初始点
    int trianglecnt;        //凸包上三角形数
    fac tri[N];             //凸包三角形可证明被创建的面不会超过6N
    int vis[N][N];          //点i到点j是属于哪个面
    double dist(point3 a) { return sqrt(a.x*a.x + a.y*a.y + a.z*a.z); }         //两点长度
    double area(point3 a, point3 b, point3 c){return dist((b - a) ^ (c - a));}  //三角形面积*2

    //返回四面体有向体积*6
    //在储存面时，保证面的法线方向朝向凸包外部，如果在某一平面和点p所组成的四面体的有向体积为正，则p点在凸包外部，并且此点可以被p点看见。
    double volume(point3 a, point3 b, point3 c, point3 d){
        return ((b - a) ^ (c - a))* (d - a);
    }
    double ptoplane(point3 &p, fac &f){     //点到平面距离，体积法
        point3 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])) > eps)
                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) {
        point3 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])) < eps
            && fabs(volume(a, b, c, ply[tri[e].b])) < eps
            && fabs(volume(a, b, c, ply[tri[e].c])) < eps;
    }
    void construct()//构造凸包
    {
        int i, j;
        trianglecnt = 0;
        if (n<4) return;
        bool tmp = true;
        for (i = 1; i < n; i++) //前两点不共点
        {
            if ((dist(ply[0] - ply[i])) > eps)
            {
                swap(ply[1], ply[i]);
                tmp = false;
                break;
            }
        }
        if (tmp)return;
        tmp = true;
        for (i = 2; i < n; i++)     //前三点不共线
        {
            if ((dist((ply[0] - ply[1]) ^ (ply[1] - ply[i]))) > eps)
            {
                swap(ply[2], ply[i]);
                tmp = false;
                break;
            }
        }
        if (tmp) return;
        tmp = true;
        for (i = 3; i < n; i++) //前四点不共面
        {
            if (fabs(((ply[0] - ply[1]) ^ (ply[1] - ply[2]))* (ply[0] - ply[i]))>eps)
            {
                swap(ply[3], ply[i]);
                tmp = false;
                break;
            }
        }
        if (tmp)return;
        fac 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;
            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 (i = 4; i < n; i++)     //构建更新凸包
        {
            for (j = 0; j < trianglecnt; j++)
            {
                if (tri[j].ok && (ptoplane(ply[i], tri[j])) > eps)
                {
                    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()
    {
        point3 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);
    }
}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("%.3lf\n", hull.area());
    }
    return 0;
}