Fnord - 1 year ago 202
Python Question

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

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)
``````

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!

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
from numpy.core.umath_tests import inner1d

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 = inner1d(e0,e0)
b = inner1d(e0,e1)
c = inner1d(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 = inner1d(e0,p)
e = inner1d(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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download