如何检查点是否在四面体中？

【问题标题】：How to check whether the point is in the tetrahedron or not?如何检查点是否在四面体中？
【发布时间】：2014-09-30 12:40:42
【问题描述】：

我知道四面体的所有坐标以及我想确定的点。那么有人知道该怎么做吗？我试图确定该点属于四面体的每个三角形，如果它对所有三角形都是正确的，那么该点就在四面体中。但这绝对是错误的。

【问题讨论】：

    标签： c# opengl tetrahedra


    【解决方案1】：

    对于四面体的每个平面，检查该点是否与剩余顶点在同一侧：

    bool SameSide(v1, v2, v3, v4, p)
{
    normal := cross(v2 - v1, v3 - v1)
    dotV4 := dot(normal, v4 - v1)
    dotP := dot(normal, p - v1)
    return Math.Sign(dotV4) == Math.Sign(dotP);
}

    你需要为每架飞机检查这个：

    bool PointInTetrahedron(v1, v2, v3, v4, p)
{
    return SameSide(v1, v2, v3, v4, p) &&
           SameSide(v2, v3, v4, v1, p) &&
           SameSide(v3, v4, v1, v2, p) &&
           SameSide(v4, v1, v2, v3, p);               
}

    【讨论】：

      【解决方案2】：

      您通过四个顶点 A B C 和 D 定义一个四面体。 因此，您也可以使用 4 个三角形来定义四面体的表面。

      您现在只需检查点 P 是否在平面的另一侧。每个平面的法线指向远离四面体的中心。 所以你只需要针对 4 架飞机进行测试。

      您的平面方程如下所示：a*x+b*y+c*z+d=0 只需填写点值 (x y z)。如果结果的符号 >0，则该点与法线位于同一侧，结果 == 0，点位于平面内，在您的情况下，您需要第三个选项：

      【讨论】：

      • 如果我猜对了，首先我必须找出每个平面的 A、B、C、D 并计算 ax+by+c*z +d。如果每个平面的结果
      • 不，命名是巧合。 A B C D 是四面体的四个点，与方程中的 a b c d 不同
      • 是的，我的意思是方程中的 a,b,c,d 而不是 A,B,C,D。对不起
      • 这个答案假设您可以知道法线的方向，而不是一般情况。
      • 如果我有三角形的 3 个点，我可以使用叉积计算法线。
      【解决方案3】：

      给定定义非退化四面体的 4 个点 A、B、C、D 和要测试的点 P，一种方法是将 P 的坐标转换为四面体坐标系，例如以 A 为原点, 向量 BA, CA, DA 为单位向量。

      在这个坐标系中，如果 P 在 P 内部，P 的坐标都在 0 和 1 之间，但它也可以在由原点和 3 个单位向量定义的变换后的立方体中的任何位置。 断言 P 在 (A,B,C,D) 内部的一种方法是依次将点 (A、B、C 和 D) 和其他三个点作为原点来定义一个新的坐标系。此测试重复 4 次有效，但可以改进。

      最有效的方法是只变换一次坐标，复用前面提出的SameSide函数，例如以A为原点，变换为(A,B,C,D)坐标系，P和A必须位于在 (B,C,D) 平面的同一侧。

      以下是该测试的 numpy/python 实现。测试表明，这种方法比 Planes 方法快 2-3 倍。

      import numpy as np

def sameside(v1,v2,v3,v4,p):
    normal = np.cross(v2-v1, v3-v1)
    return ((np.dot(normal, v4-v1)*p.dot(normal, p-v1) > 0)

def tetraCoord(A,B,C,D):
    v1 = B-A ; v2 = C-A ; v3 = D-A
    # mat defines an affine transform from the tetrahedron to the orthogonal system
    mat = np.concatenate((np.array((v1,v2,v3,A)).T, np.array([[0,0,0,1]])))
    # The inverse matrix does the opposite (from orthogonal to tetrahedron)
    M1 = np.linalg.inv(mat)
    return(M1)

def pointInsideT(v1,v2,v3,v4,p):
    # Find the transform matrix from orthogonal to tetrahedron system
    M1=tetraCoord(v1,v2,v3,v4)
    # apply the transform to P
    p1 = np.append(p,1)
    newp = M1.dot(p1)
    # perform test
    return(np.all(newp>=0) and np.all(newp <=1) and sameside(v2,v3,v4,v1,p))

      【讨论】：

        【解决方案4】：

        从Hugues' solution开始，这里有一个更简单、（甚至）更高效的：

        import numpy as np

def tetraCoord(A,B,C,D):
  # Almost the same as Hugues' function, 
  # except it does not involve the homogeneous coordinates.
  v1 = B-A ; v2 = C-A ; v3 = D-A
  mat = np.array((v1,v2,v3)).T
  # mat is 3x3 here
  M1 = np.linalg.inv(mat)
  return(M1)

def pointInside(v1,v2,v3,v4,p):
  # Find the transform matrix from orthogonal to tetrahedron system
  M1=tetraCoord(v1,v2,v3,v4)
  # apply the transform to P (v1 is the origin)
  newp = M1.dot(p-v1)
  # perform test
  return (np.all(newp>=0) and np.all(newp <=1) and np.sum(newp)<=1)

        在与四面体相关的坐标系中，与原点相对的面（此处表示为 v1）的特征是 x+y+z=1。因此，该面与v1同侧的半空间满足x+y+z

        作为比较，这里是比较Nico、Hugues和我提出的方法的完整代码：

        import numpy as np
import time

def sameside(v1,v2,v3,v4,p):
    normal = np.cross(v2-v1, v3-v1)
    return (np.dot(normal, v4-v1) * np.dot(normal, p-v1) > 0)

# Nico's solution
def pointInside_Nico(v1,v2,v3,v4,p):   
    return sameside(v1, v2, v3, v4, p) and sameside(v2, v3, v4, v1, p) and sameside(v3, v4, v1, v2, p) and sameside(v4, v1, v2, v3, p)      

# Hugues' solution
def tetraCoord(A,B,C,D):
    v1 = B-A ; v2 = C-A ; v3 = D-A
    # mat defines an affine transform from the tetrahedron to the orthogonal system
    mat = np.concatenate((np.array((v1,v2,v3,A)).T, np.array([[0,0,0,1]])))
    # The inverse matrix does the opposite (from orthogonal to tetrahedron)
    M1 = np.linalg.inv(mat)
    return(M1)
    
def pointInside_Hugues(v1,v2,v3,v4,p):
    # Find the transform matrix from orthogonal to tetrahedron system
    M1=tetraCoord(v1,v2,v3,v4)
    # apply the transform to P
    p1 = np.append(p,1)
    newp = M1.dot(p1)
    # perform test
    return(np.all(newp>=0) and np.all(newp <=1) and sameside(v2,v3,v4,v1,p))    

# Proposed solution
def tetraCoord_Dorian(A,B,C,D):
    v1 = B-A ; v2 = C-A ; v3 = D-A
    # mat defines an affine transform from the tetrahedron to the orthogonal system
    mat = np.array((v1,v2,v3)).T
    # The inverse matrix does the opposite (from orthogonal to tetrahedron)
    M1 = np.linalg.inv(mat)
    return(M1) 

def pointInside_Dorian(v1,v2,v3,v4,p):
    # Find the transform matrix from orthogonal to tetrahedron system
    M1=tetraCoord_Dorian(v1,v2,v3,v4)
    # apply the transform to P
    newp = M1.dot(p-v1)
    # perform test
    return (np.all(newp>=0) and np.all(newp <=1) and np.sum(newp)<=1)
    
npt=100000    
Pt=np.random.rand(npt,3)
A=np.array([0.1, 0.1, 0.1])
B=np.array([0.9, 0.2, 0.1])
C=np.array([0.1, 0.9, 0.2])
D=np.array([0.3, 0.3, 0.9])

inTet_Nico=np.zeros(shape=(npt,1),dtype=bool)
inTet_Hugues=inTet_Nico
inTet_Dorian=inTet_Nico

start_time = time.time()
for i in range(0,npt):
    inTet_Nico[i]=pointInside_Nico(A,B,C,D,Pt[i,:])
print("--- %s seconds ---" % (time.time() - start_time)) # https://stackoverflow.com/questions/1557571/how-do-i-get-time-of-a-python-programs-execution

start_time = time.time()
for i in range(0,npt):
    inTet_Hugues[i]=pointInside_Hugues(A,B,C,D,Pt[i,:])
print("--- %s seconds ---" % (time.time() - start_time))   

start_time = time.time()    
for i in range(0,npt):
    inTet_Dorian[i]=pointInside_Dorian(A,B,C,D,Pt[i,:])
print("--- %s seconds ---" % (time.time() - start_time))

        以下是运行时间方面的结果：

        --- 15.621951341629028 seconds ---
--- 8.97989797592163 seconds ---
--- 4.597853660583496 seconds ---

        [编辑]

        基于矢量化过程的Tom's idea，如果想要找到哪个网格元素包含给定点，这里有一个高度矢量化的解决方案：

        输入数据：

        • node_coordinates: (n_nodes,3) 包含每个节点坐标的数组
        • node_ids: (n_tet, 4) 数组，其中第 i 行给出了第 i 个四面体的顶点索引。
        def where(node_coordinates, node_ids, p):
    ori=node_coordinates[node_ids[:,0],:]
    v1=node_coordinates[node_ids[:,1],:]-ori
    v2=node_coordinates[node_ids[:,2],:]-ori
    v3=node_coordinates[node_ids[:,3],:]-ori
    n_tet=len(node_ids)
    v1r=v1.T.reshape((3,1,n_tet))
    v2r=v2.T.reshape((3,1,n_tet))
    v3r=v3.T.reshape((3,1,n_tet))
    mat = np.concatenate((v1r,v2r,v3r), axis=1)
    inv_mat = np.linalg.inv(mat.T).T    # https://stackoverflow.com/a/41851137/12056867        
    if p.size==3:
        p=p.reshape((1,3))
    n_p=p.shape[0]
    orir=np.repeat(ori[:,:,np.newaxis], n_p, axis=2)
    newp=np.einsum('imk,kmj->kij',inv_mat,p.T-orir)
    val=np.all(newp>=0, axis=1) & np.all(newp <=1, axis=1) & (np.sum(newp, axis=1)<=1)
    id_tet, id_p = np.nonzero(val)
    res = -np.ones(n_p, dtype=id_tet.dtype) # Sentinel value
    res[id_p]=id_tet
    return res

        这里的技巧是用多维数组做矩阵乘积。

        where 函数将点坐标作为第三个参数。事实上，这个函数可以同时在多个坐标上运行；输出参数的长度与p 相同。如果对应的坐标不在网格中，则返回-1。

        在由 1235 个四面体组成的网格上，这种方法比在每个四面体上循环要快 170-180 倍。这样的网格非常小，因此对于较大的网格，此间隙可能会增加。

        【讨论】：

          【解决方案5】：

          我对 Dorian 和 Hughes 解决方案进行了矢量化处理，以将整个点数组作为输入。我还将 tetraCoord 函数移到 pointsInside 函数之外并重命名了两者，因为没有必要为每个点调用它。

          在我的电脑上，@Dorian 的解决方案和示例在 2.5 秒内运行。在相同的数据上，我的运行速度快了近一千倍，为 0.003 秒。如果出于某种原因需要更快的速度，则将 GPU cupy 包导入为“np”会将其推入 100 微秒范围。

          import time
# alternatively, import cupy as np if len(points)>1e7 and GPU
import numpy as np 

def Tetrahedron(vertices):
    """
    Given a list of the xyz coordinates of the vertices of a tetrahedron, 
    return tetrahedron coordinate system
    """
    origin, *rest = vertices
    mat = (np.array(rest) - origin).T
    tetra = np.linalg.inv(mat)
    return tetra, origin

def pointInside(point, tetra, origin):
    """
    Takes a single point or array of points, as well as tetra and origin objects returned by 
    the Tetrahedron function.
    Returns a boolean or boolean array indicating whether the point is inside the tetrahedron.
    """
    newp = np.matmul(tetra, (point-origin).T).T
    return np.all(newp>=0, axis=-1) & np.all(newp <=1, axis=-1) & (np.sum(newp, axis=-1) <=1)

npt=10000000
points = np.random.rand(npt,3)
# Coordinates of vertices A, B, C and D
A=np.array([0.1, 0.1, 0.1])
B=np.array([0.9, 0.2, 0.1])
C=np.array([0.1, 0.9, 0.2])
D=np.array([0.3, 0.3, 0.9])

start_time = time.time()
vertices = [A, B, C, D]
tetra, origin = Tetrahedron(vertices)
inTet = pointInside(points, tetra, origin)
print("--- %s seconds ---" % (time.time() - start_time))

          【讨论】：

            【解决方案6】：

            感谢 Dorian 的测试用例脚本，我可以研究另一种解决方案，并快速将其与目前的解决方案进行比较。

            直觉

            对于三角形 ABC 和点 P，如果将 P 连接到角以得到向量 PA、PB、PC 并比较由 PA、PC 和 PB、PC 跨越的两个三角形 X 和 Y，则该点如果 X 和 Y 重叠，则 P 位于三角形 ABC 内。

            或者换句话说，如果P在三角形ABC中，则不可能通过仅将正系数线性组合PC和PB来构造向量PA。

            从那里我尝试将其转移到四面体案例并阅读here，可以通过检查由向量构成的矩阵的行列式是否为非零来检查向量是否线性独立. 我尝试了使用行列式的各种方法，我偶然发现了这个：

            令 PA、PB、PC、PD 为 P 与四面体点 ABCD 的连接（即 PA = A - P 等）。计算行列式 detA = det(PB PC PD)、detB、detC 和 detD（如 detA）。

            那么点 P 位于 ABCD 所跨越的四面体内，如果：

            detA > 0 和 detB 0 和 detD

            或

            detA 0 和 detC 0

            所以行列式切换符号，从负数开始，或从正数开始。

            有效吗？显然。为什么它有效？我不知道，或者至少，我无法证明。也许其他具有更好数学技能的人可以在这里帮助我们。

            （编辑：实际上可以使用这些行列式定义重心坐标，最后，重心坐标需要加起来为一。这就像比较由 P 与点的组合所跨越的四面体的体积A,B,C,D 与四面体 ABCD 本身的体积。观察行列式符号的解释案例仍然不清楚它是否普遍有效，我不推荐它）

            我将测试用例更改为不针对一个四面体 T 检查 n 个点 Pi，而是针对 n 个四面体 Ti 检查 n 个点 Pi。所有答案仍然给出正确的结果。我认为这种方法更快的原因是它不需要矩阵求逆。 我离开了 TomNorway 用一个四面体实现的方法，我将这种新方法的矢量化留给其他人，因为我对 python 和 numpy 不太熟悉。

            import numpy as np
import time

def sameside(v1,v2,v3,v4,p):
    normal = np.cross(v2-v1, v3-v1)
    return (np.dot(normal, v4-v1) * np.dot(normal, p-v1) > 0)

# Nico's solution
def pointInside_Nico(v1,v2,v3,v4,p):   
    return sameside(v1, v2, v3, v4, p) and sameside(v2, v3, v4, v1, p) and sameside(v3, v4, v1, v2, p) and sameside(v4, v1, v2, v3, p)      

# Hugues' solution
def tetraCoord(A,B,C,D):
    v1 = B-A ; v2 = C-A ; v3 = D-A
    # mat defines an affine transform from the tetrahedron to the orthogonal system
    mat = np.concatenate((np.array((v1,v2,v3,A)).T, np.array([[0,0,0,1]])))
    # The inverse matrix does the opposite (from orthogonal to tetrahedron)
    M1 = np.linalg.inv(mat)
    return(M1)

def pointInside_Hugues(v1,v2,v3,v4,p):
    # Find the transform matrix from orthogonal to tetrahedron system
    M1=tetraCoord(v1,v2,v3,v4)
    # apply the transform to P
    p1 = np.append(p,1)
    newp = M1.dot(p1)
    # perform test
    return(np.all(newp>=0) and np.all(newp <=1) and sameside(v2,v3,v4,v1,p))    

#Dorian's solution
def tetraCoord_Dorian(A,B,C,D):
    v1 = B-A ; v2 = C-A ; v3 = D-A
    # mat defines an affine transform from the tetrahedron to the orthogonal system
    mat = np.array((v1,v2,v3)).T
    # The inverse matrix does the opposite (from orthogonal to tetrahedron)
    M1 = np.linalg.inv(mat)
    return(M1) 

def pointInside_Dorian(v1,v2,v3,v4,p):
    # Find the transform matrix from orthogonal to tetrahedron system
    M1=tetraCoord_Dorian(v1,v2,v3,v4)
    # apply the transform to P
    newp = M1.dot(p-v1)
    # perform test
    return (np.all(newp>=0) and np.all(newp <=1) and np.sum(newp)<=1)

#TomNorway's solution adapted to cope with n tetrahedrons

def Tetrahedron(vertices):
    """
    Given a list of the xyz coordinates of the vertices of a tetrahedron, 
    return tetrahedron coordinate system
    """
    origin, *rest = vertices
    mat = (np.array(rest) - origin).T
    tetra = np.linalg.inv(mat)
    return tetra, origin

def pointInside(point, tetra, origin):
    """
    Takes a single point or array of points, as well as tetra and origin objects returned by 
    the Tetrahedron function.
    Returns a boolean or boolean array indicating whether the point is inside the tetrahedron.
    """
    newp = np.matmul(tetra, (point-origin).T).T
    return np.all(newp>=0, axis=-1) & np.all(newp <=1, axis=-1) & (np.sum(newp, axis=-1) <=1)



# Proposed solution
def det3x3_Philipp(b,c,d):
    return b[0]*c[1]*d[2] + c[0]*d[1]*b[2] + d[0]*b[1]*c[2] - d[0]*c[1]*b[2] - c[0]*b[1]*d[2] - b[0]*d[1]*c[2]

def pointInside_Philipp(v0,v1,v2,v3,p):
    a = v0 - p
    b = v1 - p
    c = v2 - p
    d = v3 - p
    detA = det3x3_Philipp(b,c,d)
    detB = det3x3_Philipp(a,c,d)
    detC = det3x3_Philipp(a,b,d)
    detD = det3x3_Philipp(a,b,c)
    ret0 = detA > 0.0 and detB < 0.0 and detC > 0.0 and detD < 0.0
    ret1 = detA < 0.0 and detB > 0.0 and detC < 0.0 and detD > 0.0
    return ret0 or ret1


npt=100000
Pt= np.array([ np.array([p[0]-0.5,p[1]-0.5,p[2]-0.5]) for p in np.random.rand(npt,3)])
A=np.array([ np.array([p[0]-0.5,p[1]-0.5,p[2]-0.5]) for p in np.random.rand(npt,3)])
B=np.array([ np.array([p[0]-0.5,p[1]-0.5,p[2]-0.5]) for p in np.random.rand(npt,3)])
C=np.array([ np.array([p[0]-0.5,p[1]-0.5,p[2]-0.5]) for p in np.random.rand(npt,3)])
D=np.array([ np.array([p[0]-0.5,p[1]-0.5,p[2]-0.5]) for p in np.random.rand(npt,3)])

inTet_Nico=np.zeros(shape=(npt,1),dtype=bool)
inTet_Hugues=np.copy(inTet_Nico)
inTet_Dorian=np.copy(inTet_Nico)
inTet_Philipp=np.copy(inTet_Nico)



print("non vectorized, n points, different tetrahedrons:")

start_time = time.time()
for i in range(0,npt):
    inTet_Nico[i]=pointInside_Nico(A[i,:],B[i,:],C[i,:],D[i,:],Pt[i,:])
print("Nico's:   --- %s seconds ---" % (time.time() - start_time)) # https://stackoverflow.com/questions/1557571/how-do-i-get-time-of-a-python-programs-execution

start_time = time.time()
for i in range(0,npt):
    inTet_Hugues[i]=pointInside_Hugues(A[i,:],B[i,:],C[i,:],D[i,:],Pt[i,:])
print("Hugues':  --- %s seconds ---" % (time.time() - start_time))   

start_time = time.time()    
for i in range(0,npt):
    inTet_Dorian[i]=pointInside_Dorian(A[i,:],B[i,:],C[i,:],D[i,:],Pt[i,:])
print("Dorian's: --- %s seconds ---" % (time.time() - start_time))

start_time = time.time()
for i in range(0,npt):
    inTet_Philipp[i]=pointInside_Philipp(A[i,:],B[i,:],C[i,:],D[i,:],Pt[i,:])
print("Philipp's:--- %s seconds ---" % (time.time() - start_time))   

print("vectorized, n points, 1 tetrahedron:")

start_time = time.time()
vertices = [A[0], B[0], C[0], D[0]]
tetra, origin = Tetrahedron(vertices)
inTet_Tom = pointInside(Pt, tetra, origin)
print("TomNorway's: --- %s seconds ---" % (time.time() - start_time)) 


for i in range(0,npt):
    assert inTet_Hugues[i] == inTet_Nico[i]
    assert inTet_Dorian[i] == inTet_Hugues[i]
    #assert inTet_Tom[i] == inTet_Dorian[i] can not compare because Tom implements 1 tetra instead of n
    assert inTet_Philipp[i] == inTet_Dorian[i]

'''errors = 0
for i in range(0,npt):
    if ( inTet_Philipp[i] != inTet_Dorian[i]):
        errors = errors + 1 
print("errors " + str(errors))'''

            结果：

            non vectorized, n points, different tetrahedrons:
Nico's:   --- 25.439453125 seconds ---
Hugues':  --- 28.724457263946533 seconds ---
Dorian's: --- 15.006574153900146 seconds ---
Philipp's:--- 4.389788389205933 seconds ---
vectorized, n points, 1 tetrahedron:
TomNorway's: --- 0.008165121078491211 seconds ---

            【讨论】：

              猜你喜欢
              • 1970-01-01
              • 2017-11-16
              • 2022-01-23
              • 2012-08-22
              • 2014-12-02
              • 2022-06-24
              • 1970-01-01
              • 1970-01-01
              • 2014-01-20
              相关资源
              最近更新 更多
              热门标签
              Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode