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

Question

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

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

Answer 2

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

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

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

Answer 3

从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)) # /sf/ask/109030001/

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