How to check whether the point is in the tetrahedron or not?

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

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.


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