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 At2 + 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]2
*/
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
x2 + y2 = (H - z)2 * tan2(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)2 + (v + bt)2 = (H - w - ct)2 * tan2(alpha)
Then after developping and re-factorizing this equation, we get the following:
At2 + 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.