检测立方体和圆锥体是否相互交叉？

Question

考虑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)？

注意:此处交点定义为:它存在一个位于立方体和圆锥内的非零三维体积(如果立方体位于圆锥内,或者如果圆锥位于立方体内,则它们相交).

Answer 1

1. 想象2无限的线条

   • 锥形轴
   • 穿过P与锥轴垂直的点(起点的立方体中心)的直线.

   圆锥轴是众所周知的,因此很容易,第二行定义为

   P+t*(perpendicular vector to cone axis)
   

   该矢量可以通过锥轴矢量和垂直于图像的矢量(假设Z轴)的叉积获得.该t是标量值参数...

2. 计算这两条线/轴的交点

   如果你不知道方程式派生他们或google他们.让交点成为Q

3. 如果交点Q不在锥内

   (在顶点和底点之间)然后点P不是相交锥.从交叉方程式中,您将获得参数t1和t2

   • 让我们t1为P轴线
   • 和t2锥轴线

   如果你的轴线方向矢量也是锥长,那么交点就在锥内 t2 = <0,1>

4. 如果P不是三角形内部(由这两个轴产生的切割锥形到平面)

   这也很容易,你知道的位置Q内锥(t2),所以你知道,锥体在P轴从Q到距离R*t2这里R是圆锥体底座半径.因此,您可以计算|P-Q|并检查它是否<=R*t2直接使用t1(如果P轴方向矢量是单位).

   如果距离较大则R*t2点P不与锥相交.

5. 如果#3和#4为正,则P与锥相交

   

   • 希望你不介意这里是你的形象,为了清晰起见,添加了一些东西

[笔记]

现在,当没有立方体的顶点与锥体相交但是立方体本身与锥体相交时,存在边缘情况.当||P-Q|-R*t2| = <0,half cube size>在这种情况下你应该检查更多点然后沿最近的立方体面立方体顶点时,可能会发生这种情况.

另一种方法是:

1. 为锥形创建变换矩阵

   哪里:

   • 它的顶点是原点
   • 它的轴+Z为轴
   • 和XY平面平行于它的底座

   所以任何一点都在锥内

   • Z = <0,h>
   • X*X + Y*Y <= (R*Z/h)^2 要么 X*X + Y*Y <= (R*Z*tan(angle))^2

2. 将立方体顶点转换为锥体空间

   并且检查是否有任何顶点在锥内也可以检查所有立方体边缘线,条件来自#1(代数)或沿着立方体面使用更多点,如前面的方法.

聊天讨论:http://chat.stackoverflow.com/rooms/48756/discussion-between-spektre-and-joojaa

Answer 2

对于代码

这个答案会比你的问题略胜一筹(例如我认为是一个方框而不是一个方块).适应你的情况应该非常简单.

一些定义

/*
    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一个留给读者练习.假设该线在锥形空间中.分段方向未规范化; 相反,该段具有方向向量的长度,并从该点开始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个向量和一个点.矢量未标准化,以简化计算.

/*
    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

方向矢量是(a,b,c),点(u,v,w)位于该线上.

现在,让我们把方程式放在一起:

(u + at)² + (v + bt)² = (H - w - ct)² * tan²(alpha)

然后在开发和重新分解这个等式之后,我们得到以下结果:

At² + Bt + C = 0

其中A,B和C显示在第一个交叉函数中.只需解决此问题并检查z和t上的边界条件.