Closest point projection of a 3D point to 3D triangles with numpy/scipy

Question

How would one calculate the closest projections of a point to N triangles using numpy/scipy?

Right now I would make a function to calculate a projection to a single triangle, basically this, then iterate across the entire array of triangles. But before I set off to do this I was wondering if there's already a solution built in scipy. Something like:

# DREAMY PSEUDOCODE
import numpy as np
N_TRIANGLES = 1000

point = np.random.rand(3) * 100 #random 3d point
triangles = np.random.rand(N_TRIANGLES,3,3) * 100 #array of triangles

from scipy.spatial import pointToTriangles
projections = pointToTriangles(point,triangles)

Here's a picture to help you visualize:

In the image above, the red dot in the middle is my query "point", the blue dots are the vertices of each triangles as define in the "triangles" np.array(). The green dots represent the result i want. They're the closest projections of "point" onto the defined triangles, and i wish to get this info back as an array of points.

cheers!

Answer 1

Here's the code i came up with. I couldn't find anything in scipy that could help me directly, and this solution is about 2x faster than querying CGAL. It doesn't handle collapsed triangles, but that could be fixed by checking edge lengths and return the closest point on the longest edge instead.

import numpy as np

def pointsToTriangles(points,triangles):
    
    with np.errstate(all='ignore'):
        
        # Unpack triangle points
        p0,p1,p2 = np.asarray(triangles).swapaxes(0,1)
                
        # Calculate triangle edges
        e0 = p1-p0
        e1 = p2-p0
        a = np.einsum('...i,...i', e0, e0)
        b = np.einsum('...i,...i', e0, e1)
        c = np.einsum('...i,...i', e1, e1)
        
        # Calculate determinant and denominator
        det = a*c - b*b
        invDet = 1. / det
        denom = a-2*b+c
        
        # Project to the edges
        p  = p0-points[:,np.newaxis]
        d = np.einsum('...i,...i', e0, p)
        e = np.einsum('...i,...i', e1, p)
        u = b*e - c*d
        v = b*d - a*e
        
        # Calculate numerators
        bd = b+d
        ce = c+e
        numer0 = (ce - bd) / denom
        numer1 = (c+e-b-d) / denom
        da = -d/a
        ec = -e/c
        
        
        # Vectorize test conditions
        m0 = u + v < det
        m1 = u < 0
        m2 = v < 0
        m3 = d < 0
        m4 = (a+d > b+e)
        m5 = ce > bd
        
        t0 =  m0 &  m1 &  m2 &  m3
        t1 =  m0 &  m1 &  m2 & ~m3
        t2 =  m0 &  m1 & ~m2
        t3 =  m0 & ~m1 &  m2
        t4 =  m0 & ~m1 & ~m2
        t5 = ~m0 &  m1 &  m5
        t6 = ~m0 &  m1 & ~m5
        t7 = ~m0 &  m2 &  m4
        t8 = ~m0 &  m2 & ~m4
        t9 = ~m0 & ~m1 & ~m2
        
        u = np.where(t0, np.clip(da, 0, 1), u)
        v = np.where(t0, 0, v)
        u = np.where(t1, 0, u)
        v = np.where(t1, 0, v)
        u = np.where(t2, 0, u)
        v = np.where(t2, np.clip(ec, 0, 1), v)
        u = np.where(t3, np.clip(da, 0, 1), u)
        v = np.where(t3, 0, v)
        u *= np.where(t4, invDet, 1)
        v *= np.where(t4, invDet, 1)
        u = np.where(t5, np.clip(numer0, 0, 1), u)
        v = np.where(t5, 1 - u, v)
        u = np.where(t6, 0, u)
        v = np.where(t6, 1, v)
        u = np.where(t7, np.clip(numer1, 0, 1), u)
        v = np.where(t7, 1-u, v)
        u = np.where(t8, 1, u)
        v = np.where(t8, 0, v)
        u = np.where(t9, np.clip(numer1, 0, 1), u)
        v = np.where(t9, 1-u, v)
        
        
        # Return closest points
        return (p0.T +  u[:, np.newaxis] * e0.T + v[:, np.newaxis] * e1.T).swapaxes(2,1)

    

Some test data projecting 100 points to 10k triangles:

    import numpy as np
    import cProfile
    
    N_TRIANGLES = 10**4 # 10k triangles
    N_POINTS    = 10**2 # 100 points
    points      = np.random.random((N_POINTS,3,)) * 100
    triangles   = np.random.random((N_TRIANGLES,3,3,)) * 100
    
    cProfile.run("pointsToTriangles(points,triangles)") # 54 function calls in 0.320 seconds

This will quickly become a memory hog, so when dealing with large data sets it's probably better to iterate over either points or triangles one at a time.