Ray/cone intersections

There are a lot of people out on the web who are willing to share raycasting code, but not a whole lot of them bother to actually apply a permissive license to it so that people can legally use it. What's more, a lot of them don't explain how they actually derived the math involved either.

Here, I'm going to walk through how to put together a ray/cone intersection test. I've also provided code at the bottom as a public domain work.

Parametric surface of a cone

The parametric equations (that is, equations with parameters) for a cone are as follows:

Technically, these equations describe an infinite cone — one with no endcap, which extends forever in both directions. (Yes, both. When you go above the tip, the cone turns inside out and starts widening again.) However, we can still use them to put together the tests for a finite cone with an endcap.

Implicit surface of a cone

At the top of this page, we define a cone as a parametric surface — that is, we have three parametric equations (equations that take parameters) which define all points on the surface of the cone. However, if we can come up with an implicit surface — that is, a surface defined by a single equation, wherein x, y, and z are on one side, and 0 is on the other side — then we can use vector math to deal with the cone more easily.

We know from our exploration above that at the cone's base, a point lies on the cone's surface (the edge of the base) if Equation 3 is true. We can represent that equation in terms of any arbitrary 2D vector P:

We also know that as we move further along the cone's height, to its tip, the radius shrinks proportional to the cone's height. Therefore:

And of course, any equation of the form "a = b" is implicitly also "a - b = 0".

Handling a translated and rotated cone

The cone's local z-axis is its main axis, while its local x- and y-axes are its cross axes. The above relationship establishes the distance that any point on the cone's surface must have along the cross axes, given the point's distance from the base along the main axis.

Let's define some variables real quick:

H_s is what we get when we project the hit position H_p onto the cone's axis — its "spine." In other words, we have flattened the cone from 3D into 2D, and are looking at it from the spine: this is the circle graph from above. H_s is the nearest point on the cone's "spine" to H_p, so the distance between the two points is the distance from the hit position to the "spine" as a whole, and therefore also the 2D distance from the hit position to the circle's (cone's) center.

We can use two equations to check whether the hit position is on the outer edge of the cone's base. The first equation is a check between the distance and the radius; recall Equations 3 and 4. The second equation constrains this check to 2D, by requiring that there be no distance between the projected hit position and the plane we've projected it onto (otherwise it would be a point/sphere test).

Meanwhile the distance from H_s to the center of the cone's base, C_b, is the height offset: the radius, the distance, that we want to check will be proportional to this. Recall Equation 2, which established the exact relationship between the height offset and radius; we can substitute this into our equation to go from 2D to 3D:

We can double-check this with some simple mental tests.

• If our hit position is at the bottom of the cone, then H_o is zero. The righthand side resolves to a division of the cylinder's height by itself, yielding 1, which is multiplied by the cylinder's radius. Ultimately, this resolves out to a 2D radius check, as expected.
• If our hit position is at the top of the cone, then H_o = C_h. The righthand side's numerator comes out to zero, so the radius at the top of the cone is zero: a single point. This, too, is as expected.
• If our hit position is halfway along the cone's height, then the numerator on the righthand side computes to half of the cone's height. Dividing that by the cylinder's height gives us 0.5, so we end up comparing to half of the cone's radius at its base, as expected.

Putting content on the righthand side in terms of vectors (except for the height, because leaving that as-is will make things easier later on):

Next, we'll divide both sides by ||C_t - H_s|| so that C_r is the lone numerator on the right. (Mathematical notation sort of obscures this fact, visually, but on the righthand side, because C_r is multiplied by the entire fraction, it is effectively part of the righthand-side numerator. We can divide by the rest of the righthand-side numerator to get the radius on its own.)

Let's split up that lefthand-side fraction.

This will be a little easier if we separate the cone's height-to-radius ratio out into a different variable, and if we square our equations. Let's do both at the same time:

I think our new Equation 5 would, technically, be enough to generate a ray/cone intersection solution. However, it's a bit messy: the lefthand side uses two points, one of which is the point we want to solve for (a raycast hit position) and the other of which is itself defined in terms of the point we want to solve for. We should try to fix that — make it so that one side of the equation contains only one unknown point (the hit position), possibly alongside other known points (i.e. points that "belong to" the cone).

Right now, the equation works in terms of three points in space: H_p, H_s, and C_t. These points form a triangle, and our equation describes the lengths of two edges: H_p to H_s, and C_t to H_s. We therefore need to find some property that is inherent to all triangles, that we can use to rearrange our equation; we need one side of the equation to be the H_p to C_t edge.

Let's break out the handy-dandy Pythagorean theorem, then, which describes the relationship between the lengths of a triangle's edges:

And therefore:

We can redefine ||H_p - H_s||^2 squared , which takes the role of c^2 squared  in the Pythagorean theorem:

Let's substitute this whole formula in, replacing ||H_p to H_s||^2 squared :

We can add ||H_p to H_s||^2 squared  to both sides, to reduce the lefthand side down to just one length:

This equation will lead to a cleaner solution. If we substitute the equation for a ray in place of H_p, the ray hit position, then we can solve for H_p to compute that hit position given any ray. However, there's still more cleanup we can do. Let's also try to replace H_s with something in terms of H_p. Recall that H_s is the projection of H_p onto the axis C_a.

Remember that H_s is the hit position projected onto the cylinder's "spine." To compute that, we start with H_p - C_b, to get the hit position relative to the cylinder's base. We then project that onto the cylinder's axis vector:

And then finally, we add the cylinder's base position back in, to move us back to world-relative coordinates:

With that in mind, let's see if we can simplify the righthand side of Equation 6. Let's start by expanding H_s:

Now, if we think about it, the cylinder's tip is really just the cylinder's base, plus the cylinder's height times the cylinder's axis:

Within the righthand side, we can strip out the additions of the base on both sides. (The righthand side is taking the distance between two points; removing those additions means translating both points by the same vector, which has no effect on the distance between them and therefore is a valid way to simplify the equation.)

We can also drop the multiplications by C_a. Those multiplications were of lengths by an axis vector, to scale along that axis; dropping the multiplications means that we're working with raw lengths, so we go from a ||length operator|| to parentheses:

Finally, we can flip the righthand side: instead of measuring the projected distance from H_p to the bottom of the cylinder and then subtracting that from the cylinder's height, we can instead measure the projected distance from H_p to the top of the cylinder.

Our final cone equation, substituting out ||C_t - H_s||^2 squared  from Equation 6, is as follows:

Keep in mind that this equation only includes the variable for the cone's tip, not the cone's base. We have effectively removed most of our vectors from the equation, and so just like with the parametric equations, this implicit equation tests whether a point lies on the curved outer surface of an infinitely tall cone anchored at its tip. That's fine, though. Testing a finite cone just involves a distance check between the point and the tip; and testing against the cone's flat endcap involves a point/disc intersection check. Our equation here tests the "unique" part of the cone.

Solving for the ray

A ray is commonly defined in terms of a starting point, or origin, and a normalized direction vector. We can therefore define the hit position in terms of the same.

We can substitute this equation into our cone equation above, and then attempt to solve for the hit distance H_d:

On both sides of the equation, we subtract C_t from R_o. Let's simplify the equation a bit by defining that difference as a vector:

The squared length of a vector is equal to the dot product of the vector with itself; we can use that to expand the lefthand side:

We can also expand the righthand side if we keep in mind that dot products obey FOIL. The squared parentheses on the righthand side only use three terms, nor four, but we can just add a zero term and then use FOIL.

So the ray-and-cone equation, with the righthand side expanded, is:

We can expand the dot product on the lefthand side using FOIL as well:

We can expand the squaring on the righthand side as well:

And finally: for any vector V and scalar s:

Now, let's group them based on H_d. We'll start by subtracting one side of the equation from the other:

Next, we'll unpack the multiplication by (C_q + 1), and then merge the groups that are both multiplied by 2:

Next, we want to move H_d outside of some of the dot products. Multiplying a vector by a scalar is the same as taking the dot product of that vector with another vector whose terms are all that same scalar; and the dot product between two vectors is commutative: you can swap the lefthand and righthand sides around and get the same result; therefore, for any two vectors A and B, and a scalar c,

This is a quadratic equation. Can you see it? This is a quadratic equation. Here — I'll mark the different parts in different colors and show you.

If H_d is the "x" in our quadratic equation, then...

And of course, per the quadratic formula, the quadratic equation is true whenever H_d is defined as follows:

The discrminant can be used to check how many solutions exist — that is, whether the ray intersects the cone, and if so, how many times — and is defined as follows:

If the discriminant is negative, then no solutions exist; there is no intersection. If the discriminant is positive, then there are two intersections. If the discriminant is zero, then there is only one intersection; if you compute the two intersection points (differing based on ± the square root of the discriminant), then they will be identical.

Recall that for any intersection, H_p = R_o + R_dH_d. Our process above tests for an intersection against an infinite cone, so we now need to constrain it.

The following expression will be true if a hit position is on the finite cone. We are taking the hit position's distance to the cone's tip, projecting that onto the cone's axis to get the hit position's cone-relative height offset, and testing that projected distance against the cone's height.

The following expression will be true if the hit position is below the finite cone. In that case, you should test for a ray/disc intersection with the cone's endcap.

The following expression will be true if the hit position is above the finite cone. In this situation, you have to ignore the hit position. Remember that odd thing I pointed out at the start: our equations don't describe just one cone. The cone isn't bounded on either side, so if you go above the tip, the cone turns inside out and begins to widen again. Essentially, the equations describe something akin to an hourglass shape; so, we have to ignore everything above the tip where the normal cone and the inside-out cone meet.

Lastly, we also have to ignore any hit positions that result from a negative hit distance H_d, as these represent points behind the ray.

Free code!

The following C++ code is licensed under CC0 (full legal text / summary), and so is effectively a public domain work. You can incorporate it into your programs without the need for attribution, payment, or similar. (That said, linking back here would be polite!)

This code makes use of the GLM library for vector math. It has transitive dependencies on some of my other math, in the form of ray_disc_intersection and, through it, ray_plane_intersection.

//
// Given coefficients in a quadratic equation, this function gives you the roots 
// and returns the number of roots. If there is only one root, then both root 
// variables are set to the same value.
//
int quadratic_roots(
   const float a,
   const float b,
   const float c,
   //
   float& root_lower,
   float& root_upper
) {
   constexpr auto EPSILON = 1e-8;

   float discriminant = (b * b) - (4.0 * a * c);
   if (discriminant > EPSILON) {
      float b_term = b < EPSILON ? -b + sqrt(discriminant) : -b - sqrt(discriminant);

      root_lower = b_term / (2.0 * a); // quadratic formula
      root_upper = (2.0 * c) / b_term; // citardauq formula

      if (root_lower > root_upper)
         std::swap(root_lower, root_upper); // use of both formulae, plus this, avoids catastrophic cancellation due to floating-point limits

      return 2;
   } else if (discriminant > -EPSILON && discriminant <= EPSILON) {
      root_lower = -(b / 2.0 * a);
      root_upper = root_lower;
      return 1;
   }
   root_lower = NAN;
   root_upper = NAN;
   return 0;
}
   
bool ray_cone_intersection(
   const glm::vec3& ray_origin,
   const glm::vec3& ray_direction, // must be normalized
   
   const glm::vec3& cone_tip,
   const glm::vec3& cone_base_centerpoint,
   float cone_base_radius,

   const bool hits_from_inside_count,
   
   float& out_hit_distance
) {
   //
   // First, identify intersections between a line and an infinite cone. An infinite 
   // cone has no base and extends in both directions -- downward from the tip as you 
   // would expect, but also inside-out upward from the tip.
   //
   glm::vec3 Rl = ray_origin - cone_tip;            // Ray origin local to centerpoint
   glm::vec3 Cs = cone_tip - cone_base_centerpoint; // Cone spine
   float     Ch = glm::length(Cs);                  // Cone height
   glm::vec3 Ca = Cs / Ch;                          // Cone axis
   float     Cq = (cone_base_radius * cone_base_radius) / (Ch * Ch); // Cone ratio

   auto Rd_dot_Ca = glm::dot(ray_direction, Ca);
   auto Rl_dot_Rl = glm::dot(Rl, Rl);
   auto Rl_dot_Ca = glm::dot(Rl, Ca);

   float a = glm::dot(ray_direction, ray_direction) - (Cq + 1) * (Rd_dot_Ca * Rd_dot_Ca);
   float b = 2 * (glm::dot(ray_direction, Rl) - (Cq + 1) * Rl_dot_Ca * Rd_dot_Ca);
   float c = Rl_dot_Rl - (Cq + 1) * Rl_dot_Ca * Rl_dot_Ca;

   float hit_near;
   float hit_away;
   auto  count = quadratic_roots(a, b, c, hit_near, hit_away);
   if (count == 0) {
      //
      // There is no intersection between a line (i.e. a "double-sided" ray) and the 
      // infinite cone that matches our finite cone. This means that we cannot be 
      // hitting any part of the cone: if we were hitting the base from the inside, 
      // for example, then the "back of our ray" would be hitting the upper part of 
      // the infinite cone.
      //
      return false;
   }
   //
   // Now, we need to take our intersection points and ensure that they lie on the 
   // surface of a finite cone. If one of them is below the finite cone, then we need 
   // to check for a valid intersection with the cone's base.
   //
   if (count > 2) {
      std::unreachable(); // requires C++23 or a polyfill
   }
   bool  valid1 = true;
   bool  valid2 = true;
   //
   glm::vec3 Hp1 = ray_origin + ray_direction * hit_near;
   glm::vec3 Hp2 = ray_origin + ray_direction * hit_away;
   float     Ho1 = glm::dot(cone_tip - Hp1, Ca); // height offset
   float     Ho2 = glm::dot(cone_tip - Hp2, Ca);
   //
   int valid_count = count;
   if (hit_near < 0.0 || Ho1 < 1.0e-8 || Ho1 > Ch) {
      valid1 = false;
      --valid_count;
   }
   if (hit_away < 0.0 || Ho2 < 1.0e-8 || Ho2 > Ch) {
      valid2 = false;
      if (count > 1)
         --valid_count;
   }
   //
   if (valid_count == 0) {
      if (hits_from_inside_count) {
         //
         // The ray never hits the bounded cone's curved surface. If it originates 
         // from inside the cone and points "downward," however, it could still hit 
         // the cone's endcap from inside.
         //
         if (ray_disc_intersection(ray_origin, ray_direction, cone_base_centerpoint, Ca, cone_base_radius, hit_away)) {
            out_hit_distance = hit_away;
            return true;
         }
      }
      return false;
   }
   if (valid_count == 1) {
      //
      // The ray hits the cone's curved surface only once. This can only happen under 
      // two cases: the ray originates from inside the cone, and points outward; or 
      // the ray passes through the bounded cone once and through the cone's endcap.
      //
      if (valid2) {
         Hp1 = Hp2;
         Ho1 = Ho2;
         valid1   = true;
         hit_near = hit_away;
      }
      valid2 = ray_disc_intersection(ray_origin, ray_direction, cone_base_centerpoint, Ca, cone_base_radius, hit_away);
      //
      if (!valid2) {
         if (!hits_from_inside_count) {
            //
            // The ray only hits the finite cone's curved surface once, and never hits 
            // the cone's base. This means that the ray must be originating from inside 
            // of the cone.
            //
            return false;
         }
         out_hit_distance = hit_near;
         return true;
      }
   }
   out_hit_distance = std::min(hit_near, hit_away);
   return true;
}