考虑3D中的两个几何对象:

  • 与轴对齐并由其中心位置及其范围(边长)定义的立方体
  • 一个圆锥体,它与轴不对齐,并且由其顶点位置,其底部中心的位置以及顶点处的半角定义

  • 以下是在C++中定义这些对象的小代码:
    // Preprocessor
    #include <iostream>
    #include <cmath>
    #include <array>
    
    // 3D cube from the position of its center and the side extent
    class cube
    {
        public:
            cube(const std::array<double, 3>& pos, const double ext)
            : _position(pos), _extent(ext)
            {;}
            double center(const unsigned int idim)
                {return _position[idim];}
            double min(const unsigned int idim)
                {return _position[idim]-_extent/2;}
            double max(const unsigned int idim)
                {return _position[idim]+_extent/2;}
            double extent()
                {return _extent;}
            double volume()
                {return std::pow(_extent, 3);}
        protected:
            std::array<double, 3> _position;
            double _extent;
    };
    
    // 3d cone from the position of its vertex, the base center, and the angle
    class cone
    {
        public:
            cone(const std::array<double, 3>& vert,
                 const std::array<double, 3>& bas,
                 const double ang)
            : _vertex(vert), _base(bas), _angle(ang)
            {;}
            double vertex(const unsigned int idim)
                {return _vertex[idim];}
            double base(const unsigned int idim)
                {return _base[idim];}
            double angle()
                {return _angle;}
            double height()
                {return std::sqrt(std::pow(_vertex[0]-_base[0], 2)+std::pow(
                _vertex[1]-_base[1], 2)+std::pow(_vertex[2]-_base[2], 2));}
            double radius()
                {return std::tan(_angle)*height();}
            double circle()
                {return 4*std::atan(1)*std::pow(radius(), 2);}
            double volume()
                {return circle()*height()/3;}
        protected:
            std::array<double, 3> _vertex;
            std::array<double, 3> _base;
            double _angle;
    };
    

    我想编写一个函数来检测立方体和圆锥体的交点是否为空:
    // Detect whether the intersection between a 3d cube and a 3d cone is not null
    bool intersection(const cube& x, const cone& y)
    {
        // Function that returns false if the intersection of x and y is empty
        // and true otherwise
    }
    

    这是问题的图解(图解是2D,但我的问题是3D):

    如何有效地做到这一点(我正在寻找一种算法,因此答案可以是C,C++或Python)?

    注意:这里的交集定义为:它存在于立方体和圆锥体中的非零3D体积(如果立方体在圆锥体内部,或者如果圆锥体在立方体内部,则它们相交)。

    最佳答案

    对于代码

    该答案将比您的问题更笼统(例如,我考虑使用盒子而不是立方体)。适应您的情况应该非常简单。

    一些定义

    /*
        Here is the cone in cone space:
    
                +        ^
               /|\       |
              /*| \      | H
             /  |  \     |
            /       \    |
           +---------+   v
    
        * = alpha (angle from edge to axis)
    */
    struct Cone // In cone space (important)
    {
        double H;
        double alpha;
    };
    
    /*
        A 3d plane
          v
          ^----------+
          |          |
          |          |
          +----------> u
          P
    */
    struct Plane
    {
        double u;
        double v;
        Vector3D P;
    };
    
    // Now, a box.
    // It is assumed that the values are coherent (that's only for this answer).
    // On each plane, the coordinates are between 0 and 1 to be inside the face.
    struct Box
    {
        Plane faces[6];
    };
    

    线-锥相交

    现在,让我们计算线段和圆锥之间的交点。请注意,我将在圆锥空间中进行计算。另请注意,我将Z轴设为垂直轴。将其更改为Y,留给读者练习。假定该线位于圆锥空间中。段方向未标准化;相反,该段的长度是方向 vector 的长度,从P点开始:
    /*
        The segment is points M where PM = P + t * dir, and 0 <= t <= 1
        For the cone, we have 0 <= Z <= cone.H
    */
    bool intersect(Cone cone, Vector3D dir, Vector3D P)
    {
        // Beware, indigest formulaes !
        double sqTA = tan(cone.alpha) * tan(cone.alpha);
        double A = dir.X * dir.X + dir.Y * dir.Y - dir.Z * dir.Z * sqTA;
        double B = 2 * P.X * dir.X +2 * P.Y * dir.Y - 2 * (cone.H - P.Z) * dir.Z * sqTA;
        double C = P.X * P.X + P.Y * P.Y - (cone.H - P.Z) * (cone.H - P.Z) * sqTA;
    
        // Now, we solve the polynom At² + Bt + C = 0
        double delta = B * B - 4 * A * C;
        if(delta < 0)
            return false; // No intersection between the cone and the line
        else if(A != 0)
        {
            // Check the two solutions (there might be only one, but that does not change a lot of things)
            double t1 = (-B + sqrt(delta)) / (2 * A);
            double z1 = P.Z + t1 * dir.Z;
            bool t1_intersect = (t1 >= 0 && t1 <= 1 && z1 >= 0 && z1 <= cone.H);
    
            double t2 = (-B - sqrt(delta)) / (2 * A);
            double z2 = P.Z + t2 * dir.Z;
            bool t2_intersect = (t2 >= 0 && t2 <= 1 && z2 >= 0 && z2 <= cone.H);
    
            return t1_intersect || t2_intersect;
        }
        else if(B != 0)
        {
            double t = -C / B;
            double z = P.Z + t * dir.Z;
            return t >= 0 && t <= 1 && z >= 0 && z <= cone.H;
        }
        else return C == 0;
    }
    

    矩形-圆锥形相交

    现在,我们可以检查平面图的矩形部分是否与圆锥相交(这将用于检查立方体的面是否与圆锥相交)。仍处于圆锥空间中。该计划以对我们有帮助的方式通过:2个 vector 和一个点。不对 vector 进行归一化,以简化计算。
    /*
        A point M in the plan 'rect' is defined by:
            M = rect.P + a * rect.u + b * rect.v, where (a, b) are in [0;1]²
    */
    bool intersect(Cone cone, Plane rect)
    {
        bool intersection = intersect(cone, rect.u, rect.P)
                         || intersect(cone, rect.u, rect.P + rect.v)
                         || intersect(cone, rect.v, rect.P)
                         || intersect(cone, rect.v, rect.P + rect.u);
    
        if(!intersection)
        {
            // It is possible that either the part of the plan lie
            // entirely in the cone, or the inverse. We need to check.
            Vector3D center = P + (u + v) / 2;
    
            // Is the face inside the cone (<=> center is inside the cone) ?
            if(center.Z >= 0 && center.Z <= cone.H)
            {
                double r = (H - center.Z) * tan(cone.alpha);
                if(center.X * center.X + center.Y * center.Y <= r)
                    intersection = true;
            }
    
            // Is the cone inside the face (this one is more tricky) ?
            // It can be resolved by finding whether the axis of the cone crosses the face.
            // First, find the plane coefficient (descartes equation)
            Vector3D n = rect.u.crossProduct(rect.v);
            double d = -(rect.P.X * n.X + rect.P.Y * n.Y + rect.P.Z * n.Z);
    
            // Now, being in the face (ie, coordinates in (u, v) are between 0 and 1)
            // can be verified through scalar product
            if(n.Z != 0)
            {
                Vector3D M(0, 0, -d/n.Z);
                Vector3D MP = M - rect.P;
                if(MP.scalar(rect.u) >= 0
                   || MP.scalar(rect.u) <= 1
                   || MP.scalar(rect.v) >= 0
                   || MP.scalar(rect.v) <= 1)
                    intersection = true;
            }
        }
        return intersection;
    }
    

    盒子-圆锥形相交

    现在,最后一部分:整个立方体:
    bool intersect(Cone cone, Box box)
    {
        return intersect(cone, box.faces[0])
            || intersect(cone, box.faces[1])
            || intersect(cone, box.faces[2])
            || intersect(cone, box.faces[3])
            || intersect(cone, box.faces[4])
            || intersect(cone, box.faces[5]);
    }
    

    对于数学

    仍然在圆锥空间中,圆锥方程为:
    // 0 is the base, the vertex is at z = H
    x² + y² = (H - z)² * tan²(alpha)
    0 <= z <= H
    

    现在,3D中线的参数方程为:
    x = u + at
    y = v + bt
    z = w + ct
    

    方向 vector 为(a,b,c),点(u,v,w)位于线上。

    现在,让我们将这些方程放在一起:
    (u + at)² + (v + bt)² = (H - w - ct)² * tan²(alpha)
    

    然后,在开发并重新分解了该方程后,我们得到以下结果:
    At² + Bt + C = 0
    

    其中A,B和C在第一个交集函数中显示。只需解决这个问题,然后检查z和t上的边界条件。

    10-07 20:22