The near intersection of two skew lines

Sometimes a problem arises where two lines in space should intersect but due to small measurement errors they are skew. The following figure shows the near intersection of two such skew lines. There is a point on each line that is closest to the other line. The midpoint of the segment connecting these points gives a good estimate of the desired intersection point.

The desired intersection point is labeled M and shown in red.

Line 1 is defined by points A and B, with unit vector U1 pointing along the line. Line 2 is defined by points C and D, with unit vector U2 pointing along the line.

A vector, W, perpendicular to both lines is given by W = U1 cross U2. If the lines are parallel then the length of W is 0 and there is no one point of closest approach.

Assuming non-parallel lines, a unit vector, U3, is given by
U3 = W/|W| = (U1 cross U2)/|U1 cross U2|. This unit vector may point from Line 1 to Line 2, or in the opposite direction, depending on the directions of U1 and U2.

A displacement vector from Line 2 to Line 1 may be found as follows. Let the vector from point C to point A be V = A-C. The length of V along U3 is (V dot U3) and is > 0 if U3 is parallel to V, and < 0 if U3 is anti-parallel to V. This length is also the distance between the lines. Since the sign of V reverses if U3 reverses, the product is a vector that always points the same way, from Line 2 to Line 1. So the displacement vector, E, needed to shift Line 2 to intersect Line 1 is given by
E = (V dot U3)U3 = ((A-C) dot U3)U3.

The shifted line is shown in the following figure.

Line 2 is shifted as follows:
C' = C + E
D' = D + E

Point I' is found using the algorithm for the intersection of two co-planar lines.

Then I = I' - E/2.

The distance by which the two lines miss each other is |E|, which may be used to check for abnormally large errors.

Algorithm

Unit vector U1 along Line 1:
U1 = (B-A)/| B-A |
Unit vector U2 along Line 2:
U2 = (D-C)/| D-C |
Unit vector U3 from Line 2 toward Line 1:
U3 = (U1 cross U2)/| U1 cross U2 |
Displacement vector from Line 2 to Line 1:
E = ((A-C) dot U3) U3
Line 2':
C' = C + E
D' = D + E
Unit vector U4 perpendicular to Line 2':
V = A - C'
V_perp = V - (V dot U2) U2
U4 = V_perp/| V_perp |
Similar triangle side lengths P and p:
P = (A - C') dot U4
p = (B - C') dot U4
Point I' fraction of distance from point A to point B:
Frac = H/(H+h) = P/(P+p)
Intersection point I':
I' = A + U1 (P/(P+p))
Intersection point I:
I = I' - E/2

In IDL:
  U1 = unit(B-A)
  U2 = unit(D-C)
  U3 = unit(cross(u1,u2))
  E = total((A-C)*U3)*U3
  C2 = C + E
  D2 = D + E
  V = A - C2
  U4 = unit(V-total(V*U2)*U2)
  P1 = total((A-C2)*U4)
  P2 = total((B-C2)*U4)
  I2 = A + U1*(P1/(P1+P2))
  I = I2 - E/2

  The JHU/APL/S1R IDL library routine SKEWINT gives the
  near-intersection point for two skew lines.  The call is 
  I = SKEWINT(A, B, C, D, [E]) where
  A, B, C, and D are the 3 element arrays representing the
  line endpoints, E is the optionlly returned distance between
  the lines, and I is the near-intersection point.

The near intersection of two skew lines

Terms