Calculate the distance between two line segments accroding to https://homepage.univie.ac.at/franz.vesely/notes/hard_sticks/hst/hst.html

line segment distance

import numpy as np
from scipy.optimize import minimize

def get_in_plane_dist_sqr(x, e1, e2):
    gamma, delta = x
    dist_sqr = gamma**2 + delta**2 - 2 * gamma * delta * (e1 @ e2)
    return dist_sqr

def get_segment_dist(r1, r2, e1, e2, L1, L2):
    """
    get the minimum distance between two line segments
    r1, r2: the starting points of the two segments
    e1, e2: the direction vectors of the two segments
    L1, L2: the lengths of the segments
    """
    r1, r2, e1, e2 = list(map(np.array, [r1, r2, e1, e2]))
    e1 = e1 / np.linalg.norm(e1)
    e2 = e2 / np.linalg.norm(e2)
    c1, c2 = r1 + e1 * L1/2, r2 + e2 * L2/2  # center points
    r12 = c2 - c1
    const = 1 / (1 - (e1 @ e2)**2)
    lambda_0 = const * (r12 @ e1 - (r12 @ e2) * (e1 @ e2))
    mu_0 =    -const * (r12 @ e2 - (r12 @ e1) * (e1 @ e2))
    s = r12 + mu_0 * e2 - lambda_0 * e1
    line_dist_sqr = s @ s
    if (np.abs(lambda_0) < np.abs(L1/2)) and (np.abs(mu_0) < np.abs(L2/2)):
        return np.sqrt(line_dist_sqr)
    else:
        init_guess = [0, 0]  # gamma, delta
        res = minimize(
            fun=get_in_plane_dist_sqr,
            x0=init_guess,
            args=(e1, e2),
            method='SLSQP', bounds=[
                (-lambda_0 - L1/2, -lambda_0 + L1/2),
                (-mu_0 - L2/2, -mu_0 + L2/2)
            ]
        )
        gamma, delta = res.x
        dist_in_plane_sqr = get_in_plane_dist_sqr((gamma, delta), e1, e2)
        return np.sqrt(dist_in_plane_sqr + line_dist_sqr)
    
def get_line_dist(r1, r2, e1, e2):
    n = np.cross(e1, e2)
    d = n @ (r1 - r2) / np.linalg.norm(n)
    return np.abs(d)
    

if __name__ == "__main__":
    import numpy as np
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt

    r1, r2 = np.zeros(3), np.ones(3)
    
    for j in range(100):
        e1 = np.random.random(3)
        e2 = np.random.random(3)
        L1, L2 = np.random.uniform(0, 5), np.random.uniform(0, 5)
        
        seg_dist = get_segment_dist(r1, r2, e1, e2, L1, L2)
        line_dist = get_line_dist(r1, r2, e1, e2)
        if (seg_dist < line_dist - 1e-10) or (seg_dist > np.sqrt(3)):
            print(f'{seg_dist:.4f} > {line_dist:.4f}, {seg_dist:.4f} <= 1.732')
            fig = plt.figure()
            ax = fig.add_subplot(111, projection='3d')
            ax.scatter(*r1, color='tomato')
            ax.quiver3D(*r1, *e1/np.linalg.norm(e1) * L1, color='tomato')
            #ax.quiver3D(*r1, *-e1, L1, color='tomato')
            ax.scatter(*r2, color='teal')
            ax.quiver3D(*r2, *e2/np.linalg.norm(e2) * L2, color='teal')
            #ax.quiver3D(*r2, *-e2, L2, color='teal')
            ax.set_xlim([-2, 2])
            ax.set_ylim([-2, 2])
            ax.set_zlim([-2, 2])
            plt.show()