I just did some math work on the distance between two lines in 3D space. After some boring calculation, I found the result itself is much more decent and help us to revisit the issue in an elegant way.

Line 1: p = p1 + t * v1

Line 2: p = p2 + t * v2

(t is a scalar argument, while v1 and v2 are normalized vectors)

Before considering the special cases, we can consider the general case and usually the special cases can be analyzed from some items in the general result.

SO the general case is that the two lines are not coplanar. We can first construct a plane from these two lines by translating Line 2 to pass through p1 and calculating the distance from p2 to the constructed plane.

The normal vector of the plane is:

n = cross( v1, v2 )

Thus the plane is:

dot( n, p ) = dot( n, p1 )

The closest point in the plane to p2 is constrained by the above equation and the following one as well.

p = p2 + t * n

t in the above equation is therefore the distance we want to calculated.

By combining these two equations, we can get this one:

dot( n, p2 + t * n ) = dot( n, p1 )

i.e.

t = abs( dot( n, p1 - p2 ) ) = abs( dot( cross( v1, v2 ), p1 - p2 ) )

After this result comes out, I said to myself, "aah, how stupid I am..., I don't actually need to do so much calculation!"

It is actually a decent result if you consider this question in another way. You can pick any two points respectively from these two lines, let's say, p1 and p2 is OK. Then after projecting the vector p1 - p2 onto the normal n is the distance we want.

For the special cases, we can see from cross( v1, v2 ) that the equation is incorrect when the two lines are parallel. It will be easy to analyze it.