这是我想到的代码。我在scipy中找不到任何直接帮助我的东西,而这个解决方案比查询CGAL快大约2倍。它不能处理折叠的三角形,但可以通过检查边长并返回最长边上最接近的点来修复。
import numpy as np
def pointsToTriangles(points,triangles):
with np.errstate(all='ignore'):
p0,p1,p2 = np.asarray(triangles).swapaxes(0,1)
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)
det = a*c - b*b
invDet = 1. / det
denom = a-2*b+c
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
bd = b+d
ce = c+e
numer0 = (ce - bd) / denom
numer1 = (c+e-b-d) / denom
da = -d/a
ec = -e/c
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 (p0.T + u[:, np.newaxis] * e0.T + v[:, np.newaxis] * e1.T).swapaxes(2,1)
一些测试数据将100个点投影到10k个三角形中:
import numpy as np
import cProfile
N_TRIANGLES = 10**4
N_POINTS = 10**2
points = np.random.random((N_POINTS,3,)) * 100
triangles = np.random.random((N_TRIANGLES,3,3,)) * 100
cProfile.run("pointsToTriangles(points,triangles)")
这将很快占用大量内存,因此在处理大型数据集时,最好逐个迭代点或三角形。
projections
应该包含什么内容? - ali_m