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

Question

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!

Answer 1

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 t₀, t₁, t₂, and t₃ -- which have the following property with respect to the point p:

• ∃ i,j: t_i^x ≤ p^x ≤ t_j^x
• ∃ i,j: t_i^y ≤ p^y ≤ t_j^y
• ∃ i,j: t_i^z ≤ p^z ≤ t_j^z

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.

Answer 2

If you plan to test a lot of points against same set of tetrahedra, I would definitely go with a preprocessing step and build spatial structure for tetrahedra.

In my comment I mentioned octtree, but knowing that the tetrahedra fill in the space (no holes) I think that there is no need for adaptive subdivision of space, and it is best to divide it into equal parts.

1. Divide the space into equal boxes (lets name them SpaceBoxes).
2. In each SpaceBox, keep a list of tetrahedra that collide with the box.

• To speed it up, I would test tetrahedron's bounding box, not tetrahedron itself.
• Note that this step can be done relatively cheaply - you know that SpaceBoxes have equal sizes, you know their position, so given tetrahedron's bounding box, it is easy to find SpaceBox candidates.

Now, having this spatial structure:

For point to be tested p

• find a corresponding SpaceBox O(1)
• you have all tetrahedra that collide with the SpaceBox, so these are candidates to test
• first test collision of p with the bounding box of each tetrahedron
• only then, with the tetrahedron itself

Note that the performance of test depends mostly on the amount of tetrahedra in each SpaceBox.

Assuming a space is a cube:

• subdividing each edge to 16 parts gives you 16^3 = 4096 SpaceBoxes
• having N = 7000000, there should be roughly 1709 candidate tetrahedra to test

Also, on the implementation side, both preprocessing and testing multiple points look like data-parallel problems, so multiprocessing may help.

Answer 3

Put the tetrahedra into a 3D R-tree such that their bounding boxes are the keys and the tetrahedra themselves are the values. Query the R-tree on the point. Then test each result value returned by the query with the point-in-tetrahedron check.

The nice thing about using an R-tree (or similar structure) is that if the tetrahedra do not change you can build the R-tree once and test many points.