Simple way to calculate minimum distance beween two 3D objects

Question

Suppose I have two tetrahedras and it's vertex coordinates. They cannot have parts inside each other and their distance will always be greater than 0. How to calculate the minimum distance between them?

My approach, which I did not like very much, because I'm having problems representing the line equations in code:

• For each tetrahedra, calculate 6 lines equations, one for each edge.
• Them calculate the distance between the points of tetrahedra A to the lines equations of tetrahedra B. Store those 24 (6 edges, 4 vertexes) distance in a list.
• Do the same in the other way.
• Find the minimum between the 48 distances calculated

Is there any algorithm to find minimum distance between 3D objects, given it's coordinates? I believe the above procedure solve the problem but it is very tedious to write all the equations

Answer 1

Here is my simple approach. If we have tetrahedra points in arrays t1 and t2, we initiate minimal distance to distance between any two vertices in tetrahedron A and tetrahedron B. We loop over all faces in A and B and, and define functions face_p_1(alpha,beta) which we use to get any point on that face by scanning parameters alpha and beta in range [0,1]. And finally we minimize distance (using scipy) between all points in each pair of faces in A and B.

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
import numpy as np
import itertools
from scipy.optimize import minimize

def mindist (t1,t2):
    dmin=np.linalg.norm(t1[0]-t2[0])
    p0=t1[0]
    p1=t2[0]

    for face_t1 in itertools.combinations([0,1, 2, 3], 3):
        for face_t2 in itertools.combinations([0,1, 2, 3], 3):
            pts_f_t1=t1[list(face_t1)]
            pts_f_t2=t2[list(face_t2)]

            def face_p_1(alpha1,beta1):
                v1=pts_f_t1[1]-pts_f_t1[0]
                v2=pts_f_t1[2]-pts_f_t1[1]
                return(pts_f_t1[0]+(v1*alpha1+v2*alpha1*beta1))

            def face_p_2(alpha2,beta2):
                v1=pts_f_t2[1]-pts_f_t2[0]
                v2=pts_f_t2[2]-pts_f_t2[1]
                return(pts_f_t2[0]+(v1*alpha2+v2*alpha2*beta2))

            def fpdist(a):
                alpha1,beta1,alpha2,beta2=a
                return(np.linalg.norm(face_p_1(alpha1,beta1)-face_p_2(alpha2,beta2)))

            res=minimize(fpdist,[0.5,0.5,0.5,0.5],bounds=[(0,1),(0,1),(0,1),(0,1)])
            if( res["fun"] < dmin):
                dmin=res["fun"]
                p0=face_p_1(*(res["x"][[0,1]]))
                p1=face_p_2(*(res["x"][[2,3]]))
    return dmin,p0,p1

    
def plot_tetrahedron(v,ax):
    ax.scatter3D(v[:, 0], v[:, 1], v[:, 2])
    verts = [ [v[0],v[1],v[2]], [v[0],v[1],v[3]], [v[0],v[2],v[3]] , [v[1],v[2],v[3]]  ]
    ax.add_collection3d(Poly3DCollection(verts, 
     facecolors='cyan', linewidths=1, edgecolors='r', alpha=.25))

Here are some tests with plots

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

t1=np.array([[ 0.5,  0.5, -0.2],[ 0.6, -0.1, -0.2],[-0.2,  0.6, -0.1],[ 0. ,  0. ,  0.6]])
t2=np.array([[0.6 , 0.65, 0.55],[1.4 , 0.75, 0.55],[0.6 , 1.45, 0.65],[0.8 , 0.85, 1.35]])

dmin,p1,p2=mindist (t1,t2)

plot_tetrahedron(t1,ax)
plot_tetrahedron(t2,ax)

ax.plot(*(array([ p1 , p2 ]).T), c='k')

another case

t1=np.array([[-0.15 , -0.225, -0.075],[-0.05 , -0.825, -0.775],[-0.85 , -0.125, -0.675],[-0.65 , -0.725,  0.025]])
t2=np.array([[0.6 , 0.65, 0.55],[1.4 , 0.75, 0.55],[0.6 , 1.45, 0.65],[0.8 , 0.85, 1.35]])


t1=np.array([[-0.15 , -0.225, -0.075],[-0.05 , -0.825, -0.775],[-0.85 , -0.125, -0.675],[-0.65 , -0.725,  0.025]])
t2=np.array([[0.6 , 0.65, 0.55],[1.4 , 0.75, 0.55],[0.6 , 1.45, 0.65],[0.8 , 0.85, 1.35]])

t1=np.array([[0.1, 0.1, 0.1],[0.9, 0.2, 0.1],[0.1, 0.9, 0.2],[0.3, 0.3, 0.9]])
t2=np.array([[0.6 , 0.65, 0.55],[1.4 , 0.75, 0.55],[0.6 , 1.45, 0.65],[0.8 , 0.85, 1.35]])

and case with two parallel faces

t1=np.array([[-1 , 0, 0],[1 , 0, 0],[0 , sqrt(2), 0],[0 , sqrt(2)/3, sqrt(2)]])
t2=np.array([[-1 , 0, -1],[1 , 0, -1],[0 , sqrt(2), -1],[0 , sqrt(2)/3, -sqrt(2)-1]])