Search code examples
pythonperformancevectorizationfinite-element-analysistetrahedra

Given a point and a large number of tetrahedrons, how to efficiently determine in which tetrahedron the point is

Suppose we define a point to be a tuple of three floating-point numbers, and a tetrahedron a tuple of four points.

Suppose we have a tetrahedron and a point, we can determine whether the point belongs to the tetrahedron following the solutions described in How to check whether the point is in the tetrahedron or not? The key idea there is to determine whether the point is on the inner sides of the four flanks of the tetrahedron.

My problem. Given a point, and N tetrahedrons, where N is about 7 million, I need to determine in which tetrahedron the point is. We will care about performance of doing repeated tests, with a large number of points.

Additional info:

  1. One could just check these tetrahedrons one by one using the methods mentioned above. But that could be too slow, given my large number of tetrahedrons.

  2. There is a specific point in the problem setting. These tetrahedrons are obtained from a FEM (finite element method ) problem for solving a medical imaging problem (they form the brain of patients). Perhaps FEM itself is unrelated to the question, but we could leverage the fact that those tetrahedrons are next to each other and there are no “holes” in the space that simulated by those tetrahedrons.

  3. The tetrahedrons have no intersections except on their adjacent boundary. So, this question should have a unique solution unless at the boundary, in which case it is fine to have either of the intersected tetrahedrons the answer to my problem.

  4. There are no specific orders in which the tetrahedrons are given on inputs. There are no specification on whether the shapes of the tetrahedrons are regular or not.

Any idea on an efficient solution to the problem? Python is preferred in solving this problem.

Thanks!

Solution

  • You could first filter the tetrahedrons, keeping only those for which the bounding cuboid (which is parallel with the X, Y and Z axes) contains p. This is faster to test:

    So find tetrahedrons -- with points t0, t1, t2, and t3 -- which have the following property with respect to the point p:

    • ∃ i,j: tix ≤ px ≤ tjx
    • ∃ i,j: tiy ≤ py ≤ tjy
    • ∃ i,j: tiz ≤ pz ≤ tjz

    On average this will leave you with only a few tetrahedrons (often only one or two) which you then use to apply the point-in-tetrahedron test.