POJ 3528 三维凸包模板

2014-11-24 00:59:26 · 作者: · 浏览: 3
#include 
#include 
#include 
#include 
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; }