Closest points between two lines

Question

I have two arbitrary lines in 3D space, and I want to find the distance between them, as well as the two points on these lines that are closest to each other. Naturally, this only concerns the skew case, since the parallel and intersecting cases are trivial.

I know how to find the distance, as the question was asked before and answered here. I haven't found a good explanation on how to find the two points that determine that distance, though.

So specifically, given two lines

$$L_1=P_1+t_1V_1$$ $$L_2=P_2+t_2V_2$$

I would like to find two points $X_1$ on $L_1$ and $X_2$ on $L_2$ such that the distance between $X_1$ and $X_2$ is minimal.

Answer 1

The idea is that for the line segment of the shortest length, it has to be perpendicular to both the other lines.

Let the perpendicular line start from a point $P_1+t_1V_1$ of the first line, and have tangent vector $V_3$, i.e.: $$L_3=P_1+t_1V_1+t_3V_3$$

It should be that $$V_3\cdot V_2=0$$ $$V_3 \cdot V_1=0$$ which means that you can get $V_{3}$ as $$ V_3=V_2\times V_1$$

Now for this line to meet also the second line you need to have

$$P_1+t_1V_1+t_3V_3=P_2+t_2V_2$$

With this you have 3 linear equations in 3 variables, $t_1$,$t_2$ and $t_3$.

Once you solve for them, then:
- the distance between line 1 and line 2 will be $d=t_{3}/|V_{3}|$, to be taken, in case, as absolute value ($d$ is the "algebraic" distance in the direction of vector $V_{3}$);
- the closest points will of course be $Q_{1}=P_{1}+t_{1}V_{1}$ on line 1, and $Q_{2}=P_{2}+t_{2}V_{2}$ on line 2.

Answer 2

In case anyone else was, like me, looking to implement this, here is some example code in Python that follows @user2255770's answer:

from numpy import array, cross
from numpy.linalg import solve, norm

# define lines A and B by two points
XA0 = array([1, 0, 0])
XA1 = array([1, 1, 1])
XB0 = array([0, 0, 0])
XB1 = array([0, 0, 1])

# compute unit vectors of directions of lines A and B
UA = (XA1 - XA0) / norm(XA1 - XA0)
UB = (XB1 - XB0) / norm(XB1 - XB0)
# find unit direction vector for line C, which is perpendicular to lines A and B
UC = cross(UB, UA); UC /= norm(UC)

# solve the system derived in user2255770's answer from StackExchange: https://math.stackexchange.com/q/1993990
RHS = XB0 - XA0
LHS = array([UA, -UB, UC]).T
print(solve(LHS, RHS))
# prints "[ 0. -0.  1.]"

Answer 3

As it happens, under the other question (about the distance between the lines) that this question linked to, one of the answers hints at a method to find not only the shortest distance but the line along which that shortest distance lies and the intersections of that line with the two given lines. The intersection points are the points on each of the given skew lines closest to the other.

That answer used an off-site link to give the details of the solution, however. In the interest of a more self-contained answer, I'll adapt the off-site answer to your particular problem statement.

Using the notation in this question, let's write the general form of a point on line $L_1$ as $L_1(t_1)$ and a point on $L_2$ as $L_2(t_2),$ where \begin{align} L_1(t_1)&=P_1+t_1V_1, \\ L_2(t_2)&=P_2+t_2V_2. \end{align}

Now let the distance $\lVert L_1(t_1) - L_2(t_2) \rVert$ be minimized at $t_1 = r,$ $t_2 = s.$ Then the line $L_1(r)L_2(s)$ is perpendicular to both $L_1$ and $L_2,$ so a vector in the direction of that line (for example, $L_2(s) - L_1(r)$) is perpendicular to vectors in the directions of each of those lines (such as $V_1$ and $V_2$): \begin{align} (L_2(s) - L_1(r)) \cdot V_1 &= 0, \tag1 \\ (L_2(s) - L_1(r)) \cdot V_2 &= 0. \tag2 \end{align}

By writing out all the components of the vectors in Equations $(1)$ and $(2),$ multiplying componentwise (for the dot product), and collecting similar terms, we can rearrange the equations into something of the following form: \begin{align} a_1 r + b_1 s + c_1 &= 0, \\ a_2 r + b_2 s + c_2 &= 0 \end{align} where $a_1, a_2, b_1, b_2, c_1,$ and $c_2$ are all known (computed from the known components of $P_1, P_2, V_1,$ and $V_2$). So we have two unknowns ($r$ and $s$) in two simultaneous equations; solve for $r$ and $s$ in these equations, and you have your two closest points, $L_1(r)$ and $L_2(s).$

Answer 4

I work in a language that lacks a solver for systems of equations as we see in Brian's answer, so I thought I'd try a different solution. It too relies on perpendicularity, but it works by scaling the unit vector $\hat{b}$ by the length of $\hat{b}$ projected onto the $ab$ plane, or the rejection of $\hat{a}$ and $\hat{c}$ from $\hat{b}$ .

So given lines $$\vec{a_0}+\hat{a}t$$ and $$\vec{b_0}+\hat{b}t$$ and their cross product $$\hat{c}=\langle\hat{a} \times \hat{b}\rangle$$ then $$t=-\dfrac{|\vec{r}|} {b\cdot\hat{r}} $$ where $\vec{r}$ is the rejection $$\vec{r}=\vec{d} - \hat{a}(\vec{d}\cdot \hat{a}) - \hat{c}(\vec{d}\cdot \hat{c}) $$ and $\vec{d}$ is the offset $$\vec{d}=\vec{b_0}-\vec{a_0}$$

I've put together a visualization using OpenScad, which converts pretty easily to glm/GLSL:

module line(a,b=[0,0,0], width=0.05){
    hull(){ translate(a) sphere(width); translate(b) sphere(width); }
}
function dot(a,b) = a*b;
function normalize(a) = a/norm(a);

a0=[1,1,1];
a1=[5,3,2];
a=normalize(a1-a0);
b0=[3,2,0];
b1=[4,4,5];
b=normalize(b1-b0);

%line(a0, a1);
%line(b0, b1);

cn = normalize(cross(b,a));
projection_ =         dot(b0-a0, a) * a;
rejection   = b0-a0 - dot(b0-a0, a) * a - dot(b0-a0, cn) * cn;
closest_approach = b0-b*norm(rejection)/dot(b,normalize(rejection));

color("red")   line(a0, a0+projection_);
color("green") line(b0, b0-rejection);
color("blue")  line(a0+projection_, a0+projection_+cn);
color("yellow")line(b0, closest_approach);

Answer 5

Hint: \begin{align} \frac{\partial}{\partial t_1}\left.\|L_2-L_1\|^2\right|_{X_1,X_2} &= 0\\ \frac{\partial}{\partial t_2}\left.\|L_2-L_1\|^2\right|_{X_1,X_2} &= 0 \end{align}

Answer 6

Hint: I would do it by first extending both lines $L_1$ and $L_2$ to respective planes $\pi_1$ and $\pi_2$ parallel to each other. There can only be one such set of planes.

Then, proceed the normal way to find the normal vector $\mathbf{n}$ to a plane, say $\pi_1$. Now, define a new plane $\pi_3$ containing $L_1$ and $n$. Then $\pi_3$ intersects $L_2$ at only one point. That's one of the points you're looking for. After that just backtrack using $\mathbf{n}$ to find the point needed on $L_1$.