Fast way to find the triangle inside a mesh that encloses a point

Question

I'm running into a performance problem for a task I need to accomplish. One of the bottlenecks at the moment is in getting an interpolated field value from an unstructured grid.

The slow part is, given a 2D point and an unstructured 2D grid, find the mesh points which immediately surround the point. It would be nice to just find the triangle it falls into.

Right now I'm using CGAL, but it's way too slow. With the current implementation, the whole task will take days to complete, running in parallel on a high end CPU.

I believe that the slow part is CGAL::natural_neighbor_coordinates_2.

#ifndef FIELD_INTERPOLATOR_H
#define FIELD_INTERPOLATOR_H

#include "Vec.h"

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>

#include <map>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Delaunay_triangulation_2< Kernel > Delaunay_triangulation;

typedef Kernel::FT FieldType;
typedef Kernel::Point_2 MeshType;

struct FieldInterpolator23 {

    Delaunay_triangulation m_triangulation;

    std::map< MeshType, FieldType, Kernel::Less_xy_2 > m_vX;
    std::map< MeshType, FieldType, Kernel::Less_xy_2 > m_vY;
    std::map< MeshType, FieldType, Kernel::Less_xy_2 > m_vZ;

    typedef CGAL::Data_access< std::map< MeshType, FieldType, Kernel::Less_xy_2 > > ValueAccess;

    FieldInterpolator23() {}

    FieldInterpolator23( const std::vector< TN::Vec2 > & mesh, const std::vector< TN::Vec3 > & field )
    {
        const int N = mesh.size();
        for ( int i = 0; i < N; ++i ) {

            MeshType p( mesh[i].x(), mesh[i].y() );

            m_triangulation.insert( p );
            m_vX.insert( std::make_pair( p, field[i].x() ) );
            m_vY.insert( std::make_pair( p, field[i].y() ) ); 
            m_vZ.insert( std::make_pair( p, field[i].z() ) );                        
        }       
    }

    void set( const std::vector< TN::Vec2 > & mesh, const std::vector< TN::Vec3 > & field ) {

        m_triangulation.clear();
        m_vX.clear();
        m_vY.clear();
        m_vZ.clear();

        const int N = mesh.size();
        for ( int i = 0; i < N; ++i ) {

            MeshType p( mesh[i].x(), mesh[i].y() );

            m_triangulation.insert( p );
            m_vX.insert( std::make_pair( p, field[i].x() ) );
            m_vY.insert( std::make_pair( p, field[i].y() ) );
            m_vZ.insert( std::make_pair( p, field[i].z() ) );
        }
    }

    TN::Vec3 operator() ( TN::Vec2 p ) {

        MeshType pos( p.x(), p.y() );

        std::vector< std::pair< MeshType, FieldType > > coords;

        FieldType norm =
            CGAL::natural_neighbor_coordinates_2( m_triangulation, pos, std::back_inserter( coords ) ).second;

        FieldType resX =
            CGAL::linear_interpolation(
                coords.begin(), 
                coords.end(),
                norm,
                ValueAccess( m_vX )
        );

        FieldType resY =
            CGAL::linear_interpolation(
                coords.begin(), 
                coords.end(),
                norm,
                ValueAccess( m_vY )
        );

        FieldType resZ =
            CGAL::linear_interpolation(
                coords.begin(), 
                coords.end(),
                norm,
                ValueAccess( m_vZ )
        );

        return TN::Vec3( resX, resY, resZ );
    }
};

#endif

Can anyone point me in the direction of an acceptable higher performance solution, either a different library or an algorithm?

Fatback

What's the number of vertices in your triangulation? How many query points do you need to locate in the triangulation? – Journalist

There are about 1 million mesh points, and I need to do about 2 billion queries. – Fatback

If you know the queries in advance, it helps to sort them (say along a space-filling curve, CGAL has functions for that) and then do the point location queries in this order, passing each query the previous location as a starting point. This way most queries will only check that the starting point is already the right triangle, and others will have at most 1 or 2 steps to walk to reach it. – Rosalba

If your query points are on a regular grid, a fast approach is by triangle filling: create an image where a pixel corresponds to a query point, and fill every triangle using the triangle ID as the fill color. This will make the queries constant-time. (This may be impractical in your case for such a big number of queries - unless you can work by sub-windows.) – Nan

Answer 1

CGAL contains an implementation of a Triangulation Hierarchy which

implements a triangulation augmented with a data structure to efficiently answer point location queries. [...] the data structure remains small and achieves fast point location queries on real data.

Its performance is optimal for Delaunay triangulations.

---

         
          ^{Fig. 36.8}

---

Duke

I switched over to the hierarchical triangulation, but it only resulted in a slight performance increase. As a test I tried naively iterating through the entire set of mesh nodes to find the closest one and return it's value. With this method, the whole task completes about 10 times faster than using a CGAL triangulation. – Fatback

@Fatback That sounds strange, but I don't quite understand what you are comparing. Do you have a benchmark you could share that does only the point location and shows naive iteration to be faster than triangulation's methods? – Rosalba

@MVTC: It is strange that a brute-force linear search outperforms a logarithmic search. I assume you are setting aside preprocessing and only looking at query times? – Anniceannie

Preprocessing is set aside. Also, using the locate function to find the face enclosing the point is also much much faster than the natural neighbours function in my application, which is strange because the natural_neighbors_coord_2 function seams to just uses locate internally. The time differences are so vast that it would seam something must be wrong with the triangulation? I will need more time to investigate what's going on. I'll mark this as the accepted answer, and then give an update if I figure out why this is happening. – Fatback

Answer 2

In my own experience, Axis Aligned Bouding Box trees are reasonably efficient for this problem (and I observed more performance than "walking in the triangulation", i.e. using locate()).

I am using my own implementation (for 3D mesh, but could be easily adapted to 2D): http://alice.loria.fr/software/geogram/doc/html/classGEO_1_1MeshFacetsAABB.html

AABB tree is also implemented in CGAL: http://doc.cgal.org/latest/AABB_tree/

Note: if your query points are structured (e.g. organized on a regular grid), then there are several ways of gaining a lot of performance, for instance by:

• using the locate() function of Delaunay triangulation in CGAL that takes a hint as an argument, and for the hint, using one of the triangles incident to the previous point. This ensures that the point location algorithm will not have to walk too far away. The gain is in general very significant with this approach. If you do not want to change your class API, you can have a mutable member that stores the hint (works also if the points are not structured, but requires to sort them spatially, see also Marc Glisse's comment).

• "drawing" the triangles over the points (each triangle "sends" its values to the points it covers). This can be done by "computer graphics" algorithms used to draw triangles on the screen. The gain with this approach is even more important (but requires more work, it will completely change your algorithm).

Flour

Answer 3

If the mesh grid is invariant, a hash table can be used to access any entry in O(1) time. Assuming the data is accessed by horizontal x and vertical y values, place a bounding box around the mesh and slice it vertically and horizontally into squares roughly equal to the area of the average mesh element. Set Nx to the number of vertical slices and Ny to the number of horizontal slices. If the mesh data extends from X_min to X_max and Y_min to Y_max, set Sx = Nx/(X_max—Xmin) and Sy = Ny/(Y_max—Y_min). Then incrementally number the vertical columns from left to right, starting at round(Sx*X_min) to round(Sx*X_max). Similarly incrementally number the horizontal rows from bottom to top, starting at round(Sy*Y_min) to round(Sy*Y_max).

If the mesh is roughly square there would be about 1000 columns and 1000 rows and nearly every square would be associated with a triangle. If the mesh has irregular shape or holes, many of the squares may not fall on the mesh. That makes no difference to the hash table. To set up the hash table, no more than one triangle can be associated with each square. If two or more triangles want to be associated with the same square, the squares are too large. To get more squares, set a smaller area to get larger Nx and Ny. Now that each triangle has been associated with exactly one square in some column #X and row #Y, a million keywords can be generated. For the cell in column #X and row #Y set up the keyword string "±#X±#Y" where ±#X and ±#Y mean integers with a leading "+" sign if the number is greater than —1 and a leading "—" sign if the number is less than zero. To give all keywords the same length, pad the integer values with leading zeros. Typical unique keywords might look like "-0378+0087" or "+0029-1007". In these examples every keyword would be exactly 10 characters. Set the keywords in the hash table in the order the triangles are stored.

To use the hash table, for a given x and y point, set integer ix = round(Sx*x) and integer iy = round(Sy*y) and form the keyword "±ix±iy" according to the sign of ix and iy. If x and y are in the mesh, the keyword "±ix±iy" should be in the hash table, and it should return the index of the triangle enclosing x and y or a triangle very close to the point. In some cases the generated "±ix±iy" keyword may not be in the hash table. In such cases the hash table function will return some random triangle index. When this happens the best fix is to perturb ix and/or iy by one unit and generate another keyword. Of course if the point is outside the mesh, the hash table may be sending back a lot of random triangle indexes. So make sure the point is valid to start with.

This discussion has assumed that the triangles were reasonably shaped and of approximate uniform area. If that is not the case and several triangles must fit into the same square or several squares must fit into the same triangle, the fix is to produce a redirection table on the side. When several triangles fit into the same square, assign the most central triangle to the square with the understanding that a quick search for the desired triangle will then be required. When several squares must fit into the same triangle, point one of the squares to the triangle index and all the others to an indirection table whose indices start after the triangle indices. The indirection table will then be used to point to the triangle index containing those squares.

Fabulist