我知道四面体的所有坐标和我想要确定的点。有人知道如何做吗？我试图确定该点是否属于四面体的每个三角形，如果每个三角形都符合条件，则该点在四面体内。但这是完全错误的。
6个回答
26
对于四面体的每个平面，检查该点是否在剩余顶点的同侧：
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);
}
- Nico Schertler
6
你可以通过四个顶点 A、B、C 和 D 定义一个四面体，因此你也可以有定义四面体表面的 4 个三角形。
现在你只需检查一个点 P 是否在平面的另一侧。每个平面的法线指向四面体中心的方向，因此你只需对 4 个平面进行测试。
你的平面方程如下：
现在你只需检查一个点 P 是否在平面的另一侧。每个平面的法线指向四面体中心的方向，因此你只需对 4 个平面进行测试。
你的平面方程如下：
a*x+b*y+c*z+d=0。只需填入点的值（x y z）。如果结果的符号为 >0，则该点与法线在同一侧；结果 == 0，表示该点位于平面上，在你的情况下，你需要第三个选项：<0 表示它在平面的背面。如果这对于所有 4 个平面都成立，则你的点在四面体内。- glethien
5
如果我理解正确的话，首先我必须为每个平面找到A、B、C、D并计算ax+by+c*z+d。如果对于每个平面结果<0，则该点在四面体内。是这样吗？ - Denis Lolik
1不，命名是巧合。A、B、C、D是四面体的四个点，与方程中的a、b、c、d不同。 - glethien
是的，我指的是方程中的a、b、c、d，而不是A、B、C、D。抱歉。 - Denis Lolik
这个答案假设你可以知道法线的方向，而不是一般情况。 - Ander Biguri
如果我有三角形的三个点，我可以通过叉积计算法线。 - glethien
5
从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://dev59.com/oHI_5IYBdhLWcg3wBuX3
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 ---
[编辑]
基于汤姆的想法，如果想要找到网格中包含给定点的哪个元素，这里有一个相当高度矢量化的解决方案:
输入数据:
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://dev59.com/blgR5IYBdhLWcg3wBpZT#41851137
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倍。对于更大的网格，这种差距可能会增加。
- Dorian
4
给定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实现。测试表明，此方法比平面法快2-3倍。
在这个坐标系中，如果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))
- Hugues Talbot
3
我将Dorian和Hughes的解决方案向量化以便接受整个点数组作为输入。我还将tetraCoord函数移到了pointsInside函数之外，并重新命名两者，因为每个点都调用它是没有意义的。
在我的电脑上，@Dorian的解决方案和示例运行时间为2.5秒。在相同的数据上，我的解决方案运行速度快近1000倍，只需0.003秒。如果需要更快的速度，导入GPU cupy包并将其重命名为“np”可以将其加速到100微秒范围内。
在我的电脑上，@Dorian的解决方案和示例运行时间为2.5秒。在相同的数据上，我的解决方案运行速度快近1000倍，只需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))
- TomNorway
1
多亏了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类似）。
如果：
detA > 0且detB < 0且detC > 0且detD < 0，则点P位于由ABCD张成的四面体内
或者
detA < 0且detB > 0且detC < 0且detD > 0，则行列式从负数开始交替变换符号，或从正数开始。
这能够起作用吗？显然可以。为什么会起作用？我不知道，或者至少我不能证明它。也许其他数学技能更好的人可以在这里帮助我们。
（编辑：实际上，重心坐标可以使用这些行列式来定义，在最后，重心坐标需要总和为1。就像比较由P与点A、B、C、D的组合延伸出的四面体的体积与四面体ABCD本身的体积。解释行列式符号的情况仍然不清楚是否适用于一般情况，我不建议这样做）我将测试用例更改为不检查n个点Pi与一个四面体T相比，而是检查n个点Pi与n个四面体Ti相比。所有答案仍然给出正确结果。我认为这种方法更快的原因是它不需要矩阵求逆。
我保留了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://dev59.com/oHI_5IYBdhLWcg3wBuX3
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 ---
- Philipp B.
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接