计算几何_三维凸包(3d convex hull)

const double eps = 1e-8;
typedef list<int>::iterator liit;
inline int sign(double d){
	if(d < -eps) return -1;
	return (d > eps) ? 1 : 0;
}
struct point{
	double x, y, z;
	point(double _x=0, double _y=0, double _z=0):x(_x), y(_y), z(_z) {}
	void read(){ scanf("%lf%lf%lf", &x, &y, &z); }
	bool operator==(const point& tp){
		return (!sign(x-tp.x)) && (!sign(y-tp.y)) && (!sign(z-tp.z));
	}
	bool isOrg(){
		return !sign(x) && !sign(y) && !sign(z);
	}
	point operator-(point tp){
		return point(x-tp.x, y-tp.y, z-tp.z);
	}
};
struct face{
	list<int> pos;
	point v;  //法向量
	double d; //常量
};
inline bool cmp(point p1, point p2){
	return (p1.x < p2.x) || (p1.x == p2.x && p1.y < p2.y)
		|| (p1.x == p2.x && p1.y == p2.y && p1.z < p2.z);
}
//向量st-->ed1和st-->ed2的叉积
inline point xmul3d(point st, point ed1, point ed2){
	ed1 = ed1 - st;
	ed2 = ed2 - st;
	return point(ed1.y*ed2.z-ed2.y*ed1.z, ed1.z*ed2.x-ed2.z*ed1.x,
		ed1.x*ed2.y-ed2.x*ed1.y);
}
//向量(0,0)-->p1和(0,0)-->ed2的叉积
inline double dmul3d(point p1, point p2){
	return p1.x*p2.x+p1.y*p2.y+p1.z*p2.z;
}
inline double dist3d(point p1, point p2){
	return sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y)
		+ (p1.z-p2.z)*(p1.z-p2.z));
}
//三维凸包类
struct convex3d{
	static const int N = 305;  //点的数目的最大值
	point ps[N];       //求解凸包的点集
	list<face> convex; //convex存储三维的凸包的各个面，这些面是三维空间的凸多边形
	int pn;            //点的个数
	int fsign[N][N];
	//添加由ps里的第a,b,c个点组成的面
	void addFace(int a, int b, int c){
		face tf;
		tf.pos.push_back(a);
		tf.pos.push_back(b);
		tf.pos.push_back(c);
		tf.v = xmul3d(ps[a], ps[b], ps[c]);
		tf.d = -dmul3d(ps[a], tf.v);
		convex.push_back(tf);
	}
	//插入第i个点是，对面f进行处理
	int handleFace(face f, int i){
		int s = sign(dmul3d(f.v, ps[i]) + f.d);
		liit now, nxt;
		now = f.pos.begin();
		nxt = f.pos.begin();
		nxt++;
		for(; nxt != f.pos.end(); nxt++, now++){
			fsign[*now][*nxt] = s;
		}
		fsign[*now][*f.pos.begin()] = s;
		return s;
	}
	//判断第i个点是否在f里面
	bool inFace(face f, int i){
		liit now, nxt;
		int pn, nn, s;
		now = f.pos.begin();
		nxt = f.pos.begin();
		nxt++;
		for(pn = nn = 0; nxt != f.pos.end(); nxt++, now++){
			s = sign(dmul3d(xmul3d(ps[*now], ps[*nxt], ps[i]), f.v));
			if(s == 1) pn++;
			else if(s == -1) nn++;
		}
		s = sign(dmul3d(xmul3d(ps[*now], ps[*f.pos.begin()], ps[i]), f.v));
		if(s == 1) pn++;
		else if(s == -1) nn++;
		if(pn >= 1 && nn >= 1) return false;
		return true;
	}
	//扩展面f,返回true表示需要删除当前的面
	bool extFace(face& f, int i){
		liit now, nxt;
		bool flag = false;
		now = f.pos.begin();
		nxt = f.pos.begin();
		nxt++;
		if(fsign[*now][*nxt] == 0){
			list<int> tpos;
			while(true){
				if(sign(dmul3d(xmul3d(ps[i], ps[*now], ps[*nxt]), f.v)) >= 1){
					break;
				}
				now++;
				if(now == f.pos.end()) now = f.pos.begin();
				nxt++;
				if(nxt == f.pos.end()) nxt = f.pos.begin();
			}
			tpos.push_back(*now);
			int st = *now;
			while(*nxt != st){
				if(sign(dmul3d(xmul3d(ps[*now], ps[i], ps[*nxt]), f.v)) >= 0){
					break;
				}
				tpos.push_back(*nxt);
				now++;
				if(now == f.pos.end()) now = f.pos.begin();
				nxt++;
				if(nxt == f.pos.end()) nxt = f.pos.begin();
			}
			tpos.push_back(i);
			while(true){
				now++;
				if(now == f.pos.end()) now = f.pos.begin();
				nxt++;
				if(nxt == f.pos.end()) nxt = f.pos.begin();
				if(*now == st) break;
				if(sign(dmul3d(xmul3d(ps[i], ps[*now], ps[*nxt]), f.v)) >= 1){
					tpos.push_back(*now);
				}
			}
			f.pos = tpos;

		}else if(fsign[*now][*nxt] > 0){
			for(; nxt != f.pos.end(); now++, nxt++){
				if(fsign[*nxt][*now] < 0){
					addFace(*now, *nxt, i);
				}
			}
			if(fsign[*f.pos.begin()][*now] < 0){
				addFace(*now, *f.pos.begin(), i);
			}
			flag = true;
		}
		return flag;
	}
	//对ps里的pn个点求最小包围多面体，结果放在convex中
	void initConvex(){
		sort(ps, ps+pn, cmp);
		pn = unique(ps, ps+pn) - ps;
		convex.clear();
		if(pn <= 2){
			return;
		}
		int a, b, c;
		double ab, bc, ac;
		a = 0;
		b = 1;
		for(c = 2; c < pn; c++){
			ab = dist3d(ps[a], ps[b]);
			bc = dist3d(ps[b], ps[c]);
			ac = dist3d(ps[a], ps[c]);
			if(sign(ab+bc-ac) == 0){
				b = c;
			}else if(sign(ab+ac-bc) == 0){
				a = c;
			}else if(sign(ac+bc-ab) != 0){
				break;
			}
		}
		if(c == pn){
			return;
		}
		int i, size, j;
		list<face>::iterator it;
		addFace(a, b, c);
		addFace(a, c, b);
		for(i = c+1; i < pn; i++){
			size = convex.size();
			for(it = convex.begin(), j = 0; j < size; j++, it++){
				if(handleFace(*it, i) == 0 && inFace(*it, i)){
						break;
				}
			}
			if(j < size) continue;
			for(it = convex.begin(), j = 0; j < size; j++){
				if(extFace(*it, i)){
					it = convex.erase(it);
				}else{
					it++;
				}
			}
		}
	}
};