Ray/cylinder 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/cylinder intersection test. I've also provided code at the bottom as a public domain work.

Parametric surface of a cylinder

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

Ch
=
height of the cylinder
Cr
=
radius of the cylinder
Ho
=
height offset on the range [0, Ch], where 0 is aligned with an endcap
θ
=
any angle, basically
 [Equation 1:] fx( Ho , θ ) = Cr × cos(θ) fy( Ho , θ ) = Cr × sin(θ) fz( Ho , θ ) = Ho

Technically, these equations describe an infinite cylinder — one with no endcaps, which extends forever in both directions. However, we can still use them to put together the tests for a finite clinder with endcaps.

You may notice that these are very similar to the parametric equations for a cone. Mathematically speaking, a cylinder is basically a cone with infinite height, and geometrically speaking, a cylinder is just a cone that never tapers. Moreover, if we ignore the Z-component, then we're literally just writing the equation for a 2D circle.

Implicit surface of a cylinder

We have a parametric surface — that is, three parametric equations which define all points on the surface of the cylinder. In order to solve for any given part of the equation, however, we need an implicit surface: a surface defined by a single equation, wherein x, y, and z are on one side, and 0 is on the other side.

The above equations basically define a 2D circle anchored at (0, 0) on the XY-plane, and then extend it infinitely along the Z-axis. The circle's radius is Cr. This means that we can represent the X and Y equations in terms of any arbitrary 2D vector P2D:

[Equation 2:] ||P2D|| = Cr

If we project a 3D vector P onto the Z-axis, then we'll have just its Z-component. Remove that component, and we can then fit it into the equation above.

P
=
any point on the curved surface of the infinite cylinder
Ca
=
axis of the cylinder; in this case, (0, 0, 1)
Ps
=
P projected onto Ca
 Ps = (P • Ca)Ca [Equation 3:] ||P - Ps|| = Cr

If we so desire, we can now subtract Cr from both sides of the equation in order to have one side contain only zero. As such, this is now an implicit surface, which defines a cylinder that has one endcap centered on (0, 0, 0).

Handling a translated and rotated cylinder

In general, we'll want to run raycasts against a bounded cylinder — one with endcaps, rather than one that extends infinitely along its axis. There's no equation that can represent a bounded cylinder directly. However, we can still represent an infinite cylinder in terms of where the endcaps would be, and then use some distance checks later to enforce bounds.

Let's take the implicit surface of an infinite cylinder, and represent it in terms of the endcaps.

Cb
=
centerpoint of the cylinder's bottom endcap
Ct
=
centerpoint of the cylinder's top endcap
Cs
=
"spine" of the cylinder: a vector from the top to the bottom, equal to Ct - Cb
Ca
=
axis of the cylinder: normalized spine
Hp
=
hit position: a point that we presume is on the cylinder's curved surface
Hs
=
hit spine position: the hit position projected onto Cs; this position = ((Hp - Cb) • Ca)Ca + Cb
Ho
=
the height offset from our cylinder equations; this is the distance from Hs to Cb i.e. ||Hs - Cb||.

Hp fills the role of P in Equation 3, and Hs fills the role of Pp. Hs is the hit position projected onto the cylinder's axis — its "spine." If we view the cylinder from its endcap, then Hs will be at the center of a circle; if Hp is indeed on the cylinder's curved surface, then the distance from it to Hs will equal the cylinder's radius:

||Hp - Hs|| = Cr

Let's expand Hs.

||Hp - (((Hp - Cb) • Ca)Ca + Cb)|| = Cr [Equation 4:] ||Hp - Cb - ((Hp - Cb) • Ca)Ca|| = Cr

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.

Ro
=
the ray's origin
Rd
=
the ray's direction
Hp
=
the position at which the ray hits some surface
Hd
=
the hit distance: the distance from the ray origin to the hit position
 Hp = Ro + RdHd

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

||Ro + RdHd - Cb - ((Ro + RdHd - Cb) • Ca)Ca|| = Cr

We can simplify our problem if we translate the entire coordinate system — move everything so that the base of the cylinder is (0, 0, 0). We'll do that by defining a new variable.

Rl
=
the ray's origin relative to the cylinder's base position = Ro - Cb
info
||RlRo + RdHd - Cb - ((RlRo + RdHd - Cb) • Ca)Ca|| = Cr

Now, we can continue working to extricate Hd from the other terms.

info
(Rl + RdHd - ((Rl + RdHd) • Ca)Ca) • (Rl + RdHd - ((Rl + RdHd) • Ca)Ca) = Cr squared 
info
(Rl • Rl) + (RdHd • RdHd) + (((Rl + RdHd) • Ca)Ca • ((Rl + RdHd) • Ca)Ca) + 2(Rl • RdHd) - 2(Rl • ((Rl + RdHd) • Ca)Ca) - 2(RdHd • ((Rl + RdHd) • Ca)Ca) = Cr squared 
info
(Rl • Rl) + (RdHd • RdHd) + (((Rl + RdHd) • Ca) squared ) + 2(Rl • RdHd) - 2(Rl • ((Rl + RdHd) • Ca)Ca) - 2(RdHd • ((Rl + RdHd) • Ca)Ca) = Cr squared 
info
(Rl • Rl) + (RdHd • RdHd) + (((Rl + RdHd) • Ca) squared ) + 2(Rl • RdHd) - 2(Rl • ((Rl + RdHd) • Ca)Ca) - 2(Rd • Ca)((Rl + RdHd) • Ca)Hd = Cr squared 
info
(Rl • Rl) + (RdHd • RdHd) + (((Ca • Rl) + (Ca • RdHd)) squared ) + 2(Rl • RdHd) - 2(Rl • ((Ca • Rl) + (Ca • RdHd))Ca) - 2(Rd • Ca)((Ca • Rl) + (Ca • RdHd))Hd = Cr squared 
info
(Rl • Rl) + (RdHd • RdHd) + (((Ca • Rl) + (Ca • RdHd)) squared ) + 2(Rl • RdHd) - 2(Rl • ((Ca • Rl) + (Ca • RdHd))Ca) - 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • RdHd)Hd = Cr squared 
info
(Rl • Rl) + (Rd • Rd)Hd squared  + (((Ca • Rl) + (Ca • Rd)Hd) squared ) + 2(Rd • Rl)Hd - 2(Rl • ((Ca • Rl) + (Ca • Rd)Hd)Ca) - 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • Rd)Hd squared  = Cr squared 
info
(Rl • Rl) + (Rd • Rd)Hd squared  + (((Ca • Rl) + (Ca • Rd)Hd) squared ) + 2(Rd • Rl)Hd - 2(Ca • Rl)((Ca • Rl) + (Ca • Rd)Hd) - 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • Rd)Hd squared  = Cr squared 
info
(Rl • Rl) + (Rd • Rd)Hd squared  + (Ca • Rl) squared  + (Ca • Rd) squared Hd squared  + 2(Ca • Rd)(Ca • Rl)Hd + 2(Rd • Rl)Hd - 2(Ca • Rl)(Ca • Rl) - 2(Ca • Rl)(Ca • Rd)Hd - 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • Rd)Hd squared  = Cr squared 
info
(Rd • Rd)Hd squared  + (Ca • Rd) squared Hd squared  - 2(Ca • Rd)(Ca • Rd)Hd squared  + 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • Rl)Hd + 2(Rd • Rl)Hd + (Ca • Rl) squared  - 2(Ca • Rl) squared  + (Rl • Rl) = Cr squared 
info
((Rd • Rd) + (Ca • Rd) squared  - 2(Ca • Rd) squared (Ca • Rd) squared )Hd squared  + 2(Ca • Rd)(Ca • Rl)Hd - 2(Ca • Rd)(Ca • Rl)Hd + 2((Rd • Rl) - (Ca • Rd)(Ca • Rl))Hd + (Ca • Rl) squared  - 2(Ca • Rl) squared  - (Ca • Rl) squared  + (Rl • Rl) = Cr squared  (1 - (Ca • Rd) squared )Hd squared  + 2((Rd • Rl) - (Ca • Rd)(Ca • Rl))Hd - (Ca • Rl) squared  + (Rl • Rl) = Cr squared 

If we arrange that lengthwise, it'll be easier to see that as with ray/cone intersections, we've ended up with a quadratic equation. We can use the quadratic formula to find its roots — its solutions; the hit distances.

ax squared  + bx + c = 0 [Equation 5:] (1 - (Ca • Rd) squared )Hd squared  + 2((Rd • Rl) - (Ca • Rd)(Ca • Rl))Hd + (Rl • Rl) - (Ca • Rl) squared  - Cr squared  = 0 ∴
a
=
1 - (Ca • Rd) squared 
b
=
2((Rd • Rl) - (Ca • Rd)(Ca • Rl))
c
=
(Rl • Rl) - (Ca • Rl) squared  - Cr squared 
Hd = ( -b ± √(b squared  - 4ac) ) ÷ 2a

In the quadratic formula, the expression b squared  - 4ac is called the discriminant. If its value is negative, then there is no valid solution — no intersection between the ray and the curved surface of our cylinder. If its value is zero, then there's exactly one solution: you can apply the formula either way (plus or minus) and you'll get the same result. If the discriminant is a positive non-zero number, then there are two different solutions — two different points at which we have potential intersections.

Yes. "Potential." This equation actually tests for intersections between an infinite double-sided line, and an infinitely long cylinder with no endcaps. In order to test for the intersection between a ray and a bounded cylinder, we need to add some additional inequalities, and require that they be true. Let's start with one to rule out any hit positions behind the ray origin:

Hd ≥ 0

Simple enough. Next up, we need to apply bounds to our cylinder. We'll take the vector from the hit position to one of the cylinder's endcaps, project that onto the cylinder's axis, and check the resulting length against the cylinder's height.

0 ≤ (Ct - Hp) • Ca ≤ ||Ct - Cb||

There's one more wrinkle: we need to test the ray against the cylinder's endcaps. That'll be two ray/disc intersection checks. However, since both discs exist on parallel planes, we can actually combine some of the math. Both discs have the same normal vectors, so when we take the dot product of the ray direction with the plane normal, we can reuse that across both ray/disc checks.

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_cylinder_intersection(
   const glm::vec3& ray_origin,
   const glm::vec3& ray_direction, // must be normalized

   const glm::vec3& cylinder_endcap_t, // centerpoint of a disc on one end of the cylinder
   const glm::vec3& cylinder_endcap_b, // centerpoint of a disc on one end of the cylinder
   const float      cylinder_radius,   // radius of both discs

   const bool hits_from_inside_count,

   float& hit_distance
) {
   //
   // First, identify intersections between a line and an infinite cylinder. An infinite 
   // cylinder has no base and extends in both directions.
   //
   glm::vec3 Rl = ray_origin - cylinder_endcap_b;        // Ray origin local to centerpoint
   glm::vec3 Cs = cylinder_endcap_t - cylinder_endcap_b; // Cylinder spine
   float     Ch = length(Cs);                            // Cylinder height
   glm::vec3 Ca = Cs / Ch;                               // Cylinder axis

   auto Ca_dot_Rd = dot(Ca, ray_direction);
   auto Ca_dot_Rl = dot(Ca, Rl);
   auto Rl_dot_Rl = dot(Rl, Rl);

   float a = 1 - (Ca_dot_Rd * Ca_dot_Rd);
   float b = 2 * (dot(ray_direction, Rl) - Ca_dot_Rd * Ca_dot_Rl);
   float c = Rl_dot_Rl - Ca_dot_Rl * Ca_dot_Rl - (cylinder_radius * cylinder_radius);

   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 cylinder that matches our finite cylinder. This means that we cannot 
      // be hitting any part of the cylinder: 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 cylinder.
      //
      return false;
   }
   //
   // Now, we need to take our intersection points and ensure that they lie on the 
   // surface of a finite cylinder. If one of them is past the edges of the finite 
   // cylinder, then we need to check for a valid intersection with the endcaps.
   //
   if (count > 2) {
      cobb::unreachable();
   }
   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 = dot(cylinder_endcap_t - Hp1, Ca); // height offset
   float     Ho2 = dot(cylinder_endcap_t - 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) {
      //
      // The ray never hits the bounded cylinder's curved surface. If we're looking 
      // along the cylinder's axis -- whether from inside or outside -- then the ray 
      // could still hit an endcap.
      // 
      // Let's project the ray origin onto the cylinder's axis, and figure out which 
      // endcap we're nearer to. (Well, actually, we already have that value: it's 
      // Ca_dot_Rl.)
      //
      if (Ca_dot_Rl <= 0.0) { // above
         valid1 = ray_disc_intersection(ray_origin, ray_direction, cylinder_endcap_t, Ca, cylinder_radius, hit_near);
      } else if (Ca_dot_Rl >= Ch) { // below
         valid1 = ray_disc_intersection(ray_origin, ray_direction, cylinder_endcap_b, Ca, cylinder_radius, hit_away);
      } else {
         //
         // Inside. Test both discs.
         //
         if (!hits_from_inside_count) {
            return false;
         }
         valid1 = ray_disc_intersection(ray_origin, ray_direction, cylinder_endcap_t, Ca, cylinder_radius, hit_near);
         valid2 = ray_disc_intersection(ray_origin, ray_direction, cylinder_endcap_b, Ca, cylinder_radius, hit_away);
         if (valid1) {
            if (valid2) {
               if (hit_away < hit_near)
                  hit_near = hit_away;
            }
            hit_distance = hit_near;
            return true;
         } else if (valid2) {
            hit_distance = hit_away;
            return true;
         }
      }
      if (valid1) {
         hit_distance = hit_near;
         return true;
      }
      return false;
   }
   if (valid_count == 1) {
      //
      // The ray hits the cylinder's curved surface only once. This can only happen under 
      // two cases: the ray originates from inside the cylinder, and points outward; or 
      // the ray passes through the bounded cylinder once and then through an endcap.
      //
      if (valid2) {
         Hp1 = Hp2;
         Ho1 = Ho2;
         valid1   = true;
         hit_near = hit_away;
      }

      float disc_near;
      float disc_away;
      bool  disc1 = ray_disc_intersection(ray_origin, ray_direction, cylinder_endcap_t, Ca, cylinder_radius, disc_near);
      bool  disc2 = ray_disc_intersection(ray_origin, ray_direction, cylinder_endcap_b, Ca, cylinder_radius, disc_away);
      if (disc1) {
         if (disc2) {
            if (disc_away < disc_near)
               disc_near = disc_away;
         }
      } else if (disc2) {
         disc_near = hit_away;
         disc1 = disc2;
      }

      if (disc1) {
         if (disc_near < hit_near) {
            hit_distance = disc_near;
            return true;
         }
      } else {
         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;
         }
         hit_distance = hit_near;
         return true;
      }
   }
   hit_distance = std::min(hit_near, hit_away);
   return true;
}