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.
 A: 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).$
A: 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);


A: 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}
A: 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$.
A: 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.
A: 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.]"

A: 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).

