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) { //判点是否在面的同侧，> 同侧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;}