CGAL - 带信息的 Delaunay 三角剖分层次结构
CGAL - Delaunay Triangulation Hierarchy with Info
我有一个 PCL 点云（3D），我想将其转换为地面 TIN（2.5D），然后采样点（2D）以在投影到锡。为此，我一直在使用 CGAL Delaunay 三角剖分 classes，它在大多数情况下都运行良好！
我能够使用基于 Triangulation_vertex_base_with_info_2 构建的 Delaunay_triangulation_2 来实现这一点，并创建一个漂亮的 TIN。我还编写了一个函数，它使用 CGAL locate() 函数提取 2D space 中任意点的面和顶点，这样我就可以在投影到 TIN 上时插入点的高度。我需要信息字段来保存一个索引，让我可以将三角剖分中的顶点关联回 PCL 点云结构中的点。
然而，当使用基本三角剖分时 class locate() 函数很慢（从三角剖分中的任意顶点开始的随机游走），并且因为我必须为要插值的云（以估计预计高度）这是目前我整个管道中最慢的部分。所以我研究了使用三角剖分层次结构 class 来提高效率。
我不知道如何使 Triangulation_hierarchy class 与带有信息的顶点基一起工作，我认为我只是在做一些愚蠢的错误。这是一个最小的示例，显示了我使用简单的三角剖分结构（无层次结构）的缓慢解决方案，它确实有效：
#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;
}
如果我转而使用 Triangulation_hierarchy_2 对象，构建在我的 Delaunay_Triangulation_2 类型之上，它不会让我插入包含信息字段的点对 - 它只会如果我自己使用点向量构建对象，则编译：
#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;
}
CGAL 中的所有 typedef 和模板层对我来说都有点不透明 - 有没有什么方法可以设置这些，而不是让我可以基于 Triangulation_hierarchy_2 构建类似 Triangulation_hierarchy_2 的东西=28=]_vertex_base_with_info_2?第二种看起来不存在。
谢谢！
好的 - 我最终尝试了几条路线，并对每条路线进行了一些基本的执行时间基准测试。这些是基于我问题中代码中的设置：
- 具有 25,000 个点的 TIN，主要分布在 XY 中从 0 到 500 的网格上，XY 中有一点噪声，z 值是随机的
- 在 X/Y
中包含 3,000,000 个随机值介于 0 和 500 之间的点的测试云
我尝试了什么：
- CGAL locate() 使用简单的三角结构，任意起始面 --> 79.3 s
- CGAL locate() 使用 Triangulation Hierarchy 结构，任意起始面 --> 4.74 s
- CGAL locate() 使用由最近邻顶点给出的起始面，在输入的顶点云上使用 K-D 树找到 --> 3.41 s
为此，在创建三角剖分后，我遍历面并在输入云顶点索引和三角剖分面句柄之间建立映射：
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;
}
然后我们运行点定位搜索的时候可以这样：
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
}
其中 pcl::Point2DGround 必须是一些自定义点类型，其中表示仅为 2D（以使搜索树正常工作）。
我最终没有尝试对要测试的传入点进行任何类型的排序，因为上面的方法有效并且因为我要定位（）的测试点总是比 TIN 顶点多得多他们自己，所以我认为对他们强加命令可能会更昂贵。
我认为这是针对我的案例的 simplest/fastest 解决方案，所以我将采用它！感谢@marcglisse 和@andreasfabri 的评论。
我有一个 PCL 点云（3D），我想将其转换为地面 TIN（2.5D），然后采样点（2D）以在投影到锡。为此，我一直在使用 CGAL Delaunay 三角剖分 classes，它在大多数情况下都运行良好！
我能够使用基于 Triangulation_vertex_base_with_info_2 构建的 Delaunay_triangulation_2 来实现这一点，并创建一个漂亮的 TIN。我还编写了一个函数，它使用 CGAL locate() 函数提取 2D space 中任意点的面和顶点，这样我就可以在投影到 TIN 上时插入点的高度。我需要信息字段来保存一个索引，让我可以将三角剖分中的顶点关联回 PCL 点云结构中的点。
然而，当使用基本三角剖分时 class locate() 函数很慢（从三角剖分中的任意顶点开始的随机游走），并且因为我必须为要插值的云（以估计预计高度）这是目前我整个管道中最慢的部分。所以我研究了使用三角剖分层次结构 class 来提高效率。
我不知道如何使 Triangulation_hierarchy class 与带有信息的顶点基一起工作，我认为我只是在做一些愚蠢的错误。这是一个最小的示例，显示了我使用简单的三角剖分结构（无层次结构）的缓慢解决方案，它确实有效：
#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;
}
如果我转而使用 Triangulation_hierarchy_2 对象，构建在我的 Delaunay_Triangulation_2 类型之上，它不会让我插入包含信息字段的点对 - 它只会如果我自己使用点向量构建对象，则编译：
#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;
}
CGAL 中的所有 typedef 和模板层对我来说都有点不透明 - 有没有什么方法可以设置这些，而不是让我可以基于 Triangulation_hierarchy_2 构建类似 Triangulation_hierarchy_2 的东西=28=]_vertex_base_with_info_2?第二种看起来不存在。
谢谢！
好的 - 我最终尝试了几条路线，并对每条路线进行了一些基本的执行时间基准测试。这些是基于我问题中代码中的设置：
- 具有 25,000 个点的 TIN，主要分布在 XY 中从 0 到 500 的网格上，XY 中有一点噪声，z 值是随机的
- 在 X/Y 中包含 3,000,000 个随机值介于 0 和 500 之间的点的测试云
我尝试了什么：
- CGAL locate() 使用简单的三角结构，任意起始面 --> 79.3 s
- CGAL locate() 使用 Triangulation Hierarchy 结构，任意起始面 --> 4.74 s
- CGAL locate() 使用由最近邻顶点给出的起始面，在输入的顶点云上使用 K-D 树找到 --> 3.41 s
为此，在创建三角剖分后，我遍历面并在输入云顶点索引和三角剖分面句柄之间建立映射：
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;
}
然后我们运行点定位搜索的时候可以这样：
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
}
其中 pcl::Point2DGround 必须是一些自定义点类型，其中表示仅为 2D（以使搜索树正常工作）。
我最终没有尝试对要测试的传入点进行任何类型的排序，因为上面的方法有效并且因为我要定位（）的测试点总是比 TIN 顶点多得多他们自己，所以我认为对他们强加命令可能会更昂贵。
我认为这是针对我的案例的 simplest/fastest 解决方案，所以我将采用它！感谢@marcglisse 和@andreasfabri 的评论。