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

Question

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

score 24 · Accepted Answer

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

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);               
}

score 5 · Accepted Answer

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

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

您的平面方程如下所示：a*x+b*y+c*z+d=0 只需填写点值 (xyz)。如果结果的符号 >0，则该点与法线在同一侧，结果 == 0，点位于平面内，在您的情况下，您需要第三个选项：<0 表示它位于飞机。如果这对所有 4 个平面都满足，则您的点位于四面体内部。

score 4 · Accepted Answer

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

在这个坐标系中，如果 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 实现。测试表明这种方法比平面方法快 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))

score 4 · Accepted Answer

从Hugues 的解决方案开始，这里有一个更简单且（甚至）更有效的解决方案：

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<1。

作为比较，这里是比较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对过程进行矢量化的想法，如果想找出网格的哪个元素包含给定点，这里有一个高度矢量化的解决方案：

输入数据：

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 倍。这样的网格非常小，因此对于较大的网格，此间隙可能会增加。

score 3 · Accepted Answer

我已经对 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))

score 1 · Accepted Answer

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

直觉

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

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

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

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

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

detA > 0 且 detB < 0 且 detC > 0 且 detD < 0

或者

detA < 0 且 detB > 0 且 detC < 0 且 detD > 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 ---

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

6 回答 6

Related

Reference