title: Line segments intersection, numerically stable test
body: <p>I need a precise and numerically stable test for 2 line segments intersection in 2D. There is one possible solution detecting 4 postions, see bellow the code.</p>
  
  <pre><code>getInters ( double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double &amp; x_int, double &amp; y_int  )
  {
      3: Intersect in two end points,
      2: Intersect in one end point,
      1: Intersect (but not in end points)
      0: Do not intersect 
  
  unsigned short code = 2;
  
  //Initialize intersections
  x_int = 0, y_int = 0;
  
  //Compute denominator
      double denom =  x1 * ( y4 - y3 ) + x2 * ( y3 - y4 ) + x4 * ( y2 - y1 ) + x3 * ( y1 - y2 ) ;
  
      //Segments are parallel
  if ( fabs ( denom ) &lt; eps)
      {
              //Will be solved later
      }
  
  //Compute numerators
      double numer1 =     x1 * ( y4 - y3 ) + x3 * ( y1 - y4 ) + x4 * ( y3 - y1 );
  double numer2 = - ( x1 * ( y3 - y2 ) + x2 * ( y1 - y3 ) + x3 * ( y2 - y1 ) );
  
  //Compute parameters s,t
      double s = numer1 / denom;
      double t = numer2 / denom;
  
      //Both segments intersect in 2 end points: numerically more accurate than using s, t
  if ( ( fabs (numer1) &lt; eps)  &amp;&amp; ( fabs (numer2) &lt; eps) || 
       ( fabs (numer1) &lt; eps)  &amp;&amp; ( fabs (numer2 - denom) &lt; eps) ||
       ( fabs (numer1 - denom)  &lt; eps)  &amp;&amp; ( fabs (numer2) &lt; eps) || 
       ( fabs (numer1 - denom) &lt; eps) &amp;&amp;  ( fabs (numer2 - denom) &lt; eps) )
      {
              code =  3;
      }
  
  //Segments do not intersect: do not compute any intersection
      else if ( ( s &lt; 0.0 ) || ( s &gt; 1 ) || 
        ( t &lt; 0.0 ) || ( t &gt; 1 ) )
      {
              return  0;
      }
  
      //Segments intersect, but not in end points
      else if ( ( s &gt; 0 ) &amp;&amp; ( s &lt; 1 ) &amp;&amp; ( t &gt; 0 ) &amp;&amp; ( t &lt; 1 ) )
      {
              code =  1;
      }
  
  //Compute intersection
  x_int = x1 + s * ( x2 - x1 );
  y_int = y1 + s * ( y2 - y1 );
  
  //Segments intersect in one end point
  return code;
   }
  </code></pre>
  
  <p>I am not sure whether all proposed conditions are designed properly (to avoid roundness errors). </p>
  
  <p>Does it make sense to use the parameters s, t for testing or use it only for the computation of an intersection?</p>
  
  <p>I am afraid that position 2 (segment intersect in one end point) may not be correctly detected (last remaining situation without any condition)...</p>