使用numpy/scipy进行三维点到三角形的最近点投影

Question

使用numpy/scipy进行三维点到三角形的最近点投影

3

如何使用numpy/scipy计算一个点到N个三角形的最近投影？

目前，我会编写一个函数来计算单个三角形的投影，基本上是这样的，然后遍历整个三角形数组。但在开始之前，我想知道是否已经在scipy中内置了解决方案。类似于：

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

这里有一张图来帮助您直观地理解：

在上面的图像中，中间的红色点是我的查询“点”，蓝色点是每个三角形的顶点，如在“三角形”np.array()中定义。绿色点代表我想要的结果。它们是“点”到定义的三角形的最近投影点，并且我希望将此信息作为点数组返回。

谢谢！

- Fnord

你能解释一下你需要什么吗？在你的例子中，projections应该包含什么内容？ - ali_m

我仍然不明白你所说的“将‘点’投影到三角形表面”的意思。如果我有一个点和一个三角形在三维空间中，那么我可以画出无数条连接点和三角形表面的线。你能否找到/画出一个图示，或者至少给出一个带有实际数字的例子？ - ali_m

1

@ali_m 抱歉，我正在寻找三角形上最接近某点的投影。 - Fnord

你可以将该点投影到每个三角形的平面上，这是你想要的吗？或者你还想找到距离该投影点最近的三角形内的点吗？ - Jaime

1

通常情况下，投影沿着一个明确定义的方向进行。如果不是三角形平面的法线，则需要定义沿哪个方向进行投影，否则您的问题将是未定义的。 - Jaime

显示剩余7条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Fnord · Accepted Answer

这是我想到的代码。我在scipy中找不到任何直接帮助我的东西，而这个解决方案比查询CGAL快大约2倍。它不能处理折叠的三角形，但可以通过检查边长并返回最长边上最接近的点来修复。

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)

一些测试数据将100个点投影到10k个三角形中：

    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

这将很快占用大量内存，因此在处理大型数据集时，最好逐个迭代点或三角形。