POJ 3528 三维凸包模板

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const double eps = 1e-8;
const int maxn = 505;
struct Point {
	double x, y, z;
	Point(double x=0, double y=0, double z=0):
		x(x), y(y), z(z){}
	Point operator+(const Point &t) const {
		return Point(x+t.x, y+t.y, z+t.z);
	}
	Point operator-(const Point &t) const {
		return Point(x-t.x, y-t.y, z-t.z);
	}
	Point operator*(const Point &t) const {
		return Point(y*t.z-z*t.y, z*t.x-x*t.z, x*t.y-y*t.x);
	}
	double operator^(const Point &t) const {
		return  x*t.x+y*t.y+z*t.z;
	}
	Point operator*(const double &t) const {
		return Point(x*t, y*t, z*t);
	}
	Point operator/(const double &t) const {
		return Point(x/t, y/t, z/t);
	}
	void in() {
		scanf("%lf%lf%lf", &x, &y, &z);
	}
};
double vlen(Point t) {
	return sqrt(t.x*t.x+t.y*t.y+t.z*t.z);
}
struct CH3D {
	struct face {
		int a, b, c;
		bool ok;
		face(int a, int b, int c):
			a(a), b(b), c(c), ok(1){}
		face(){}
	}f[8*maxn];
	int cnt, n, to[maxn][maxn];
	Point p[maxn];
	bool fix() { //为了保证前四个点不公面, 若有保证就可以不写fix函数
		int i;
		bool ok = true;
		for(i = 1; i < n; i++) //使前两点不公点
			if(vlen(p[i]-p[0]) > eps) {
				swap(p[i], p[1]);
				ok = false;
				break;
			}
		if(ok) return 0; ok = true;
		for(i = 2; i < n; i++)  //使前三点不公线
			if(vlen( (p[1]-p[0])*(p[i]-p[0]) ) > eps) {
				swap(p[i], p[2]);
				ok = false;
				break;
			}
		if(ok) return 0; ok = true;
		for(i = 3; i < n; i++)  //使前四点不共面
			if(fabs( (p[1]-p[0])*(p[2]-p[0])^(p[i]-p[0]) ) > eps) {
				swap(p[i], p[3]);
				ok = false;
				break;
			}
		if(ok) return 0;
		return 1;
	}
	double ptoface(Point &t, face &f) { //判点是否在面的同侧，>0 同侧
		return volume(p[f.a], p[f.b], p[f.c], t);
	}
	void dfs(int i, int j) {
		f[j].ok = 0;
		deal(i, f[j].b, f[j].a);
		deal(i, f[j].c, f[j].b);
		deal(i, f[j].a, f[j].c);
	}
	void deal(int i, int a, int b) {
		int j = to[a][b];
		if(f[j].ok) {
			if(ptoface(p[i], f[j]) > eps)
				dfs(i, j);
			else {
				face add(b, a, i);
				to[b][a] = to[a][i] = to[i][b] = cnt;
				f[cnt++] = add;
			}
		}
	}
	void creat() { //构造凸包
		if(n < 4 || !fix()) return;
		int i, j;
		cnt = 0;
		for(i = 0; i < 4; i++) {
			face add((i+1)%4, (i+2)%4, (i+3)%4);
			if(ptoface(p[i], add) > eps)
				swap(add.b, add.c);
			to[add.a][add.b] = to[add.b][add.c] = to[add.c][add.a] = cnt;
			f[cnt++] = add;
		}
		for(i = 4; i < n; i++) {
			for(j = 0; j < cnt; j++)
				if(f[j].ok && ptoface(p[i], f[j]) > eps) {
					dfs(i, j);
					break;
				}
		}
		int t = cnt; cnt = 0;
		for(i = 0; i < t; i++)
			if(f[i].ok) f[cnt++] = f[i];
	}
	double area() { //凸包表面积
		double ret = 0.0;
		for(int i = 0; i < cnt; i++)
			ret += area(p[f[i].a], p[f[i].b], p[f[i].c]);
		return ret*0.5;
	}
	double volume() { //凸包体积
		Point o;  //o是原点
		double ret = 0.0;
		for(int i = 0; i < cnt; i++)
			ret += volume(o, p[f[i].a], p[f[i].b], p[f[i].c]);
		return fabs(ret/6.0);
	}
	bool same(int s, int t) { // 判两个平面是否为同一平面
        Point &a = p[f[s].a], &b = p[f[s].b], &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;
    }
	int faceCnt() {  //凸包的面数（除去相同的平面）
		int i, j;
	    int ans = 0;
	    for(i = 0; i < cnt; i++) {
	        bool ok = 1;
	        for(j = 0; j < i; j++) {
	            if(same(i, j)) {
	                ok = 0;
	                break;
	            }
	        }
	            ans += ok;
	    }
	    return ans;
	}
	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);
	}
}hull;
int main() {
	int i;
	while(~scanf("%d", &hull.n)) {
		for(i = 0; i < hull.n; i++) hull.p[i].in();
		hull.creat();
		printf("%.3f\n", hull.area());
	}
	return 0;
}