https://pcl.readthedocs.io/projects/tutorials/en/latest/#surface

https://pointcloudlibrary.github.io/documentation/group__surface.html

此模块个人代码：https://github.com/HuangCongQing/pcl-learning/tree/master/11surface

博客

http://robot.czxy.com/docs/pcl/chapter04/resampling/
pcl几种表面重建_3D的博客-CSDN博客

PCL点云曲面重建（1）
PCL法线估计

基本概念

曲面重建技术在逆向工程、数据可视化、机器视觉、虚拟现实、医疗技术等领域中得到了广泛的应用 。 例如，在汽车、航空等工业领域中，复杂外形产品的设计仍需要根据手工模型，采用逆向工程的手段建立产品的数字化模型，根据测量数据建立人体以及骨骼和器官的计算机模型，在医学、定制生产等方面都有重要意义 。
除了上述传统的行业，随着新兴的廉价 RGBD 获取设备在数字娱乐行业的病毒式扩展，使得更多人开始使用点云来处理对象并进行工程应用 。
根据重建曲面和数据点云之间的关系，可将曲面重建分为两大类：插值法和逼近法

  • 插值法得到的重建曲面完全通过原始数据点
  • 逼近法是用分片线性曲面或其他形式的曲面来逼近原始数据点,从而使得得到的重建曲面是原始点集的一个逼近曲面。

代码实践

1 基于多项式重构的平滑和法线估计

本教程说明如何使用移动最小二乘（MLS）曲面重构方法来平滑和重采样噪声数据。
使用统计分析很难消除某些数据不规则性（由较小的距离测量误差引起）。要创建完整的模型，必须考虑光滑的表面以及数据中的遮挡。在无法获取其他扫描的情况下，一种解决方案是使用重采样算法，该算法尝试通过周围数据点之间的高阶多项式插值来重新创建表面的缺失部分。通过执行重采样，可以纠正这些小的错误，并且可以将多个扫描记录在一起执行平滑操作合并成同一个点云。
11surface表面 - 图1

在上图的左侧，我们在包含两个配准点云的数据集中看到了配准后的效果及表面法线估计。由于对齐错误，所产生的法线有噪声。在右侧，使用移动最小二乘法对表面法线估计进行平滑处理后，在同一数据集中看到了该效果。绘制每个点的曲率，作为重新采样前后特征值关系的度量，我们得到：

11surface表面 - 图2

函数对象

class pcl::MovingLeastSquares< PointInT, PointOutT >
MovingLeastSquares represent an implementation of the MLS (Moving Least Squares) algorithm for data smoothing and improved normal estimation. More…

代码实现

创建文件：resampling.cpp
准备资源：../ism_train_cat.pcd
编译执行：./resampling

个人代码：11surface表面

  1. /*
  2.  * @Description: 基于多项式重构的平滑和法线估计¶
  3. http://robot.czxy.com/docs/pcl/chapter04/resampling/#_2   推荐
  4. https://www.cnblogs.com/li-yao7758258/p/6497446.html
  5.  * @Author: HCQ
  6.  * @Company(School): UCAS
  7.  * @Email: 1756260160@qq.com
  8.  * @Date: 2020-10-21 18:08:55
  9.  * @LastEditTime: 2020-10-21 20:15:10
  10.  * @FilePath: /pcl-learning/11surface表面 /1基于多项式重构的平滑和法线估计/resampling.cpp
  11.  */
  12. #include <pcl/point_types.h>
  13. #include <pcl/io/pcd_io.h>
  14. #include <pcl/kdtree/kdtree_flann.h>
  15. #include <pcl/surface/mls.h>
  16. #include <pcl/visualization/cloud_viewer.h>
  18. int
  19. main(int argc, char **argv) {
  20.     // Load input file into a PointCloud<T> with an appropriate type
  21.     pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>());
  22.     // Load bun0.pcd -- should be available with the PCL archive in test
  23.     // pcl::io::loadPCDFile(argv[1], *cloud);  // 获取pcd文件
  24.     pcl::io::loadPCDFile("../ism_train_cat.pcd", *cloud);  // 获取pcd文件
  26.     // Create a KD-Tree
  27.     pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>);
  29.     // Output has the PointNormal type in order to store the normals calculated by MLS
  30.     pcl::PointCloud<pcl::PointNormal> mls_points;
  32.     // Init object (second point type is for the normals, even if unused)
  33.     pcl::MovingLeastSquares<pcl::PointXYZ, pcl::PointNormal> mls;
  35.     mls.setComputeNormals(true);
  37.     // Set parameters
  38.     mls.setInputCloud(cloud);
  39.     mls.setPolynomialOrder(2);
  40.     mls.setSearchMethod(tree);
  41.     mls.setSearchRadius(0.03);
  43.     // Reconstruct
  44.     mls.process(mls_points);  // void     process (PointCloudOut &output) override
  45.     // Save output
  46.     if(mls_points.size() > 0){
  47.         pcl::io::savePCDFileASCII("../target-mls.pcd", mls_points); //保存出错
  48.     }else{
  49.             std::cout << "保存数据为空." << std::endl;
  50.     }
  52.     pcl::visualization::CloudViewer viewer("Cloud Viewer");;
  53.     viewer.showCloud(cloud);
  54.     while (!viewer.wasStopped()) {
  55.     }
  57. }

实现效果

2 在平面模型上提取凸（凹）多边形

本例子

  1. 先从点云中提取平面模型
  2. 再通过该估计的平面模型系数从滤波后的点云投影一组点集形成点云
  3. 最后为投影后的点云计算其对应的二维凸多边形

代码实现

创建文件：hull_2d.cpp
准备资源：../table_scene_mug_stereo_textured.pcd
编译执行：./hull_2d

个人代码：11surface表面

/*
 * @Description: 在平面模型上提取凸（凹）多边形     https://www.cnblogs.com/li-yao7758258/p/6497446.html
 * @Author: HCQ
 * @Company(School): UCAS
 * @Email: 1756260160@qq.com
 * @Date: 2020-10-22 10:06:12
 * @LastEditTime: 2020-10-22 10:15:10
 * @FilePath: /pcl-learning/11surface表面 /2在平面模型上提取凸（凹）多边形/hull_2d.cpp
 */
#include <pcl/ModelCoefficients.h>           //采样一致性模型相关类头文件
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/sample_consensus/method_types.h>
#include <pcl/sample_consensus/model_types.h>
#include <pcl/filters/passthrough.h>
#include <pcl/filters/project_inliers.h>          //滤波相关类头文件
#include <pcl/segmentation/sac_segmentation.h>   //基于采样一致性分割类定义的头文件
#include <pcl/surface/concave_hull.h>                 //创建凹多边形类定义头文件

int
main (int argc, char** argv)
{
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>), 
                                      cloud_filtered (new pcl::PointCloud<pcl::PointXYZ>), 
                                      cloud_projected (new pcl::PointCloud<pcl::PointXYZ>);
  pcl::PCDReader reader;

  reader.read ("../table_scene_mug_stereo_textured.pcd", *cloud);
  // 建立过滤器消除杂散的NaN
  pcl::PassThrough<pcl::PointXYZ> pass;
  pass.setInputCloud (cloud);                  //设置输入点云
  pass.setFilterFieldName ("z");             //设置分割字段为z坐标
  pass.setFilterLimits (0, 1.1);             //设置分割阀值为(0, 1.1)
  pass.filter (*cloud_filtered);              
  std::cerr << "PointCloud after filtering has: "
            << cloud_filtered->points.size () << " data points." << std::endl;

// 分割
  pcl::ModelCoefficients::Ptr coefficients (new pcl::ModelCoefficients);
  pcl::PointIndices::Ptr inliers (new pcl::PointIndices);   //inliers存储分割后的点云
  // 创建分割对象
  pcl::SACSegmentation<pcl::PointXYZ> seg;
  // 设置优化系数，该参数为可选参数
  seg.setOptimizeCoefficients (true);
  // Mandatory
  seg.setModelType (pcl::SACMODEL_PLANE);
  seg.setMethodType (pcl::SAC_RANSAC);
  seg.setDistanceThreshold (0.01);

  seg.setInputCloud (cloud_filtered);
  seg.segment (*inliers, *coefficients);
  std::cerr << "PointCloud after segmentation has: "
            << inliers->indices.size () << " inliers." << std::endl;

  // Project the model inliers点云投影滤波模型
  pcl::ProjectInliers<pcl::PointXYZ> proj;//点云投影滤波模型
  proj.setModelType (pcl::SACMODEL_PLANE); //设置投影模型
  proj.setIndices (inliers);             
  proj.setInputCloud (cloud_filtered);
  proj.setModelCoefficients (coefficients);      //将估计得到的平面coefficients参数设置为投影平面模型系数
  proj.filter (*cloud_projected);            //得到投影后的点云
  std::cerr << "PointCloud after projection has: "
            << cloud_projected->points.size () << " data points." << std::endl;

  // 存储提取多边形上的点云
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_hull (new pcl::PointCloud<pcl::PointXYZ>);
  pcl::ConcaveHull<pcl::PointXYZ> chull;        //创建多边形提取对象
  chull.setInputCloud (cloud_projected);       //设置输入点云为提取后点云
  chull.setAlpha (0.1);
  chull.reconstruct (*cloud_hull);   //创建提取创建凹多边形

  std::cerr << "Concave hull has: " << cloud_hull->points.size ()
            << " data points." << std::endl;

  pcl::PCDWriter writer;
  writer.write ("../table_scene_mug_stereo_textured_hull.pcd", *cloud_hull, false);

  return (0);
}

实现效果

image.png

原始 处理后（红边）
image.png image.png

3 无序点云的快速三角化

使用贪婪投影三角化算法对有向点云进行三角化，
具体方法是：
（1）先将有向点云投影到某一局部二维坐标平面内
（2）在坐标平面内进行平面内的三角化
（3）根据平面内三位点的拓扑连接关系获得一个三角网格曲面模型.
贪婪投影三角化算法原理：
是处理一系列可以使网格“生长扩大”的点（边缘点）延伸这些点直到所有符合几何正确性和拓扑正确性的点都被连上，该算法可以用来处理来自一个或者多个扫描仪扫描到得到并且有多个连接处的散乱点云但是算法也是有很大的局限性，它更适用于采样点云来自表面连续光滑的曲面且点云的密度变化比较均匀的情况

代码实现

创建文件：greedy_projection.cpp
准备资源：../bun0.pcd
编译执行：./greedy_projection

个人代码：11surface表面

/*
 * @Description: 无序点云的快速三角化  https://www.cnblogs.com/li-yao7758258/p/6497446.html
 * @Author: HCQ
 * @Company(School): UCAS
 * @Email: 1756260160@qq.com
 * @Date: 2020-10-22 10:35:16
 * @LastEditTime: 2020-10-22 10:50:31
 * @FilePath: /pcl-learning/11surface表面 /3无序点云的快速三角化/greedy_projection.cpp
 */
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/features/normal_3d.h>
#include <pcl/surface/gp3.h>      //贪婪投影三角化算法
int
main (int argc, char** argv)
{
  // 将一个XYZ点类型的PCD文件打开并存储到对象中
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
  pcl::PCLPointCloud2 cloud_blob;
  pcl::io::loadPCDFile ("../bun0.pcd", cloud_blob);
  pcl::fromPCLPointCloud2 (cloud_blob, *cloud);
  //* the data should be available in cloud

  // Normal estimation*
  pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> n;      //法线估计对象
  pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal>);   //存储估计的法线
  pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);  //定义kd树指针
  tree->setInputCloud (cloud); //用cloud构建tree对象
  n.setInputCloud (cloud);
  n.setSearchMethod (tree);
  n.setKSearch (20);
  n.compute (*normals); //估计法线存储到其中
  //* normals should not contain the point normals + surface curvatures

  // Concatenate the XYZ and normal fields*
  pcl::PointCloud<pcl::PointNormal>::Ptr cloud_with_normals (new pcl::PointCloud<pcl::PointNormal>);
  pcl::concatenateFields (*cloud, *normals, *cloud_with_normals);    //连接字段
  //* cloud_with_normals = cloud + normals

  //定义搜索树对象
  pcl::search::KdTree<pcl::PointNormal>::Ptr tree2 (new pcl::search::KdTree<pcl::PointNormal>);
  tree2->setInputCloud (cloud_with_normals);   //点云构建搜索树

  // Initialize objects
  pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp3;   //定义三角化对象
  pcl::PolygonMesh triangles;                //存储最终三角化的网络模型

  // Set the maximum distance between connected points (maximum edge length)
  gp3.setSearchRadius (0.025);  //设置连接点之间的最大距离，（即是三角形最大边长）

  // 设置各参数值
  gp3.setMu (2.5);  //设置被样本点搜索其近邻点的最远距离为2.5，为了使用点云密度的变化
  gp3.setMaximumNearestNeighbors (100);    //设置样本点可搜索的邻域个数
  gp3.setMaximumSurfaceAngle(M_PI/4); // 设置某点法线方向偏离样本点法线的最大角度45
  gp3.setMinimumAngle(M_PI/18); // 设置三角化后得到的三角形内角的最小的角度为10
  gp3.setMaximumAngle(2*M_PI/3); // 设置三角化后得到的三角形内角的最大角度为120
  gp3.setNormalConsistency(false);  //设置该参数保证法线朝向一致

  // Get result
  gp3.setInputCloud (cloud_with_normals);     //设置输入点云为有向点云
  gp3.setSearchMethod (tree2);   //设置搜索方式
  gp3.reconstruct (triangles);  //重建提取三角化

  // 附加顶点信息
  std::vector<int> parts = gp3.getPartIDs();
  std::vector<int> states = gp3.getPointStates();


  // Finish
  return (0);
}

实现效果

image.pngimage.png

4 将修剪的B样条曲线拟合到无序点云

pcl/doc/tutorials/content/sources/bspline_fitting/bspline_fitting.cpp

代码实现

创建文件：bspline_fitting.cpp
准备资源：../bunny.pcd
编译执行：./bspline_fitting ../bunny.pcd

个人代码：11surface表面

/*
 * @Description: 将修剪的B样条曲线拟合到无序点云 https://pcl.readthedocs.io/projects/tutorials/en/latest/bspline_fitting.html#compiling-and-running-the-program
 * @Author: HCQ
 * @Company(School): UCAS
 * @Email: 1756260160@qq.com
 * @Date: 2020-10-22 11:08:24
 * @LastEditTime: 2020-10-22 11:35:49
 * @FilePath: /pcl-learning/11surface表面 /4将修剪的B样条曲线拟合到无序点云/bspline_fitting.cpp
 */

#include <pcl/point_cloud.h>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>

#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/surface/on_nurbs/fitting_surface_tdm.h> // fatal error: pcl/surface/on_nurbs/fitting_surface_tdm.h: 没有那个文件或目录
#include <pcl/surface/on_nurbs/fitting_curve_2d_asdm.h>
#include <pcl/surface/on_nurbs/triangulation.h>

typedef pcl::PointXYZ Point;

void
PointCloud2Vector3d (pcl::PointCloud<Point>::Ptr cloud, pcl::on_nurbs::vector_vec3d &data);

void
visualizeCurve (ON_NurbsCurve &curve,
                ON_NurbsSurface &surface,
                pcl::visualization::PCLVisualizer &viewer);

int
main (int argc, char *argv[])
{
  std::string pcd_file, file_3dm;

  if (argc < 3)
  {
    printf ("\nUsage: pcl_example_nurbs_fitting_surface pcd<PointXYZ>-in-file 3dm-out-file\n\n");
    exit (0);
  }
  pcd_file = argv[1];
  file_3dm = argv[2];

  pcl::visualization::PCLVisualizer viewer ("B-spline surface fitting");
  viewer.setSize (800, 600);

  // ############################################################################
  // load point cloud

  printf ("  loading %s\n", pcd_file.c_str ());
  pcl::PointCloud<Point>::Ptr cloud (new pcl::PointCloud<Point>);
  pcl::PCLPointCloud2 cloud2;
  pcl::on_nurbs::NurbsDataSurface data;

  if (pcl::io::loadPCDFile (pcd_file, cloud2) == -1)
    throw std::runtime_error ("  PCD file not found.");

  fromPCLPointCloud2 (cloud2, *cloud);
  PointCloud2Vector3d (cloud, data.interior);
  pcl::visualization::PointCloudColorHandlerCustom<Point> handler (cloud, 0, 255, 0);
  viewer.addPointCloud<Point> (cloud, handler, "cloud_cylinder");
  printf ("  %lu points in data set\n", cloud->size ());

  // ############################################################################
  // fit B-spline surface

  // parameters
  unsigned order (3);
  unsigned refinement (5);
  unsigned iterations (10);
  unsigned mesh_resolution (256);

  pcl::on_nurbs::FittingSurface::Parameter params;
  params.interior_smoothness = 0.2;
  params.interior_weight = 1.0;
  params.boundary_smoothness = 0.2;
  params.boundary_weight = 0.0;

  // initialize
  printf ("  surface fitting ...\n");
  ON_NurbsSurface nurbs = pcl::on_nurbs::FittingSurface::initNurbsPCABoundingBox (order, &data);
  pcl::on_nurbs::FittingSurface fit (&data, nurbs);
  //  fit.setQuiet (false); // enable/disable debug output

  // mesh for visualization
  pcl::PolygonMesh mesh;
  pcl::PointCloud<pcl::PointXYZ>::Ptr mesh_cloud (new pcl::PointCloud<pcl::PointXYZ>);
  std::vector<pcl::Vertices> mesh_vertices;
  std::string mesh_id = "mesh_nurbs";
  pcl::on_nurbs::Triangulation::convertSurface2PolygonMesh (fit.m_nurbs, mesh, mesh_resolution);
  viewer.addPolygonMesh (mesh, mesh_id);

  // surface refinement
  for (unsigned i = 0; i < refinement; i++)
  {
    fit.refine (0);
    fit.refine (1);
    fit.assemble (params);
    fit.solve ();
    pcl::on_nurbs::Triangulation::convertSurface2Vertices (fit.m_nurbs, mesh_cloud, mesh_vertices, mesh_resolution);
    viewer.updatePolygonMesh<pcl::PointXYZ> (mesh_cloud, mesh_vertices, mesh_id);
    viewer.spinOnce ();
  }

  // surface fitting with final refinement level
  for (unsigned i = 0; i < iterations; i++)
  {
    fit.assemble (params);
    fit.solve ();
    pcl::on_nurbs::Triangulation::convertSurface2Vertices (fit.m_nurbs, mesh_cloud, mesh_vertices, mesh_resolution);
    viewer.updatePolygonMesh<pcl::PointXYZ> (mesh_cloud, mesh_vertices, mesh_id);
    viewer.spinOnce ();
  }

  // ############################################################################
  // fit B-spline curve

  // parameters
  pcl::on_nurbs::FittingCurve2dAPDM::FitParameter curve_params;
  curve_params.addCPsAccuracy = 5e-2;
  curve_params.addCPsIteration = 3;
  curve_params.maxCPs = 200;
  curve_params.accuracy = 1e-3;
  curve_params.iterations = 100;

  curve_params.param.closest_point_resolution = 0;
  curve_params.param.closest_point_weight = 1.0;
  curve_params.param.closest_point_sigma2 = 0.1;
  curve_params.param.interior_sigma2 = 0.00001;
  curve_params.param.smooth_concavity = 1.0;
  curve_params.param.smoothness = 1.0;

  // initialisation (circular)
  printf ("  curve fitting ...\n");
  pcl::on_nurbs::NurbsDataCurve2d curve_data;
  curve_data.interior = data.interior_param;
  curve_data.interior_weight_function.push_back (true);
  ON_NurbsCurve curve_nurbs = pcl::on_nurbs::FittingCurve2dAPDM::initNurbsCurve2D (order, curve_data.interior);

  // curve fitting
  pcl::on_nurbs::FittingCurve2dASDM curve_fit (&curve_data, curve_nurbs);
  // curve_fit.setQuiet (false); // enable/disable debug output
  curve_fit.fitting (curve_params);
  visualizeCurve (curve_fit.m_nurbs, fit.m_nurbs, viewer);

  // ############################################################################
  // triangulation of trimmed surface

  printf ("  triangulate trimmed surface ...\n");
  viewer.removePolygonMesh (mesh_id);
  pcl::on_nurbs::Triangulation::convertTrimmedSurface2PolygonMesh (fit.m_nurbs, curve_fit.m_nurbs, mesh,
                                                                   mesh_resolution);
  viewer.addPolygonMesh (mesh, mesh_id);


  // save trimmed B-spline surface
  if ( fit.m_nurbs.IsValid() )
  {
    ONX_Model model;
    ONX_Model_Object& surf = model.m_object_table.AppendNew();
    surf.m_object = new ON_NurbsSurface(fit.m_nurbs);
    surf.m_bDeleteObject = true;
    surf.m_attributes.m_layer_index = 1;
    surf.m_attributes.m_name = "surface";

    ONX_Model_Object& curv = model.m_object_table.AppendNew();
    curv.m_object = new ON_NurbsCurve(curve_fit.m_nurbs);
    curv.m_bDeleteObject = true;
    curv.m_attributes.m_layer_index = 2;
    curv.m_attributes.m_name = "trimming curve";

    model.Write(file_3dm.c_str());
    printf("  model saved: %s\n", file_3dm.c_str());
  }

  printf ("  ... done.\n");

  viewer.spin ();
  return 0;
}

void
PointCloud2Vector3d (pcl::PointCloud<Point>::Ptr cloud, pcl::on_nurbs::vector_vec3d &data)
{
  for (unsigned i = 0; i < cloud->size (); i++)
  {
    Point &p = cloud->at (i);
    if (!std::isnan (p.x) && !std::isnan (p.y) && !std::isnan (p.z))
      data.push_back (Eigen::Vector3d (p.x, p.y, p.z));
  }
}

void
visualizeCurve (ON_NurbsCurve &curve, ON_NurbsSurface &surface, pcl::visualization::PCLVisualizer &viewer)
{
  pcl::PointCloud<pcl::PointXYZRGB>::Ptr curve_cloud (new pcl::PointCloud<pcl::PointXYZRGB>);

  pcl::on_nurbs::Triangulation::convertCurve2PointCloud (curve, surface, curve_cloud, 4);
  for (std::size_t i = 0; i < curve_cloud->size () - 1; i++)
  {
    pcl::PointXYZRGB &p1 = curve_cloud->at (i);
    pcl::PointXYZRGB &p2 = curve_cloud->at (i + 1);
    std::ostringstream os;
    os << "line" << i;
    viewer.removeShape (os.str ());
    viewer.addLine<pcl::PointXYZRGB> (p1, p2, 1.0, 0.0, 0.0, os.str ());
  }

  pcl::PointCloud<pcl::PointXYZRGB>::Ptr curve_cps (new pcl::PointCloud<pcl::PointXYZRGB>);
  for (int i = 0; i < curve.CVCount (); i++)
  {
    ON_3dPoint p1;
    curve.GetCV (i, p1);

    double pnt[3];
    surface.Evaluate (p1.x, p1.y, 0, 3, pnt);
    pcl::PointXYZRGB p2;
    p2.x = float (pnt[0]);
    p2.y = float (pnt[1]);
    p2.z = float (pnt[2]);

    p2.r = 255;
    p2.g = 0;
    p2.b = 0;

    curve_cps->push_back (p2);
  }
  viewer.removePointCloud ("cloud_cps");
  viewer.addPointCloud (curve_cps, "cloud_cps");
}

实现效果 ERROR

ERROR: fatal error: pcl/memory.h: 没有那个文件或目录
image.png