4

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!

zell
  • 9,830
  • 10
  • 62
  • 115
  • 2
    1. Is there any order in which the tetrahedra are given on input? (Perhaps the machine scans the brain in layers?) 2. Do I understand correctly that we care only about performance of one test, not repeated tests? (my first thought was octtree, but it may be too costly to build for one test) 3. How about shape of tetrahedra - most are regular or skewed? 4. I would start with trivial (minx, miny, minz) <= (px, py, pz) <= (maxx, maxy, maxz) to prune out candidates and only then do a full test. – Lesiak Sep 15 '20 at 20:45
  • Thanks for thinking about this. To answer your question: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. And we do care about performance of repeated tests. – zell Sep 16 '20 at 20:17

3 Answers3

3

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.

trincot
  • 317,000
  • 35
  • 244
  • 286
  • Thanks a lot! The condition of your filter makes sense to me. I am just wondering how did you get. "On average this will leave you with only a few tetrahedrons (often only one or two) ". Did you get that by experience or by theory? – zell Sep 17 '20 at 20:10
  • When thinking this over, the first idea was to sort all tetrahedrons by their least x-coordinate, and then it became clear that the point's x coordinate had to be at least the x that belonged to its encompassing tetrahedron. By extension you find this must be true also for y and z, and in both directions. As to the average of "one or two": that was just an educated guess while visualising such a space. Sure there can be extreme cases where many tetrahedrons share the same point, maybe even all, and the point to locate could be *that* point. But that is just extreme, and not the average case. – trincot Sep 17 '20 at 20:16
  • I see. Also, I suppose that the filtering constraint "(minx, miny, minz) <= (px, py, pz) <= (maxx, maxy, maxz) " that Lesiak proposed in his comments actually do the same as yours. Right? – zell Sep 17 '20 at 20:31
  • Yes, that is the cuboid test. The way I put it, you can sometimes save time: as soon as you find one x among the corners of the tetrahedron that is not greater than px, you can move to the next requirement, while for calculating minx, you have to always look at all four corners. Obvously you can preprocess the tetrahedrons and store their min/max values for x,y,z, and then it will be faster at the cost of some space. – trincot Sep 17 '20 at 21:01
  • I wouldn't be so quick to assume that a bounding box test is worthwhile. Each face of a tetrahedron can be converted to a normal _n_ plus an offset _d_ from the origin. Testing whether _p_ is in the right half-space for that face is then a single dot product (vectorizable!), followed by a comparison. And you need to test only 4 tetrahedron faces instead of 6 cube faces. Since branching is often the bottleneck due to branch (mis)prediction, this might actually be faster despite the additional arithmetic. – Thomas Aug 17 '22 at 12:25
3

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.

Lesiak
  • 22,088
  • 2
  • 41
  • 65
  • Thank you very much! How would you compare your filtering method with that trincot proposed below? I feel trincot's seems to be easier to implement, and trincot mentions that after filtering one would only get 1-2 tetrahedrons (you would get."1709 candidates "). How would you compare the efficiency of the two approaches? – zell Sep 17 '20 at 20:24
  • Besides, by "box", did you mean 3D-cube? Also, I am not sure what "edge" did you refer in "subdividing each edge to 16 parts"? But I think I get your idea -- preprocess first. – zell Sep 17 '20 at 20:25
  • BTW, I just realize that trincot's filtering constraint should be the same as what you proposed in your comments earlier, that is "(minx, miny, minz) <= (px, py, pz) <= (maxx, maxy, maxz) ". Do you think your SpaceBox idea would be more efficient? – zell Sep 17 '20 at 20:33
  • 1. By box i mean 3D cube (or rectangular cuboid, depending on dimensions of entire space). 2. subdividing each edge to 16 parts - I assumed that the space is 3D . You subdivide each dimension into 16 parts, so you get 16^3 `SpaceBoxes` 3. (minx, miny, minz) <= (px, py, pz) <= (maxx, maxy, maxz) is indeed the test for being in a bounding box (or bounding cuboid, if you prefer that name) – Lesiak Sep 17 '20 at 22:02
  • 4. trincots answer is exactly what I proposed in initial comment - make a pass through entire 7 million tetrahedra, but instead of testing collision with tetrahedron (which requires some numerical calculations) do a bounding box collision test first. Thus, hopefully, you can test with actual tetrahedra for a couple that passed bounding box test. The downside - for the next point you need to repeat 7 million bounding box checks. – Lesiak Sep 17 '20 at 22:10
  • 1
    5. My idea is to reduce the number of bounding box checks for each repeated test - so for each point you do only ~1700 (with 4096 boxes) bounding box tests and arrive to same filtered handful of candidates that passed bounding box test. 6. Definitely start with bounding box test first, before adding preprocessings step. As you noticed (and I mentioned in initial comment) this is trivial to implement. maybe you get desired speedup with this simple modification. – Lesiak Sep 17 '20 at 22:11
2

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.

jwezorek
  • 8,592
  • 1
  • 29
  • 46