Detect if a cube and a cone intersect each other?

Question

Consider two geometrical objects in 3D:

• a cube aligned with the axes and defined by the position of its center and its extent (edge length)
• a cone not aligned with the axes and defined by the position of its vertex, the position of the center of its base, and the half-angle at the vertex

Here is a small code to define these objects in 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;
};

I would like to write a function to detect whether the intersection of a cube and a cone is empty or not:

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

Here is an illustration of the problem (the illustration is in 2D, but my problem is in 3D):

How to do that efficiently (I am searching for an algorithm, so the answer can be in C, C++ or Python) ?

Note: Here intersection is defined as: it exists a non-null 3D volume that is in the cube and in the cone (if the cube is inside the cone, or if the cone is inside the cube, they intersect).

Answer 1

For the code

This answer will be slightly more general than your problem (I consider a box instead of a cube for example). Adapting to your case should be really straightforward.

Some definitions

/*
    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];
};

Line - cone intersection

Now, let's compute the intersection between a segment and our cone. Note that I will do calculations in cone space. Note also that I take the Z axis to be the vertical one. Changing it to the Y one is left as an exercise to the reader. The line is assumed to be in cone space. The segment direction is not normalized; instead, the segment is of the length of the direction vector, and starts at the point 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;
}

Rect - cone intersection

Now, we can check whether a rectangular part of a plan intersect the cone (this will be used to check whether a face of the cube intersects the cone). Still in cone space. The plan is passed in a way that will help us: 2 vectors and a point. The vectors are not normalized, to simplify the computations.

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

Box - cone intersection

Now, the final part : the whole cube:

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]);
}

For the maths

Still in cone space, the cone equations are:

// 0 is the base, the vertex is at z = H
x² + y² = (H - z)² * tan²(alpha)
0 <= z <= H

Now, the parametric equation of a line in 3D is:

x = u + at
y = v + bt
z = w + ct

The direction vector is (a, b, c), and the point (u, v, w) lie on the line.

Now, let's put the equations together:

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

Then after developping and re-factorizing this equation, we get the following:

At² + Bt + C = 0

where A, B and C are shown in the first intersection function. Simply resolve this and check the boundary conditions on z and t.

Answer 2

1. imagine 2 infinite lines

• axis of a cone
• line going through a point P (cube center for starters) which is perpendicular to cone axis.

Cone axis is known to you so that is easy, second line is defined as

    P+t*(perpendicular vector to cone axis)

This vector can be obtained by cross product of cone axis vector and vector perpendicular to your image (assuming Z axis). The t is scalar value parameter ...

2. compute intersection of these 2 lines/axises

if you do not know the equations derive them or google them. Let the intersection point be Q

3. if intersection point Q does not lie inside cone

(between vertex and base) then point P is not intersecting cone. From intersection equations you will obtain parameters t1 and t2

• let t1 be for P axis line
• and t2 for cone axis line

if your axis line direction vector is also the cone length then intersection is inside cone if t2 = <0,1>

4. if P is not inside triangle (cutted cone to plane generated by those 2 axises)

this is also easy you know the position of Q inside cone (t2) so you know that the cone is in P-axis from Q to distance of R*t2 where R is base radius of cone. So you can compute |P-Q| and check if it is <=R*t2 or use directly t1 (if P axis direction vector is unit).

if the distance is bigger then R*t2 point P does not intersect cone.

5. if #3 and #4 are positive then P intersects cone

• hope you dont mind here is your image with few things added for clarity

[notes]

Now the hard part there are edge cases when no vertex of cube is intersecting cone but the cube itself is intersecting the cone anyway. This can occur when the ||P-Q|-R*t2| = <0,half cube size> In this case you should check more points then just cube vertexes along the closest cube face.

Another approach is:

1. create a transformation matrix for cone

Where:

• its vertex as origin
• its axis as +Z axis
• and XY plane is parallel to its base

so any point is inside cone if

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

2. convert cube vertexes into cone space

and check if any vertex is inside cone also you can check all cube edge lines with the conditions from #1 (algebraically) or use more points along cube faces as in previous method.

Chat discussion: https://chat.stackoverflow.com/rooms/48756/discussion-between-spektre-and-joojaa

Answer 3

Info: I don't know if this idea is already a patented intellectual property (in your region), or not, or how to find out, or whatever that means. I do this for fun. :)

But, here is the beef:

• Step 1: Approximation: For efficiency, treat both objects as spheres (use outer spheres). Calculate their distance (between their two center points), to find out if they are close enough to intersect at all. Quickly return false, if they can not possibly intersect (because their distance is greater than the sum of the radius of both spheres).

• Step 2: Exact Computation: This is the easy way to do it: Interpret the cone as a batch of 3-D pixels called voxels (or legos): Choose whatever resolution (granularity) you find acceptable (maybe 0.01). Create a vector pointing from (0,0,0) to any voxel point inside the volume of your cone (starting at the point you already named "vertex"). Return true, if the coordinate of that voxel exists inside your given cube. Repeat this for every voxel that you can compute for your cone object, based on the chosen granularity.

• Step 3: If nothing matches, return false.

Optimization: The function that determines if any given 3-D point is inside the cube might be optimizable, by considering the cube's inner sphere. Meaning, any given 3-D point is within the cube, if it is close enough to the center of the cube's inner sphere to be within that sphere. The fun begins, when you start filling the empty cube corners with additional spheres for more optimization (this is totally optional).

I'm sure there are further optimizations possible for step 2. However, the nice thing about this approach is that you can freely adjust the granularity, to fine tune between computation time and computation accuracy.

You can also create a solver that automatically reduces the granularity over multiple iterations. Meaning the accuracy increases over time (for better resolution).