Intersection of tetrahedra

Question

I hope to be as clear as possible. I'm trying to implement a function that, given two tetrahedra, checks if they intersect with each other. I am working with python and the only library I am using is NumPy. To describe a tetrahedron I use its 4 vertices which are each described by coordinates [x, y, z].

vertex = [x, y, z]

tetrahedra = [vertex 1,vertex 2,vertex 3,vertex 4]

This is the reasoning I want to use:

• A tetrahedron is nothing more than a region defined by inequalities.
• These inequalities are described by the planes containing one face of the tetrahedron.
• So given the inequalities of the two tetrahedra, and putting them in a system, if this system admits a solution then there is an intersection.

This is my function:

def IsInterpenetrated(self, tetrahedra):
    A= []
    B= []
    sol= 0
    for tr in [self, tetrahedra]:
        print("Plane of tetrahedra")
        vertexList = tr.vertices
        i=0
        while i<4:
            if handedness(vertexList)>0:
                n= numpy.cross(vertexList[1].coords - vertexList[0].coords, vertexList[2].coords - vertexList[0].coords)
            else:
                n= numpy.cross(vertexList[2].coords - vertexList[0].coords, vertexList[1].coords - vertexList[0].coords)
            
            p0= vertexList[0].coords
            d= -(n[0]*p0[0] + n[1]*p0[1] + n[2]*p0[2])
            
            print("normal: ", n , end="      ")
            print("termine noto: ",(d))

            if len(A) > 3:
                j=0
                while j<=3:
                    if numpy.all(-n == A[j]) and -d == B[j]:
                        sol = 1
                    j= j+1

            A.append(n)
            B.append(d)

            p0= vertexList[0]
            vertexList[0] = vertexList[1]
            vertexList[1] = vertexList[2]
            vertexList[2] = vertexList[3]
            vertexList[3] = p0

            i=i+1

    A= numpy.array(A)
    B= numpy.array(B)
    print("\n")

    print("Disequazioni:\n")
    i=0
    for n in A:
        print("({0})x + ({1})y + ({2})z + ({3}) > 0".format(n[0],n[1],n[2],B[i]))
        i=i+1
    print("\n")

    x = cvxpy.Variable(3)
    prob = cvxpy.Problem(cvxpy.Minimize(0),[A @ x + B >= 0])
    prob.solve()
    if prob.value == 0 and sol != 1:
        return 1
    return 0

In this case I have solved the system of inequalities using cvxpy and I have verified the particular case in which the two tetrahedra have a common face. I was wondering if you think the following reasoning is correct to avoid working with systems of inequalities. Each plane that identifies the face of the tetrahedron belongs to the family of bundles of parallel planes which are described in the following way; ax + by + cz + k = 0 where k is the term that indicates the exact position of the plane on space. Then I can describe the tetrahedron in the following way:

System:
a'x + b'y + c'z = k '
a "x + b" y + c "z = k"
a '"x + b'" y + c '"z = k'"
a "" x + b "" y + c "" z = k ""

with k > d where d is the known term of the plane that identifies the face. Thanks to the Rouché-Capelli theorem I know that this system admits solution if Rg (A) = Rg (A | B) where Rg stands for rank. To ensure that this equality is respected then Det (A | B) = 0 where Det stands for determinant. Since B in my case consists of variables:

(k ', k ", k"', ......, kᵐ)

then to make Det (A | B) = 0 I have to solve the equation that is created by this calculation. Having carried out this reasoning for both tetrahedra, I find myself with two equations with 3 unknowns. One for each tetrahedron. By putting these two equations into a system I have to see for which values of k it admits solutions. If there are values of k for which the system is respected then I have intersection, otherwise no. I don't know how feasible it is but I preferred to share my idea, in order to discuss it together.

Thanks in advance.

Answer 1

Why not formulate a convex optimization problem, or precisely a feasibility problem using the plane inequalities that you have? Let's say, the two tetrahedra can be represented as A1.X + d1 <= 0 and A2.X + d2 <= 0 where the 4 rows of A1 and A2 store the a, b, c of four planes corresponding to the two tetrahedra in ax + by + cz + d <= 0, and column vectors d1 and d2 store the constants ie d. And also note that A1.X is matrix multiplication.

Represent (x, y, z) as the vector X.

Now basically you want to solve a feasibility problem for X like this:

minimize 0
subject to A1.X + d1 <= 0
           A2.X + d2 <= 0

Note that if the solver returns inf, that means there is no X which satisfies the above constraints. If the solver returns 0 (which is the value of the constant objective function), that means there is atleast one X which satisfies the constraints.

You can use cvxpy library for this. Here is a nice tutorial. Also cvxpy library goes well with numpy.

And I don't think solving equations would work in this case as a tetrahedron is basically composed of four linear inequalities. So you have to solve inequalities in order to find a solution in their intersection region.