How to check whether the point is in the tetrahedr

问题:

I know all coordinates of tetrahedron and the point I would like to determine. So does anyone know how to do it? I've tried to determine the point's belonging to each triangle of tetrahedron, and if it's true to all triangles then the point is in the tetrahedron. But it's absolutely wrong.

回答1:

You define a tetrahedron by four vertices, A B C and D. Therefore you also can have the 4 triangles defining the surface of the tetrahedron.

You now just check if a point P is on the other side of the plane. The normal of each plane is pointing away from the center of the tetrahedron. So you just have to test against 4 planes.

Your plane equation looks like this: a*x+b*y+c*z+d=0 Just fill in the point values (x y z). If the sign of the result is >0 the point is of the same side as the normal, result == 0, point lies in the plane, and in your case you want the third option: <0 means it is on the backside of the plane. If this is fulfilled for all 4 planes, your point lies inside the tetrahedron.

回答2:

For each plane of the tetrahedron, check if the point is on the same side as the remaining vertex:

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

And you need to check this for each plane:

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

回答3:

Given 4 points A,B,C,D defining a non-degenerate tetrahedron, and a point P to test, one way would be to transform the coordinates of P into the tetrahedron coordinate system, for example taking A as the origin, and the vectors B-A, C-A, D-A as the unit vectors.

In this coordinate system, the coordinates of P are all between 0 and 1 if it is inside P, but it could also be in anywhere in the transformed cube defined by the origin and the 3 unit vectors. One way to assert that P is inside (A,B,C,D) is by taking in turn as origin the points (A, B, C, and D) and the other three points to define a new coordinate system. This test repeated 4 times is effective but can be improved.

It is most efficient to transform the coordinates only once and reuse the SameSide function as proposed earlier, for example taking A as the origin, transforming into the (A,B,C,D) coordinates system, P and A must lie on the same side of the (B,C,D) plane.

Following is a numpy/python implementation of that test. Tests indicate this method is 2-3 times faster than the Planes method.

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

回答4:

Starting from Hugues' solution, here is a simpler and (even) more efficient one:

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)

In the coordinate system associated to the tetrahedron, the opposite face from the origin (denoted v1 here) is characterized by x+y+z=1. Thus, the half space on the same side as v1 from this face satisfies x+y+z<1.

As a comparison, here is the full code for comparing the methods proposed by Nico, Hugues and me:

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

And here are the results, in terms of running time:

--- 15.621951341629028 seconds ---
--- 8.97989797592163 seconds ---
--- 4.597853660583496 seconds ---