Fnord - 6 months ago 37

Python 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

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.