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; }
补充:软件开发 , C++ ,