计算几何_三维凸包(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++;}}}}};?