c++ - Trying to optimize line vs cylinder intersection

Question

My brain has been melting over a line segment-vs-cylinder intersection routine I've been working on.

 /// Line segment VS <cylinder>
 // - cylinder (A, B, r) (start point, end point, radius)
 // - line has starting point (x0, y0, z0) and ending point (x0+ux, y0+uy, z0+uz) ((ux, uy, uz) is "direction")
 // => start = (x0, y0, z0)
 //   dir = (ux, uy, uz)
 //   A
 //   B
 //   r
 //   optimize? (= don't care for t > 1)
 // <= t  = "time" of intersection
 //   norm = surface normal of intersection point
 void CollisionExecuter::cylinderVSline(const Ogre::Vector3& start, const Ogre::Vector3& dir, const Ogre::Vector3& A, const Ogre::Vector3& B, const double r,
             const bool optimize, double& t, Ogre::Vector3& normal) {
  t = NaN;

  // Solution : http://www.gamedev.net/community/forums/topic.asp?topic_id=467789
  double cxmin, cymin, czmin, cxmax, cymax, czmax;
  if (A.z < B.z) { czmin = A.z - r; czmax = B.z + r; } else { czmin = B.z - r; czmax = A.z + r; }
  if (A.y < B.y) { cymin = A.y - r; cymax = B.y + r; } else { cymin = B.y - r; cymax = A.y + r; }
  if (A.x < B.x) { cxmin = A.x - r; cxmax = B.x + r; } else { cxmin = B.x - r; cxmax = A.x + r; }
  if (optimize) {
   if (start.z >= czmax && (start.z + dir.z) > czmax) return;
   if (start.z <= czmin && (start.z + dir.z) < czmin) return;
   if (start.y >= cymax && (start.y + dir.y) > cymax) return;
   if (start.y <= cymin && (start.y + dir.y) < cymin) return;
   if (start.x >= cxmax && (start.x + dir.x) > cxmax) return;
   if (start.x <= cxmin && (start.x + dir.x) < cxmin) return;
  }

  Ogre::Vector3 AB = B - A;
  Ogre::Vector3 AO = start - A;
  Ogre::Vector3 AOxAB = AO.crossProduct(AB);
  Ogre::Vector3 VxAB  = dir.crossProduct(AB);
  double ab2 = AB.dotProduct(AB);
  double a = VxAB.dotProduct(VxAB);
  double b = 2 * VxAB.dotProduct(AOxAB);
  double c = AOxAB.dotProduct(AOxAB) - (r*r * ab2);
  double d = b * b - 4 * a * c;
  if (d < 0) return;
  double time = (-b - sqrt(d)) / (2 * a);
  if (time < 0) return;

  Ogre::Vector3 intersection = start + dir * time;        /// intersection point
  Ogre::Vector3 projection = A + (AB.dotProduct(intersection - A) / ab2) * AB; /// intersection projected onto cylinder axis
  if ((projection - A).length() + (B - projection).length() > AB.length()) return; /// THIS IS THE SLOW SAFE WAY
  //if (projection.z > czmax - r || projection.z < czmin + r ||
  // projection.y > cymax - r || projection.y < cymin + r ||
  // projection.x > cxmax - r || projection.x < cxmin + r ) return; /// THIS IS THE FASTER BUGGY WAY

  normal = (intersection - projection);
  normal.normalise();
  t = time; /// at last
 }

I have thought of this way to speed up the discovery of whether the projection of the intersection point lies inside the cylinder's length. However, it doesn't work and I can't really get it because it seems so logical : if the projected point's x, y or z co-ordinates are not within the cylinder's limits, it should be considered outside. It seems though that this doesn't work in practice.

Any help would be greatly appreciated!

Cheers,

Bill Kotsias

Edit : It seems that the problems rise with boundary-cases, i.e when the cylinder is parallel to one of the axis. Rounding errors come into the equation and the "optimization" stops working correctly.

Maybe, if the logic is correct, the problems will go away by inserting a bit of tolerance like :

  if (projection.z > czmax - r + 0.001 || projection.z < czmin + r - 0.001 || ... etc...

Answer 1

The cylinder is circular, right? You could transform coordinates so that the center line of the cylinder functions as the Z axis. Then you have a 2D problem of intersecting a line with a circle. The intersection points will be in terms of a parameter going from 0 to 1 along the length of the line, so you can calculate their positions in that coordinate system and compare to the top and bottom of the cylinder.

You should be able to do it all in closed form. No tolerances. And sure, you will get singularities and imaginary solutions. You seem to have thought of all this, so I guess I'm not sure what the question is.