三维几何

来源：互联网发布：生男生女预测表算法编辑：程序博客网时间：2024/04/28 03:39
#define eps 1e-8#define Pi acos(-1.0)#define sqr(x) ((x) * (x))int dcmp(double x){    if(fabs(x) < eps) return 0;    return x < 0 ? -1 : 1;}// 点struct point3{    double x, y, z;    point3(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}    point3 operator + (const point3 &b) const{        return point3(x + b.x, y + b.y, z + b.z);    }    point3 operator - (const point3 &b) const{        return point3(x - b.x, y - b.y, z - b.z);    }    point3 operator * (const double &k) const{        return point3(x * k, y * k, z * k);    }    point3 operator / (const double &k) const{        return point3(x / k, y / k, z / k);    }    bool operator == (const point3 &b) const{        return dcmp(x - b.x) == 0 && dcmp(y - b.y) == 0 && dcmp(z - b.z) == 0;    }    double len(){        return sqrt(x * x + y * y + z * z);    }    point3 adjust(double L){        double t = len();        L /= t;        return point3(x * L, y * L, z * L);    }};// 点积double dot(point3 a, point3 b){    return a.x * b.x + a.y * b.y + a.z * b.z;}// 叉积point3 cross(point3 a, point3 b){    return point3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y-a.y*b.x);}// 三角形abc面积的两倍double area2(point3 a, point3 b, point3 c){    return cross(b - a, c - a).len();}// 返回ab，ac，ad的混合积。它等于四面体abcd的有向体积的6倍double volume6(point3 a, point3 b, point3 c, point3 d){    return dot(d - a, cross(b - a, c - a));}//点到平面的距离(n为单位法向量)double point3_to_plane(point3 p, point3 p0, point3 n){    return fabs(dot(p - p0, n));}// 点在平面上的投影(n为单位法向量)point3 point3_plane_projection(point3 p, point3 p0, point3 n){    return p - n * dot(p - p0, n);}// p1和p2是否在线段a-b的同侧bool SameSide(point3 p1, point3 p2, point3 a, point3 b) {  return dcmp(dot(cross(b - a, p1 - a), cross(b - a, p2 - a))) >= 0;}// 点在三角形P0, P1, P2中bool point_in_tri(point3 p, point3 a, point3 b, point3 c){    return SameSide(p, a, b, c) && SameSide(p, b, a, c) && SameSide(p, c, a, b);}// 点P到直线AB的距离double point3_to_line(point3 p, point3 a, point3 b){    point3 v1 = b - a, v2 = p - a;    return cross(v1, v2).len() / v1.len();}// 直线p1 - p2到平面p0 - n的交点bool line_plane_intersection(point3 p1, point3 p2, point3 p0, point3 n, point3 &s){    point3 v = p2 - p1;    if(dcmp(dot(n, p2 - p1)) == 0) return 0;   //直线和平面平行或者直线在平面上    double t = (dot(n, p0 - p1) / dot(n, p2 - p1));    //if( t < 0 || t > 1 ) return 0;       如果是线段，判断t是不是在0 和 1之间    s = p1 + v * t;    return 1;}// 求异面直线p1＋su 和 p2+tv 的公垂线对应的s，如果平行/重合，return 0;bool line_distance(point3 p1, point3 u, point3 p2, point3 v, double &s){    double b = dot(u, v) * dot(v, u) - dot(u, u) * dot(v, v);    if(dcmp(b) == 0) return 0;    double a = dot(p2 - p1, v) * dot(v, u) - dot(p2 - p1, u) * dot(v, v);    s = a / b;    return true;}//求异面直线的距离 和 距离最近的两点double line_distance(point3 a, point3 b, point3 c, point3 d, point3 &ans1, point3 &ans2){    point3 p1 = b - a, p2 = c - d;    point3 p = cross(p1, p2);    double ret = dot(p, c - a) / p.len(), s, t;    s = (dot(cross(c - a, p2), p));    s /= p.len() * p.len();    t = (dot(cross(c - a, p1), p));    t /= p.len() * p.len();    ans1 = a + p1 * s, ans2 = c + p2 * t;    return ret;}struct face{    int a, b, c;    //表示凸包一个面上的三个点的编号    bool ok;    //表示该面是否属于最终凸包上的面};struct CH3D{    int n;    //初始顶点数    point3 P[N];    //初始顶点    int num;    //凸包表面的三角形数    face F[8*N];    //凸包表面的三角形    int g[N][N];    //凸包表面的三角形    //向量长度    double vlen(point3 a){        return a.len();    }    // 三角形面积*2    double area2(point3 a, point3 b, point3 c){        return vlen(cross(b - a, c - a));    }    //四面体有向体积*6    double volume6(point3 a, point3 b, point3 c, point3 d){        return dot(d - a, cross(b - a, c - a));    }    //正：点在面同向    double dblcmp(point3 &p, face &f){        point3 m = P[f.b] - P[f.a];        point3 n = P[f.c] - P[f.a];        point3 t = p - P[f.a];        double ret = dot(cross(m, n), t);        return ret;    }    void deal(int p, int a, int b){        int f = g[a][b];        //搜索与该边相邻的另一个平面        face add;        if(F[f].ok){            if(dblcmp(P[p], F[f]) > eps)                dfs(p, f);            else{                add.a = b;                add.b = a;                add.c = p; //这里注意顺序，要成右手系                add.ok = 1;                g[p][b] = g[a][p] = g[b][a] = num;                F[num++] = add;            }        }    }    void dfs(int p, int now){ //递归搜索所有应该从凸包内删除的面        F[now].ok = 0;        deal(p, F[now].b, F[now].a);        deal(p, F[now].c, F[now].b);        deal(p, F[now].a, F[now].c);    }    bool same(int s, int t){        point3 &a = P[F[s].a];        point3 &b = P[F[s].b];        point3 &c = P[F[s].c];        return dcmp(volume6(a, b, c, P[F[t].a])) == 0 &&               dcmp(volume6(a, b, c, P[F[t].b])) == 0 &&               dcmp(volume6(a, b, c, P[F[t].c])) == 0;    }    void create(){        face add;        num = 0;        if(n < 4) return;        //下面一段是为了保证前四个点不共面        bool flag = 1;        for(int i = 1; i < n; ++i){            if(vlen(P[0] - P[i]) > eps){                swap(P[1], P[i]);                flag = 0;                break;            }        }        if(flag) return;        flag = 1;        for(int i = 2; i < n; ++i){            if(vlen(cross(P[1] - P[0], P[i] - P[0])) > eps){                swap(P[2], P[i]);                flag = 0;                break;            }        }        if(flag) return;        flag = 1;        for(int i = 3; i < n; ++i){            if(fabs(volume6(P[0], P[1], P[2], P[i])) > eps){                swap(P[3], P[i]);                flag = 0;                break;            }        }        if(flag) return;        for(int i = 0; i < 4; ++i){            add.a = (i + 1) % 4;            add.b = (i + 2) % 4;            add.c = (i + 3) % 4;            add.ok = 1;            if(dblcmp(P[i], add) > 0)                swap(add.b, add.c);            g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;            F[num++] = add;        }        for(int i = 4; i < n; ++i){            for(int j = 0; j < num; ++j){                if(F[j].ok && dblcmp(P[i], F[j]) > eps){                    dfs(i, j);                    break;                }            }        }        int tmp = num; num = 0;        for(int i = 0; i < tmp; ++i)            if(F[i].ok)                F[num++] = F[i];    }    //凸包表面积    double area(){        double ret = 0;        if(n == 3)            return area2(P[0], P[1], P[2]) / 2.0;        for(int i = 0; i < num; ++i)            ret += area2(P[F[i].a], P[F[i].b], P[F[i].c]);        return ret / 2.0;    }    //凸包体积    double volume(){        double ret = 0;        point3 tmp(0, 0, 0);        for(int i = 0; i < num; ++i)            ret += volume6(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);        return fabs(ret) / 6.0;    }    //凸包表面三角形个数    int triangle(){        return num;    }    //凸包表面多边形个数    int polygon(){        int ret = 0, flag = 0;        for(int i = 0; i < num; ++i){            flag = 1;            for(int j = 0; j < i; ++j)                if(same(i, j)){                    flag = 0;                    break;                }            ret += flag;        }        return ret;    }    //三维凸包重心    point3 masscenter(){        point3 ans(0, 0, 0), o(0, 0, 0);        double all = 0;        for(int i = 0; i < num; ++i){            double vol = volume6(o, P[F[i].a], P[F[i].b], P[F[i].c]);            ans = ans + (o + P[F[i].a] + P[F[i].b] + P[F[i].c]) / 4 * vol;            all += vol;        }        ans = ans / all;        return ans;    }    //点到面的距离    double point3_to_face(point3 p, int i){        return fabs(volume6(P[F[i].a], P[F[i].b], P[F[i].c], p) / area2(P[F[i].a], P[F[i].b], P[F[i].c]));    }};// 把法向量为t的平面投影到x-y平面上，并旋转所有点(把平面P:a*x+b*y+c*z=d旋转平移到平面z=0)void change(point *p, point3 *pp, point3 t, int n){double a = t.x, b = t.y, c = t.z;if(dcmp(a * a + b * b) == 0){for(int i = 0;  i < n; ++i)p[i].x = pp[i].x, p[i].y = pp[i].y;return ;}for(int i = 0; i < n; ++i){double cosC = b / sqrt(a * a + b * b);double sinC = a / sqrt(a * a + b * b);point3 temp;temp.x = pp[i].x * cosC - pp[i].y * sinC;temp.y = pp[i].x * sinC + pp[i].y * cosC;temp.z = pp[i].z;pp[i] = temp;cosC = c / sqrt(a * a + b * b + c * c);sinC = sqrt(a * a + b * b) / sqrt(a * a + b * b + c * c);temp.x = pp[i].x;temp.y = pp[i].y * cosC - pp[i].z * sinC;temp.z = pp[i].y * sinC + pp[i].z * cosC;pp[i] = temp;p[i].x = pp[i].x , p[i].y = pp[i].y;}}// 任意四面体体积公式// 已知任意四面体P-ABC，记PA=a，PB=b，PC=c，cos角APB=x，cos角BPC=y，cos角CPA=z，// 则：V=1/6*abc*sqrt(1+2xyz-x^2-y^2-z^2)
0 0