# Closest points between two lines

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.

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.

• Is it not $d = t_3 |V_3|$ rather? Feb 27, 2018 at 20:21
• Going through the answer from @user2277550, I found an error. Distance between the lines should be computed as: d = ||t3 * **V3**||. Otherwise, the equation ![enter image description here](i.stack.imgur.com/NNyA9.png) makes no sense.
– Mike
Oct 21, 2019 at 11:27

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.]"

• How can this be done without numpy? Jun 15 at 20:16

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.

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).$

• details for this approach are given here. Jan 3 at 0:38

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}|} {\hat{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);


• could you please explain what "the rejection" is? Dec 29, 2021 at 11:16
• the "rejection" is the standard vector rejection described on wikipedia: en.wikipedia.org/wiki/Vector_projection. In this case I'm rejecting both a and c. You can think of it as expressing the offset d in terms of a coordinate system that uses a, b, and c as basis vectors, then isolating the component that's defined using b. Jan 1 at 19:42
• Checking out this code in OpenSCAD was way easier than I expected! If anyone else is struggling with this I highly recommend playing around with the vectors in OpenSCAD. Feb 26 at 18:44
• I found an edge case where the rejection is 0 and it causes a divide by 0 error: a=[0,0,0] a1=[1,0,0] b0=[1,1,0] b1=[1,1,1] Feb 26 at 20:38
• What is $\vec{a_0}$? Jun 15 at 20:22

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}

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$.

Here is a GDScript implementation similar to @16807's response.

func line_intersection(line_position, line_normal, ray_position, ray_normal):
var pos_diff = line_position - ray_position
var cross_normal = line_normal.cross(ray_normal).normalized()
var rejection = pos_diff - pos_diff.project(ray_normal) - pos_diff.project(cross_normal)
var distance_to_line_pos = rejection.length() / line_normal.dot(rejection.normalized())
var closest_approach = line_position - line_normal * distance_to_line_pos
return closest_approach


Some notes:

• All arguments are Vector3, and the function returns Vector3.
• The lines are defined by position vectors and normal vectors.
• The resulting closest_approach point will be on the lined defined by the line_ parameters.
• In my case, I am working with a ray from the mouse into the screen, so I used the term ray_ as the naming convention for the secondary line for readability. However, this is treated the same as an infinite line in the math.
• The name of the variable distance_to_line_pos is slightly misleading, it's actually the signed distance.
• If you are looking to implement this in another language, you may not have .project() available. a.project(b) is equivalent to b * (a.dot(b) * b.length_squared()), however for the uses here, b is normalized so b.length_squared() is always 1, so you can simplify this to b * a.dot(b).