title: Trying to optimize line vs cylinder intersection
body: <p>My brain has been melting over a line segment-vs-cylinder intersection routine I've been working on.</p>
  
  <pre><code> /// Line segment VS &lt;cylinder&gt;
   // - 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")
   // =&gt; start = (x0, y0, z0)
   //   dir = (ux, uy, uz)
   //   A
   //   B
   //   r
   //   optimize? (= don't care for t &gt; 1)
   // &lt;= t  = "time" of intersection
   //   norm = surface normal of intersection point
   void CollisionExecuter::cylinderVSline(const Ogre::Vector3&amp; start, const Ogre::Vector3&amp; dir, const Ogre::Vector3&amp; A, const Ogre::Vector3&amp; B, const double r,
               const bool optimize, double&amp; t, Ogre::Vector3&amp; 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 &lt; B.z) { czmin = A.z - r; czmax = B.z + r; } else { czmin = B.z - r; czmax = A.z + r; }
    if (A.y &lt; B.y) { cymin = A.y - r; cymax = B.y + r; } else { cymin = B.y - r; cymax = A.y + r; }
    if (A.x &lt; B.x) { cxmin = A.x - r; cxmax = B.x + r; } else { cxmin = B.x - r; cxmax = A.x + r; }
    if (optimize) {
     if (start.z &gt;= czmax &amp;&amp; (start.z + dir.z) &gt; czmax) return;
     if (start.z &lt;= czmin &amp;&amp; (start.z + dir.z) &lt; czmin) return;
     if (start.y &gt;= cymax &amp;&amp; (start.y + dir.y) &gt; cymax) return;
     if (start.y &lt;= cymin &amp;&amp; (start.y + dir.y) &lt; cymin) return;
     if (start.x &gt;= cxmax &amp;&amp; (start.x + dir.x) &gt; cxmax) return;
     if (start.x &lt;= cxmin &amp;&amp; (start.x + dir.x) &lt; 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 &lt; 0) return;
    double time = (-b - sqrt(d)) / (2 * a);
    if (time &lt; 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() &gt; AB.length()) return; /// THIS IS THE SLOW SAFE WAY
    //if (projection.z &gt; czmax - r || projection.z &lt; czmin + r ||
    // projection.y &gt; cymax - r || projection.y &lt; cymin + r ||
    // projection.x &gt; cxmax - r || projection.x &lt; cxmin + r ) return; /// THIS IS THE FASTER BUGGY WAY
  
    normal = (intersection - projection);
    normal.normalise();
    t = time; /// at last
   }
  </code></pre>
  
  <p>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.</p>
  
  <p>Any help would be greatly appreciated!</p>
  
  <p>Cheers,</p>
  
  <p>Bill Kotsias</p>
  
  <p><strong>Edit :</strong> 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.</p>
  
  <p>Maybe, if the logic <em>is</em> correct, the problems will go away by inserting a bit of tolerance like :</p>
  
  <pre><code>  if (projection.z &gt; czmax - r + 0.001 || projection.z &lt; czmin + r - 0.001 || ... etc...
  </code></pre>