Introduction

Finding the shortest distance between two lines in 3D is a standard geometry problem, but the segment version is the one that usually matters in code. Infinite lines extend forever, while segments stop at endpoints, so a formula that works for lines is not enough for collision tests, CAD tools, or simulation code.

The reliable way to think about the problem is to find the two closest points, one on each object, and then measure the distance between them. For infinite lines, those points lie on a common perpendicular when the lines are not parallel. For segments, the closest points may instead sit on one or both endpoints.

Parametric Setup

Write the two objects in parametric form:

$$P(s) = P_0 + s \mathbf{u}$$

$$Q(t) = Q_0 + t \mathbf{v}$$

Here, $P_0$ and $Q_0$ are starting points, and $\mathbf{u}$ and $\mathbf{v}$ are direction vectors. For infinite lines, $s$ and $t$ can be any real values. For line segments, they are restricted to the interval $[0, 1]$.

The squared distance between points on the two objects is:

$$f(s, t) = \|P(s) - Q(t)\|^2$$

Using the squared distance keeps the algebra simpler and avoids an unnecessary square root while solving the minimization problem.

Shortest Distance for Infinite Lines

Let:

$$\mathbf{w} = P_0 - Q_0$$

and define the dot products:

$$a = \mathbf{u} \cdot \mathbf{u},\quad b = \mathbf{u} \cdot \mathbf{v},\quad c = \mathbf{v} \cdot \mathbf{v},\quad d = \mathbf{u} \cdot \mathbf{w},\quad e = \mathbf{v} \cdot \mathbf{w}$$

For non-parallel lines, the minimizing parameters are:

$$s = \frac{be - cd}{ac - b^2}, \quad t = \frac{ae - bd}{ac - b^2}$$

The denominator $ac - b^2$ becomes zero only when the direction vectors are parallel or numerically very close to parallel.

There is also a compact distance-only formula. If $\mathbf{n} = \mathbf{u} \times \mathbf{v}$, then the distance between two skew lines is:

$$\text{distance} = \frac{|(Q_0 - P_0) \cdot \mathbf{n}|}{\|\mathbf{n}\|}$$

That expression is useful for quick reasoning, but many implementations also need the closest points, so the parameter-based form is usually the better starting point.

Why Segments Need Extra Work

For segments, the unconstrained solution for $s$ and $t$ may fall outside $[0, 1]$. When that happens, the closest point on a segment is not an interior point anymore; it is an endpoint.

So the segment problem is a constrained minimization problem:

$$0 \le s \le 1,\quad 0 \le t \le 1$$

A practical routine usually follows this logic:

1. Solve for the closest points on the infinite lines.
2. Clamp the parameters into $[0, 1]$ if they fall outside the segment bounds.
3. Recompute the paired parameter after clamping, because fixing one endpoint changes the best point on the other segment.

That last step is easy to miss. Clamping one parameter and leaving the other untouched can produce the wrong answer.

Runnable Python Example

The following function returns the shortest distance and the actual closest points between two 3D line segments:

1from math import sqrt
2
3
4def sub(a, b):
5    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])
6
7
8def add(a, b):
9    return (a[0] + b[0], a[1] + b[1], a[2] + b[2])
10
11
12def mul(a, s):
13    return (a[0] * s, a[1] * s, a[2] * s)
14
15
16def dot(a, b):
17    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
18
19
20def clamp(x, lo, hi):
21    return max(lo, min(hi, x))
22
23
24def norm(a):
25    return sqrt(dot(a, a))
26
27
28def segment_distance(p0, p1, q0, q1):
29    u = sub(p1, p0)
30    v = sub(q1, q0)
31    w = sub(p0, q0)
32
33    a = dot(u, u)
34    b = dot(u, v)
35    c = dot(v, v)
36    d = dot(u, w)
37    e = dot(v, w)
38
39    eps = 1e-12
40    denom = a * c - b * b
41
42    if denom < eps:
43        s = 0.0
44        t = clamp(e / c if c > eps else 0.0, 0.0, 1.0)
45    else:
46        s = clamp((b * e - c * d) / denom, 0.0, 1.0)
47        t = (a * e - b * d) / denom
48
49        if t < 0.0:
50            t = 0.0
51            s = clamp(-d / a if a > eps else 0.0, 0.0, 1.0)
52        elif t > 1.0:
53            t = 1.0
54            s = clamp((b - d) / a if a > eps else 0.0, 0.0, 1.0)
55
56    cp = add(p0, mul(u, s))
57    cq = add(q0, mul(v, t))
58    return norm(sub(cp, cq)), cp, cq
59
60
61distance, cp, cq = segment_distance(
62    (1.0, 0.0, 0.0),
63    (1.0, 2.0, 0.0),
64    (0.0, 1.0, 1.0),
65    (2.0, 1.0, 1.0),
66)
67
68print(distance)
69print(cp)
70print(cq)

For this example, the distance is 1.0. The closest points are (1.0, 1.0, 0.0) and (1.0, 1.0, 1.0).

Interpreting Special Cases

Three cases matter.

First, intersecting lines or segments have distance 0. In floating-point code, compare against a tolerance instead of exact equality.

Second, parallel objects make the denominator very small. That does not mean the algorithm failed; it means there is no unique common perpendicular direction, so the implementation has to fall back to endpoint and projection logic.

Third, degenerate segments can happen when the two endpoints are identical. In that case, one segment is just a point, and the problem reduces to point-to-segment distance.

Common Pitfalls

• Using the infinite-line formula directly for segments and forgetting to clamp the parameters.
• Comparing floating-point values to exact zero instead of using a tolerance for parallel or intersecting cases.
• Returning only the distance when downstream code also needs the closest points for visualization or contact resolution.
• Forgetting that a zero-length segment is a point and needs special handling.
• Computing the square root too early instead of minimizing squared distance first.

Summary

• Model each line or segment with a parametric equation.
• For infinite lines, solve for the minimizing parameters with dot products.
• For skew lines, the cross-product formula gives a compact distance check.
• For segments, clamp parameters into $[0, 1]$ and recompute when an endpoint becomes active.
• A robust implementation should handle parallel lines, intersections, and zero-length segments.