Как проверить, находится ли точка в тетраэдре или нет?

Question

Я знаю все координаты тетраэдра и точку, которую я хотел бы определить. Так кто-нибудь знает, как это сделать? Я попытался определить принадлежность точки к каждому треугольнику тетраэдра, и если это верно для всех треугольников, то точка находится в тетраэдре. Но это абсолютно неправильно.

Answer 1

Вы определяете тетраэдр четырьмя вершинами. A B C и D Следовательно, вы также можете иметь 4 треугольника, определяющих поверхность тетраэдра.

Вы знаете, черт возьми, если точка P находится на другой стороне плоскости. Нормаль каждой плоскости направлена ​​в сторону от центра тетраэдра. Так что вам просто нужно протестировать против 4 самолетов.

Ваше уравнение плоскости выглядит так: a*x+b*y+c*z+d=0 Просто заполните значения точек (x y z). Если знак результата>0, точка находится на той же стороне, что и нормаль, результат == 0, точка лежит на плоскости, и в вашем случае вам нужен третий вариант. <0 означает, что он находится на задней стороне самолета. Если это выполнено для всех 4-х плоскостей, ваша точка находится внутри тетраэдра

Answer 2

Для каждой плоскости тетраэдра проверьте, находится ли точка на той же стороне, что и остальная вершина:

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 3

Начиная с решения Хьюга, вот более простой и (даже) более эффективный вариант:

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.

Для сравнения, вот полный код для сравнения методов, предложенных Нико, Хьюгом и мной:

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://stackru.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 ---

Answer 4

Учитывая 4 точки A,B,C,D, определяющие невырожденный тетраэдр, и точку P для проверки, одним из способов будет преобразование координат P в систему координат тетраэдра, например, взяв A в качестве начала координат, и векторы BA, CA, DA как единичные векторы.

В этой системе координат все координаты P находятся в диапазоне от 0 до 1, если он находится внутри P, но он также может находиться в любом месте преобразованного куба, определенного в начале координат и 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))

Answer 5

Я векторизовал решения Дориана и Хьюза, чтобы взять на вход весь массив точек. Я также переместил функцию 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)) 

Answer 6

Благодаря сценарию тестового примера Дориана я мог работать над еще одним решением и быстро сравнивать его с уже существующими.

интуиция

для треугольника ABC и точки P, если соединить P с углами, чтобы получить векторы PA, PB,PC, и сравнить два треугольника X и Y, которые охватывают PA,PC и PB,PC, то точка P находится внутри треугольник ABC, если X и Y перекрываются.

Или, другими словами, невозможно построить вектор PA путем линейного объединения PC и PB только с положительными коэффициентами, если P находится в треугольнике ABC.

С этого момента я попытался перенести его на случай тетраэдра и прочитал здесь, что можно проверить, являются ли векторы линейно независимыми, проверив, что определитель матрицы, построенной из векторов в виде столбцов, не равен нулю. Я пробовал разные подходы с использованием детерминантов и наткнулся на один:

Пусть 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. Объясненный случай с соблюдением детерминантных знаков до сих пор неясно, работает ли он в целом и я не рекомендую его)

Я изменил тестовый пример, чтобы не проверять 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://stackru.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 ---