# Finding which triangles points are in

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!

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.

• Victor -- do you know of any open-source code that impliments the edge-hopping approach? IT look slike it coudl be a very good solution for my case. (particle tracking model driven by current fields in a traingualr mesh grid) -Thanks Jun 28, 2012 at 23:28
• I have code for this and I can send it to you; it's in C/C++. Haven't had time to clean it up and post it on Github yet. I've had to write this at least twice in my life, once with a halfedge data structure, again with a quadedge, but it can easily be used when those are not available and you need to build topological structure yourself. Look on my profile page for my website, where you can find contact info. We can discuss this further offline. Jun 29, 2012 at 4:45
• I'm close to finished implementing this in matlab using Hilbert curve ordering and randomized triangle walk. It's research code: not optimized, not documented, etc, but still pretty fast - I can give you the code if you're interested. Jun 29, 2012 at 18:13
• About : """ edge hopping should be O(logN) """ I don't see that. For example, in the pathological case of a big long triangle strip (like a narrow channel only on triangle wide), on the worst case, you'd need to hop from one triangle to the next all the way to the end. In the average case, half way. So if you double the number of triangles, it would be O(N) In the more normal case of a square-ish arrangement of triangles, I'd expect O( sqrt(N) ). Or am I missing something? -Chris Jul 2, 2012 at 18:58
• @Chris - Welcome to scicomp! As part of scicomp's housekeeping, I have migrated your answers and the ensuing conversation over as comments on Victor's answer. We're looking forward to your participation on the site. Jul 2, 2012 at 19:03

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.

• Hmm you're right. On the other hand, if the triangulation was Delaunay, I think nearest neighbor would work. It's too restrictive for what I'm trying to do, but in the Delaunay case consider the dual Voronoi diagram - the Voronoi cells are the set of points closest to a node, and the edges of the delaunay triangles all meet the edges of the Voronoi cells at right angles, so any point must be in a triangle connected to its closest node. I wonder if this is how matlab's pointLocation function works under the hood..? Jun 26, 2012 at 11:19

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