# 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}\cdot|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.

Edit: The distance between the lines will be $$t_3 \cdot V_3$$ not $$d=t_{3}/|V_{3}|$$ as originally written.

• Is it not $d = t_3 |V_3|$ rather? Commented 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.sstatic.net/NNyA9.png) makes no sense.
– Mike
Commented Oct 21, 2019 at 11:27
• @Eskapp Mike Corrected. Commented May 27, 2023 at 7:02

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? Commented Jun 15, 2022 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. Commented Jan 3, 2022 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? Commented 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. Commented Jan 1, 2022 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. Commented Feb 26, 2022 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] Commented Feb 26, 2022 at 20:38
• What is $\vec{a_0}$? Commented Jun 15, 2022 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}

Okay, my answer from the book Mathematical structures for computer graphics by Steven I. Janke (ex.3.13). Practically, it uses Gaussian elimination to remove t2/t1 unknowns from system. So the reasoning is very simple and requires only 2 point and 2 direction vectors of lines

1. As we wish to find closest/intersection point/s we set two point descriptions equal to each other

$$P = P_1 + t_1\vec{v}_1$$

$$P = P_2 + t_2\vec{v}_2$$

$$P_1 + t_1\vec{v}_1 = P_2 + t_2\vec{v}_2 \Rightarrow t_1\vec{v}_1 - t_2\vec{v}_2 = P_2 - P_1$$

1. Then we have to eliminate $$t_2$$, as we don't like the possibility to solve system of equations on our no-solver language systems. We see that $$t_2$$ is a scalar multiplying the vector $$\vec{v}_2$$, we can cross product both sides of equation with vector $$\vec{v}_2$$ to eliminate $$t_2\vec{v}_2$$ completely (since the cross product of a vector with itself is zero).

$$(t_1 \vec{v}_1 - t_2 \vec{v}_2) \times \vec{v}_2 = t_1 (\vec{v}_1 \times \vec{v}_2) = (P_2 - P_1) \times \vec{v}_2$$

1. To turn the vector equation into an equation with real numbers, we take the dot product of each side with the vector $$\vec{v}_1 \times \vec{v}_2$$

$$t_1(\vec{v}_1 \times \vec{v}_2) \cdot (\vec{v}_1 \times \vec{v}_2) = ((P_2 - P_1) \times \vec{v}_2) \cdot (\vec{v}_1 \times \vec{v}_2)$$

1. And finally we solve for $$t_1$$

$$t_1 = \frac{((P_2 - P_1) \times \vec{v}_2)) \cdot (\vec{v}_1 \times \vec{v}_2)}{|\vec{v}_1 \times \vec{v}_2|^2}$$

If in section 2 we cross product both sides with $$v_1$$, we would get $$t_2$$ at the end by the same logic. After solving for $$t_1$$ and $$t_2$$ you can easily find intersection / closest points on both lines.

• I like this, and this can be easily implemented. The calculation can be simplified more by using vector triple product identity,
– xue
Commented Mar 30 at 0:02

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