三维凸包 CH3D

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;

const int N = 505;  // 点数量
const double eps = 1e-8;

struct point {
    double x, y, z;
    point() {}
    point(double xx, double yy, double zz) : x(xx), y(yy), z(zz) {}

    point operator+(const point p1) {  //两向量之和
        return point(x + p1.x, y + p1.y, z + p1.z);
    }  // 求重心的时候用到

    point operator-(const point p1) {  //两向量之差
        return point(x - p1.x, y - p1.y, z - p1.z);
    }

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

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

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

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

struct CH3D {
    struct face {
        int a, b, c;  //表示凸包一个面上三个点的编号
        bool ok;      //表示该面是否属于最终凸包中的面
    };

    int n;       //初始顶点数
    point P[N];  //初始顶点

    int num;        //凸包表面的三角形数
    face F[8 * N];  //凸包表面的三角形
    int g[N][N];    //凸包表面的三角形

    double vlen(point a) {  //向量长度
        return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
    }

    point cross(const point &a, const point &b, const point &c) {  //叉乘
        return point((b.y - a.y) * (c.z - a.z) - (b.z - a.z) * (c.y - a.y),
                     -((b.x - a.x) * (c.z - a.z) - (b.z - a.z) * (c.x - a.x)),
                     (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x));
    }  // 三点确定一个平面 叉乘返回该平面的法向量

    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) {  //四面体有向体积*6
        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];  //搜索与该边相邻的另一个平面
        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 = 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;
    }

    void create() {  //构建三维凸包
        face add;
        bool flag = true;
        num = 0;
        if (n < 4) return;
        //此段是为了保证前四个点不共面,若以保证,则可去掉
        /*********************************/
        for (int i = 1; i < n; i++) {
            if (vlen(P[0] - P[i]) > eps) {
                swap(P[1], P[i]);
                flag = false;
                break;
            }
        }
        if (flag) return;
        flag = true;
        for (int i = 2; i < n; i++) {  //使前三点不共线
            if (vlen((P[0] - P[1]) * (P[1] - P[i])) > eps) {
                swap(P[2], P[i]);
                flag = false;
                break;
            }
        }
        if (flag) return;
        flag = true;
        for (int i = 3; i < n; i++) {  //使前四点不共面
            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;
        /*********************************/
        for (int i = 0; i < 4; i++) {
            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;
        }
        for (int i = 4; i < n; i++) {
            for (int j = 0; j < num; j++) {
                if (F[j].ok && dblcmp(P[i], F[j]) > eps) {
                    dfs(i, j);
                    break;
                }
            }
        }
        int tmp = num;
        num = 0;
        for (int i = 0; i < tmp; i++)
            if (F[i].ok) F[num++] = F[i];
    }

    double area() {  //表面积
        double res = 0.0;
        if (n == 3) {
            point p = cross(P[0], P[1], P[2]);
            res = vlen(p) / 2.0;
            return res;
        }
        for (int i = 0; i < num; i++)
            res += area(P[F[i].a], P[F[i].b], P[F[i].c]);
        return res / 2.0;
    }

    double volume() {  //体积
        double res = 0.0;
        point tmp(0, 0, 0);
        for (int i = 0; i < num; i++)
            res += volume(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);
        return fabs(res / 6.0);
    }

    int triangle() {  //表面三角形个数
        return num;
    }

    int polygon() {  //表面多边形个数
        int res = 0, flag;
        for (int i = 0; i < num; i++) {
            flag = 1;
            for (int j = 0; j < i; j++)
                if (same(i, j)) {
                    flag = 0;
                    break;
                }
            res += flag;
        }
        return res;
    }

    point barycenter() {
        point ans(0, 0, 0), o(0, 0, 0);
        double all = 0;
        for (int i = 0; i < num; i++) {
            double vol = volume(o, P[F[i].a], P[F[i].b], P[F[i].c]);
            ans = ans + (o + P[F[i].a] + P[F[i].b] + P[F[i].c]) / 4.0 * vol;
            all += vol;
        }
        ans = ans / all;
        return ans;
    }
    //点到面的距离
    double ptoface(point p, int i) {
        return fabs(volume(P[F[i].a], P[F[i].b], P[F[i].c], p) /
                    vlen((P[F[i].b] - P[F[i].a]) * (P[F[i].c] - P[F[i].a])));
    }
} hull;

int main() {
    while (scanf("%d", &hull.n) != EOF) {
        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.create();
        printf("%.3f\n", hull.area());
    }
    return 0;
}