Line-line distance

From CGAFaq

Let be points on the respective lines and unit direction vectors along the lines. Then the distance is

Often one wants the points of closest approach as well as the distance. The following is robust C code from Seth Teller that computes the these points on two 3D lines. It also classifies the lines as parallel, intersecting, or (the generic case) skew. What is listed below shows the main ideas; the full code is at http://graphics.lcs.mit.edu/~seth/geomlib/linelinecp.c.

   // computes pB ON line B closest to line A
   // computes pA ON line A closest to line B
   // return: 0 if parallel; 1 if coincident; 2 if generic (i.e., skew)
   int
   line_line_closest_points3d (
       POINT *pA, POINT *pB,           // computed points
       const POINT *a, const VECTOR *adir,     // line A, point-normal form
       const POINT *b, const VECTOR *bdir )    // line B, point-normal form
   {
   static VECTOR Cdir, *cdir = &Cdir;
   static PLANE Ac, *ac = &Ac, Bc, *bc = &Bc;
   
       // connecting line is perpendicular to both
       vcross ( cdir, adir, bdir );
   
       // check for near-parallel lines
       if ( !vnorm( cdir ) )   {
           *pA = *a;   // all points are closest
           *pB = *b;
           return 0;   // degenerate: lines parallel
           }
   
       // form plane containing line A, parallel to cdir
       plane_from_two_vectors_and_point ( ac, cdir, adir, a );
   
       // form plane containing line B, parallel to cdir
       plane_from_two_vectors_and_point ( bc, cdir, bdir, b );
   
       // closest point on A is line A ^ bc
       intersect_line_plane ( pA, a, adir, bc );
   
       // closest point on B is line B ^ ac
       intersect_line_plane ( pB, b, bdir, ac );
   
       // distinguish intersecting, skew lines
       if ( edist( pA, pB ) < 1.0E-5F )
           return 1;     // coincident: lines intersect
       else
           return 2;     // distinct: lines skew
   }