Finding which triangles points are in

Question

Suppose I have a 2D mesh consisting of nonoverlapping triangles $\{T_k\}_{k=1}^N$, and a set of points $\{p_i\}_{i=1}^M \subset \cup_{k=1}^N T_K$. What is the best way to determine which triangle each of the points lies in?

For example, in the following image we have $p_1 \in T_2$, $p_2 \in T_4$, $p_3 \in T_2$, so I would like a function $f$ that returns the list $f(p_1,p_2,p_3)=[2,4,2]$.

Matlab has the function pointlocation which does what I want for Delaunay meshes, but it fails for general meshes.

My first (dumb) thought is, for all nodes $p_i$, loop through all the triangles to find out which triangle $p_i$ is in. However, this is is extremely inefficient - you might have to loop through every triangle for every point, so it could take $O(N \cdot M)$ work.

My next thought is, for all points $p_i$, find the nearest mesh node via nearest-neighbor search, then look through triangles attached to that nearest node. In this case, the work would be $O(a\cdot M\cdot log(N))$, where $a$ is the maximum number of triangles attached to any node in the mesh. There are a couple solvable but annoying issues with this approach,

• It requires implementing an efficient nearest-neighbor search (or finding a library that has it), which could be a nontrivial task.
• It requires storing a list of which triangles are attached to each node, which my code is currently not set up for - right now there is just a list of node coordinates and a list of elements.

Altogether it seems inelegant, and I think there should be a better way. This must be a problem that arises a lot, so I was wondering if anyone could recommend the best way to approach finding what triangles the nodes are in, either theoretically or in terms of available libraries.

Thanks!

Answer 1

The usual randomized edge hopping method should work. Basically, start with any triangle of the mesh, then determine which of the edges the target point lies on the opposite side of. That is, determine which of the edges, when extended out to a line, separate the point from the interior of the triangle. When there are two possibilities, choose one at random, and consider the triangle that is adjacent to that shared edge, and repeat. The randomization should make this method converge with probability 1 for Delaunay triangulations, and I can't think of a reason it would not work for arbitrary triangulations.

Edit: I should add that edge hopping should be $O(\log N)$ with a reasonable constant for a single point, so it would be $O(M \log N)$ for $M$ points. However, if you sort your points by locality (like using a Hilbert curve ordering first), you can initialize each new query with the triangle of the previous query, to further reduce the runtime (I'm not a CS theorist so I can't tell you what the big-O would be there).

Edit2: Found this PDF describing such a "walking" scheme that is guaranteed to terminate, and reviews the more naive approaches.

Another alternative to using quadtrees is using a Triangulation Hierarchy. See Olivier Devillers. Improved incremental randomized Delaunay triangulation. In Proc. 14th Annu. ACM Sympos. Comput. Geom., pages 106-115, 1998. It works best for Delaunay triangulations, but can also work for non-Delaunay.

Basically whatever you do to speed up point location will require the construction of an auxiliary data structure. In the case of quadtrees or some other spatial subdivision, you need to build the subdivision tree. In the case of edge-hopping, you need to build the triangle adjacent topological structure. The triangulation hierarchy also requires building a tree of coarser triangulations.

Answer 2

I'm not convinced your solution is actually correct. Consider the situation where you have these nodes:

• A: (-3, 1)
• B: (0, 2)
• C: (3, 1)
• D: (0, -5)

There are triangles ABC and ACD. Now B is the closest point to the origin, but the origin is in triangle ACD, which doesn't contain B.

For the rest, this answer provides a solution that may in degenerate cases be as slow as $O(N \cdot M)$, but will typically be faster - especially for Delaunay triangulations and triangulations close to that in some sense.

I would consider the option of building a quadtree that contains the triangles themselves. I.e., you have a quaternary tree that stores in each node (which corresponds to a bounding box):

• The coordinates at which the box is being split, or alternatively, the bounding boxes of the four subtrees;
• Pointers to the subtrees;
• The set of triangles that fall completely within the bounding box of this rectangle, but not completely within any of the four subtrees. In other words, the triangles that intersect with any of the two subdividing line segments of the quadtree.

When given a point P, traverse all nodes on the path from the root of the quadtree to the smallest box containing P. You'll need to examine all triangles that you encounter in those nodes. For a 'well-behaved' triangulation, there should only be something like $\sqrt{n}$ triangles that need to be examined at the level of a node that contains $n$ triangles in its subtree, and the depth should be limited by $\log n$. For a 'badly behaved' triangulation, you may get the $O(N\cdot M)$ work.

Answer 3

CGAL handles triangulations and point location (and a great deal more). In particular, the 2D triangulation module can do what you want.