# The shortest line between two lines in 3D

Two lines in 3 dimensions generally don’t intersect at a point, they may be parallel (no intersections) or they may be coincident (infinite intersections) but most often only their projection onto a plane intersect.. When they don’t exactly intersect at a point they can be connected by a line segment, the shortest line segment is unique and is often considered to be their intersection in 3D.

The following will show how to compute this shortest line segment that joins two lines in 3D, it will as a biproduct identify parrallel lines. In what follows a line will be defined by two points lying on it, a point on line “a” defined by points P_{1} and P_{2} has an equation

P_{a} = P_{1} + mu_{a} (P_{2} –

P_{1})

similarly a point on a second line “b” defined by points

P_{4} and P_{4} will be written as

P_{b} = P_{3} + mu_{b}

(P_{4} – P_{3})

The values of mu_{a} and mu_{b} range from negative to positive infinity. The line segments between P_{1}

P_{2} and P_{3} P_{4} have their corresponding mu

between 0 and 1.

There are two approaches to finding the shortest line segment between lines “a” and “b”. The first is to write down the length of the line segment joining the two lines and then find the minimum. That is, minimize the following

|| P_{b} – P_{a} ||^{2}

Substituting the equations of the lines gives

|| P_{1} – P_{3} + mu_{a} (P_{2} –

P_{1}) – mu_{b} (P_{4} – P_{3}) ||^{2}

The above can then be expanded out in the (x,y,z) components. There are conditions to be met at the minimum, the derivative with respect to mu_{a} and mu_{b} must be zero. Note: it is easy to convince

oneself that the above function only has one minima and no other minima or maxima. These two equations can then be solved for mu_{a} and mu_{b}, the actual intersection points found by substituting the values of mu into the original equations of the line.

An alternative approach but one that gives the exact same equations is to realize that the shortest line segment between the two lines will be perpendicular to the two lines. This allows us to write two equations for the dot product as

(P_{a} – P_{b}) dot (P_{2} – P_{1}) = 0

(P_{a} – P_{b}) dot (P_{4} – P_{3}) = 0

Expanding these given the equation of the lines

( P_{1} – P_{3} + mu_{a} (P_{2} –

P_{1}) – mu_{b} (P_{4} – P_{3}) ) dot

(P_{2} – P_{1}) = 0

( P_{1} – P_{3} + mu_{a} (P_{2} –

P_{1}) – mu_{b} (P_{4} – P_{3}) ) dot

(P_{4} – P_{3}) = 0

Expanding these in terms of the coordinates (x,y,z) is a nightmare but the

result is as follows

d_{1321} + mu_{a} d_{2121} – mu_{b}

d_{4321} = 0

d_{1343} + mu_{a} d_{4321} – mu_{b}

d_{4343} = 0

where

d_{mnop} = (x_{m} – x_{n})(x_{o} –

x_{p}) + (y_{m} – y_{n})(y_{o} – y_{p})

+ (z_{m} – z_{n})(z_{o} – z_{p})

Note that d_{mnop} = d_{opmn}

Finally, solving for mu_{a} gives

mu_{a} = ( d_{1343} d_{4321} – d_{1321}

d_{4343} ) / ( d_{2121} d_{4343} – d_{4321}

d_{4321} )

and backsubstituting gives mu_{b}

mu_{b} = ( d_{1343} + mu_{a} d_{4321} )

/ d_{4343}

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 |
typedef struct { double x,y,z; } XYZ; /* Calculate the line segment PaPb that is the shortest route between two lines P1P2 and P3P4. Calculate also the values of mua and mub where Pa = P1 + mua (P2 - P1) Pb = P3 + mub (P4 - P3) Return FALSE if no solution exists. */ int LineLineIntersect( XYZ p1,XYZ p2,XYZ p3,XYZ p4,XYZ *pa,XYZ *pb, double *mua, double *mub) { XYZ p13,p43,p21; double d1343,d4321,d1321,d4343,d2121; double numer,denom; p13.x = p1.x - p3.x; p13.y = p1.y - p3.y; p13.z = p1.z - p3.z; p43.x = p4.x - p3.x; p43.y = p4.y - p3.y; p43.z = p4.z - p3.z; if (ABS(p43.x) < EPS && ABS(p43.y) < EPS && ABS(p43.z) < EPS) return(FALSE); p21.x = p2.x - p1.x; p21.y = p2.y - p1.y; p21.z = p2.z - p1.z; if (ABS(p21.x) < EPS && ABS(p21.y) < EPS && ABS(p21.z) < EPS) return(FALSE); d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z; d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z; d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z; d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z; d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z; denom = d2121 * d4343 - d4321 * d4321; if (ABS(denom) < EPS) return(FALSE); numer = d1343 * d4321 - d1321 * d4343; *mua = numer / denom; *mub = (d1343 + d4321 * (*mua)) / d4343; pa->x = p1.x + *mua * p21.x; pa->y = p1.y + *mua * p21.y; pa->z = p1.z + *mua * p21.z; pb->x = p3.x + *mub * p43.x; pb->y = p3.y + *mub * p43.y; pb->z = p3.z + *mub * p43.z; return(TRUE); } |