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 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.