CGAL - Delaunay Triangulation Hierarchy with Info

Question

I have a PCL Point Cloud (in 3D) which I would like to turn into a ground surface TIN (2.5D) and then sample points (in 2D) to find their elevation when projected onto the TIN. To do this I've been using the CGAL Delaunay Triangulation classes, which has mostly been working pretty well!

I was able to implement this using a Delaunay_triangulation_2 built on a Triangulation_vertex_base_with_info_2 and create a nice-looking TIN. I also wrote a function which extracts the face and vertices for an arbitrary point in 2D space using the CGAL locate() function so that I can interpolate the point's height if it were projected onto the TIN. I need the info field to hold an index that lets me associate vertices within the triangulation back to points in the PCL point cloud structure.

However, when using the base Triangulation class the locate() function is slow (stochastic walk starting from an arbitrary vertex in the triangulation), and because I have to call this for every query point in the cloud to be interpolated (to estimate projected heights) this is currently the slowest part of my whole pipeline. So I looked into using the Triangulation Hierarchy class to make this more efficient.

I can't figure out how to make the Triangulation_hierarchy class work with a vertex base with info, and I think I'm just doing something silly wrong. Here is a miminal example showing my slow solution with a simple triangulation structure (no hierarchy), which does work:

#include <chrono>

#include <pcl/point_types.h>
#include <pcl/filters/passthrough.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Delaunay_triangulation_2.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel             K;
typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, K>    Vb;
typedef CGAL::Triangulation_data_structure_2<Vb>                        Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds>                          Delaunay;
typedef Delaunay::Point_2                                               CGALPoint;
typedef Delaunay::Face_handle                                           Face_handle;


// This is templated on a PCL PointCloud Ptr - input cloud is basically a 3D vector of XYZ points
template <typename CloudType> 
void delaunayTriangulation(CloudType input_cloud, Delaunay& triangulation)
{ 
    std::cout << "Performing Delaunay triangulation on cloud of size " << input_cloud->points.size() << std::endl;
    // Convert ground minima cloud to CGAL vector of points
    std::vector< std::pair<CGALPoint, unsigned> > minima_vec;
    for(std::size_t i=0; i<input_cloud->points.size(); i++)
    {
        minima_vec.push_back(std::make_pair(CGALPoint(input_cloud->points[i].x,input_cloud->points[i].y), i));
    }
    // Generate Delaunay Triangulation for ground minima 
    triangulation = Delaunay(minima_vec.begin(), minima_vec.end());
    std::cout << "  Number of vertices in Delaunay: " << triangulation.number_of_vertices() << std::endl;
    std::cout << "  Number of faces in Delaunay: " << triangulation.number_of_faces() << std::endl; 
}

int main()
{
    // Generate a starting point cloud with random points
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
    for(int i=0; i<500; i++)
        for(int j=0; j<500; j++)
        {
            // Generate points which are gridded + a bit of noise in XY, and random Z 
            pcl::PointXYZ point;
            point.x = i + (std::rand()%100)/100.0;
            point.y = j + (std::rand()%100)/100.0; 
            point.z = std::rand();
            cloud->points.push_back(point);
        }
    // Get the ground triangulation
    Delaunay triangulation;
    delaunayTriangulation(cloud, triangulation);
    // Locate the containing face for a bunch of random points
    std::cout << "Starting to search for faces..." << std::endl;
    auto start_time = std::chrono::high_resolution_clock::now();
    for(int i=0; i<3000000; i++)
    {
        // Random point with X and Y between 0 and 500
        CGALPoint test_point((std::rand() % 5000)/10.0, (std::rand() % 5000)/10.0);
        Face_handle face = triangulation.locate(test_point);
        // here we would do some math using the vertices located above
    }
    auto stop_time = std::chrono::high_resolution_clock::now();
    float duration = std::chrono::duration_cast<std::chrono::microseconds>( stop_time - start_time ).count();
    std::cout << "Finished searching for faces - took " << duration/1000000 << std::endl;
}

If I switch to using a Triangulation_hierarchy_2 object instead, built on top of my Delaunay_Triangulation_2 type, it won't let me insert point pairs which contain the info field - it will only compile if I build the object using a vector of points by themselves, instead:

#include <chrono>

#include <pcl/point_types.h>
#include <pcl/filters/passthrough.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_hierarchy_2.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel             K;
typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, K>    Vbb;
typedef CGAL::Triangulation_hierarchy_vertex_base_2<Vbb>                Vb;
typedef CGAL::Triangulation_data_structure_2<Vb>                        Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds>                          Delaunay;
typedef Delaunay::Point_2                                               CGALPoint;
typedef Delaunay::Face_handle                                           Face_handle;
typedef CGAL::Triangulation_hierarchy_2<Delaunay>                       Delaunay_hierarchy;

// This is templated on a PCL PointCloud Ptr - input cloud is basically a 3D vector of XYZ points
template <typename CloudType> 
void delaunayTriangulation(CloudType input_cloud, Delaunay_hierarchy& triangulation)
{ 
    std::cout << "Performing Delaunay triangulation on cloud of size " << input_cloud->points.size() << std::endl;
    // Convert ground minima cloud to CGAL vector of points
    std::vector<CGALPoint> minima_vec_simple;
    for(std::size_t i=0; i<input_cloud->points.size(); i++)
    {
        minima_vec_simple.push_back(CGALPoint(input_cloud->points[i].x,input_cloud->points[i].y));
    }
    // Generate Delaunay Triangulation for ground minima 
    triangulation = Delaunay_hierarchy(minima_vec_simple.begin(), minima_vec_simple.end());
    std::cout << "  Number of vertices in Delaunay: " << triangulation.number_of_vertices() << std::endl;
    std::cout << "  Number of faces in Delaunay: " << triangulation.number_of_faces() << std::endl; 
}

int main()
{
    // Generate a starting point cloud with random points
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
    for(int i=0; i<500; i++)
        for(int j=0; j<500; j++)
        {
            // Generate points which are gridded + a bit of noise in XY, and random Z 
            pcl::PointXYZ point;
            point.x = i + (std::rand()%100)/100.0;
            point.y = j + (std::rand()%100)/100.0; 
            point.z = std::rand();
            cloud->points.push_back(point);
        }
    // Get the ground triangulation
    Delaunay_hierarchy triangulation;
    delaunayTriangulation(cloud, triangulation);
    // Locate the containing face for a bunch of random points
    std::cout << "Starting to search for faces..." << std::endl;
    auto start_time = std::chrono::high_resolution_clock::now();
    for(int i=0; i<3000000; i++)
    {
        // Random point with X and Y between 0 and 500
        CGALPoint test_point((std::rand() % 5000)/10.0, (std::rand() % 5000)/10.0);
        Face_handle face = triangulation.locate(test_point);
        // here we would do some math using the vertices located above
    }
    auto stop_time = std::chrono::high_resolution_clock::now();
    float duration = std::chrono::duration_cast<std::chrono::microseconds>( stop_time - start_time ).count();
    std::cout << "Finished searching for faces - took " << duration/1000000 << std::endl;
}

All the layers of typedefs and templates in CGAL are a bit opaque to me - is there some way I can set these up instead which would allow me to build something like a Triangulation_hierarchy_2 based on a Triangulation_hierarchy_vertex_base_with_info_2? That second type doesn't look like it exists.

Thanks!

Answer 1

Ok - I ended up trying a few routes, and did some basic benchmarking of execution time for each. These are based on the settings in the code in my question:

• TIN with 25,000 points spaced mostly on a grid from 0 to 500 in XY, with a bit of noise in XY, and random z values
• Test cloud with 3,000,000 points with random values between 0 and 500 in X/Y

What I tried:

1. CGAL locate() using simple Triangulation structure, arbitrary starting face --> 79.3 s
2. CGAL locate() using a Triangulation Hierarchy structure, arbitrary starting face --> 4.74 s
3. CGAL locate() using a starting face given by the nearest neighbor vertex, found using a K-D tree on the input cloud of vertices --> 3.41 s

To do this, after creating the triangulation I iterated over faces and made a mapping between input cloud vertex indices and triangulation face handles:

std::vector<Face_handle> face_mapping(cloud->points.size());
std::vector<bool> faces_filled(cloud->points.size(), false);
auto start_time_list = std::chrono::high_resolution_clock::now();
// Iterate over all faces in triangulation
for (Face_handle face : triangulation.finite_face_handles())
    // Iterate over 3 vertices for each face
    for(int i=0; i<3; i++)
    {
        int index = uint32_t(face->vertex(i)->info());
        if(!faces_filled[index])
            face_mapping[index] = face;
    }

Then when we run the point locating search we can do something like this:

std::vector<int> nearest_indices;
std::vector<float> nearest_dists;
pcl::KdTreeFLANN<pcl::Point2DGround> tree;
pcl::PointCloud<pcl::Point2DGround>::Ptr cloud_2d(new pcl::PointCloud<pcl::Point2DGround>);
copyPointCloud3D(cloud, cloud_2d);
tree.setInputCloud(cloud_2d);
for(int i=0; i<3000000; i++)
{
    // Random point with X and Y between 0 and 500
    CGALPoint test_point((std::rand() % 5000)/10.0, (std::rand() % 5000)/10.0);
    
    // Get closest vertex in triangulation
    pcl::Point2DGround pcl_point;
    pcl_point.x = test_point.x();
    pcl_point.y = test_point.y();
    tree.nearestKSearch(pcl_point, 1, nearest_indices, nearest_dists);

    Face_handle face = triangulation.locate(test_point, face_mapping[nearest_indices[0]]);

    // here we would do some math using the vertices located above
}

Where pcl::Point2DGround has to be some custom point type where the representation is only 2D (to make the search tree work properly).

I didn't end up trying any kind of sort on the incoming points to be tested, because the above works and because my test points to be locate()-ed will always be much more numerous than the TIN vertices themselves, so I think it would probably be more expensive to impose order on them instead.

I think this is the simplest/fastest solution for my case, so I'll go with this! Thanks @marcglisse and @andreasfabri for your comments.