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

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

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

用四个顶点a、B、C、D来定义一个四面体。因此你也可以用4个三角形来定义四面体的表面。

现在检查点p是否在平面的另一边。每个平面的法线都指向远离四面体中心的方向。所以你只需要测试4个平面。

你的平面方程看起来是这样的:a*x+b*y+c*z+d=0只要填入点值(x y z)。如果结果的符号>0，则点与法线在同一侧，result == 0，点位于平面内，在你的情况下，你想要第三个选项:<0表示它在平面的背面。如果所有4个平面都满足这个条件，你的点就在四面体内。

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

给定定义非简并四面体的4个点A,B,C,D，和一个测试点p，一种方法是将p的坐标变换为四面体坐标系，例如以A为原点，以向量B-A, C-A, D-A为单位向量。

在这个坐标系中，如果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))

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

多亏了Dorian的测试用例脚本，我可以在另一个解决方案上工作，并将其与目前为止的解决方案进行快速比较。

直觉

对于三角形ABC和点p，如果将p与角连接得到向量PA, PB,PC，并比较由PA,PC和PB,PC张成的两个三角形X和Y，则如果X和Y重叠，则点p位于三角形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 and detB <0和detC> 0和detD <0

或

detA & lt;0 and detB> 0 and detC <</p>

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

它工作吗?显然。为什么它会起作用?我不知道，至少我无法证明。也许其他数学更好的人可以帮我们。

(编辑:实际上重心坐标可以用这些行列式来定义最后，重心坐标的和需要等于1。这就像将P与点A,B,C,D组合成的四面体的体积与四面体ABCD本身的体积进行比较。观察决定符号的解释情况仍不清楚它是否一般有效，我不建议它)

i改变了测试用例，不再针对一个四面体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 ---