点云PCL学习笔记-分割segmentation-RANSAC随机采样一致性算法&&欧式聚类提取

随机采样一致性算法RANSAC

程序实例参考网址：
https://pcl.readthedocs.io/projects/tutorials/en/latest/random_sample_consensus.html?highlight=console#the-explanation
PCLAPI文档：
https://pointclouds.org/documentation/
应用主要是：
1、对点云进行分割，根据设定不同的几何模型（PCL中主要支持的模型有：空间平面、直线、二维或三维圆周、圆球、锥体），估计对应的几何模型参数，在一定允许误差范围内分割出在模型上的点云。
2、点云的配准对的剔除，见点云配准内容

简介

从样本中随即抽选一个样本子集——最小方差估计算法对这个子集计算模型参数——计算所有样本与该模型的偏差——使用一个设定好的阈值与偏差比较——偏差<阈值，样本点属于模型内样本点（inliers）;偏差>阈值，模型外样本点（outliers）——记录下当前inliers的个数，重复这一过程，记录当前最佳模型参数（即inliers个数最多，对应的inliers称为best_ninliers)——每次迭代后计算：期望的误差率、best_ninliers、总样本个数、当前迭代次数，——计算迭代结束评判因子，决定是否结束迭代——结束迭代，最佳模型参数为最终模型参数的估计值
问题：
1、预先设定的阈值，模型如果抽象不好设定
2、RANSAC的迭代次数未知

RANSAC优缺点

RANSAC 的一个优点是它能够对模型参数进行鲁棒估计，即，即使数据集中存在大量异常值，它也可以高度准确地估计参数。 RANSAC 的一个缺点是计算这些参数所需的时间没有上限。当计算的迭代次数有限时，获得的解可能不是最佳的，甚至可能不是很好地拟合数据的解。通过这种方式，RANSAC 提供了一种权衡；通过计算更多的迭代次数，可以增加生成合理模型的概率。 RANSAC 的另一个缺点是它需要设置特定于问题的阈值。

RANSAC 只能为特定数据集估计一个模型。对于存在两个（或更多）模型时的任何一种模型方法，RANSAC 可能无法找到任何一个。

PCL中RANSAC算法实现类

class pcl::RandomSampleConsensus

应用实例

学习如何用RandomSampleConsensus类获得点云的拟合平面模型

初始化点云

  // initialize PointClouds
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);//存储源点云
  pcl::PointCloud<pcl::PointXYZ>::Ptr final (new pcl::PointCloud<pcl::PointXYZ>);//存储提取的局内点inliers

//填充点云数据，作为处理前的原始点云
  cloud->width    = 5000;//设置点云数目
  cloud->height   = 1;//设置为无序点云
  cloud->is_dense = false;
  cloud->points.resize (cloud->width * cloud->height);

for (size_t i = 0; i < cloud->points.size (); ++i)
  {
    if (pcl::console::find_argument (argc, argv, "-s") >= 0 || pcl::console::find_argument (argc, argv, "-sf") >= 0)//如果命令行参数为-s或者-sf，则用x^2+y^2+z^2=1设置一部分点云数据
    {
      cloud->points[i].x =  rand () / (RAND_MAX + 1.0);
      cloud->points[i].y =  rand () / (RAND_MAX + 1.0);
      if (i % 5 == 0)//1/5的点云数据被随机放置作为局外点
        cloud->points[i].z =  rand () / (RAND_MAX + 1.0);
      else if(i % 2 == 0)//另外4/5的点云数据z=+或-(1-x^2-y^2)开根号
        cloud->points[i].z =  sqrt( 1 - (cloud->points[i].x * cloud->points[i].x)
                                      - (cloud->points[i].y * cloud->points[i].y));
      else
        cloud->points[i].z =  - sqrt( 1 - (cloud->points[i].x * cloud->points[i].x)
                                        - (cloud->points[i].y * cloud->points[i].y));
    }
    else//命令行参数未指定，则用x+y+z=1设置一部分点云数据，点云组成的菱形平面作为内点
    {
      cloud->points[i].x =  rand () / (RAND_MAX + 1.0);
      cloud->points[i].y =  rand () / (RAND_MAX + 1.0);
      if( i % 5 == 0)
        cloud->points[i].z = rand () / (RAND_MAX + 1.0);//此处点为局外点
      else
        cloud->points[i].z = -1 * (cloud->points[i].x + cloud->points[i].y);//z=-(x+y)此处为局内点
    }
  }

创建 RandomSampleConsensus 对象并计算合适的模型

 std::vector<int> inliers;//存储局内点集合的点的索引向量

  // created RandomSampleConsensus object and compute the appropriated model
  pcl::SampleConsensusModelSphere<pcl::PointXYZ>::Ptr
    model_s(new pcl::SampleConsensusModelSphere<pcl::PointXYZ> (cloud));//针对球模型的对象
  pcl::SampleConsensusModelPlane<pcl::PointXYZ>::Ptr
    model_p (new pcl::SampleConsensusModelPlane<pcl::PointXYZ> (cloud));//针对平面模型的对象
  if(pcl::console::find_argument (argc, argv, "-f") >= 0)//命令行输入的参数为-f时，随机估算对应的平面模型，并存储估计的局内点
  {
    pcl::RandomSampleConsensus<pcl::PointXYZ> ransac (model_p);//RANSAC，对象为平面
    ransac.setDistanceThreshold (.01);//设置与平面距离小于0.01的点作为局内点考虑
    ransac.computeModel();//执行随机参数估计
    ransac.getInliers(inliers);//存储估计所得的局内点
  }
  else if (pcl::console::find_argument (argc, argv, "-sf") >= 0 )//命令行输入的参数为-sf时，随机估算对应的球面模型，并存储估计的局内点
  {
    pcl::RandomSampleConsensus<pcl::PointXYZ> ransac (model_s);//RANSAC，对象为球面
    ransac.setDistanceThreshold (.01);//与球面距离小于0.01的点作为局内点考虑
    ransac.computeModel();//执行随机参数估计
    ransac.getInliers(inliers);//存储估计所得的局内点
  }

最后，将局内点复制到final点云中，在三维窗口中显示原始点云或者估计得到的局内点组成的点云

 pcl::copyPointCloud<pcl::PointXYZ>(*cloud, inliers, *final);

  // creates the visualization object and adds either our orignial cloud or all of the inliers
  // depending on the command line arguments specified.
  boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer;
  if (pcl::console::find_argument (argc, argv, "-f") >= 0 || pcl::console::find_argument (argc, argv, "-sf") >= 0)
    viewer = simpleVis(final);//命令行参数为-f或者-sf时，显示最终结果
  else
    viewer = simpleVis(cloud);//显示原始点云

segmentation

segmentation平面分割

使用到的头文件

#include <pcl/sample_consensus/method_types.h>//随机参数估计方法头文件
#include <pcl/sample_consensus/model_types.h>//模型定义头文件
#include <pcl/segmentation/sac_segmentation.h>//基于采样一致性分割的类的头文件

预先准备数据

 	pcl::PointCloud<pcl::PointXYZ> cloud;
//填充点云数据
	cloud.width  = 15;
	cloud.height = 1;
	cloud.points.resize (cloud.width * cloud.height);
//生成数据
for (size_t i = 0; i < cloud.points.size (); ++i)
  {
    cloud.points[i].x = 1024 * rand () / (RAND_MAX + 1.0f);
    cloud.points[i].y = 1024 * rand () / (RAND_MAX + 1.0f);
    cloud.points[i].z = 1.0;
  }
//设置几个局外点
cloud.points[0].z = 2.0;
cloud.points[3].z = -2.0;
cloud.points[6].z = 4.0;
//打印点云坐标值
  std::cerr << "Point cloud data: " << cloud.points.size () <<" points" << std::endl;
for (size_t i = 0; i < cloud.points.size (); ++i)
    std::cerr << "    " << cloud.points[i].x << " " 
<< cloud.points[i].y << " " 
<< cloud.points[i].z << std::endl;

分割

//创建分割时所需要的模型系数对象coefficients及存储内点的点索引集合对象inliers
pcl::ModelCoefficients::Ptr coefficients (new pcl::ModelCoefficients);
pcl::PointIndices::Ptr inliers (new pcl::PointIndices);
//创建分割对象
  pcl::SACSegmentation<pcl::PointXYZ> seg;
 //可选设置，设置模型系数需要优化
  seg.setOptimizeCoefficients (true);
//必须设置
	seg.setModelType (pcl::SACMODEL_PLANE);//设置分割的模型类型为平面
  seg.setMethodType (pcl::SAC_RANSAC);//所用参数估计方法为RANSAC
  seg.setDistanceThreshold (0.01);//设置距离阈值小于0.01m为局内点
  seg.setInputCloud (cloud.makeShared ());//输入原始点云
//引发分割实现，并存储分割结果到点集合inliers及存储平面模型系数
  seg.segment (*inliers, *coefficients);
if (inliers->indices.size () == 0)
  {
    PCL_ERROR ("Could not estimate a planar model for the given dataset.");
return (-1);
  }
 //打印平面模型参数ax+by+cz+d=0形式
  std::cerr << "Model coefficients: " << coefficients->values[0] << " " 
<<coefficients->values[1] << " "
<<coefficients->values[2] << " " 
<<coefficients->values[3] <<std::endl;
  std::cerr << "Model inliers: " << inliers->indices.size () << std::endl;
for (size_t i = 0; i < inliers->indices.size (); ++i)
    std::cerr << inliers->indices[i] << "    " <<cloud.points[inliers->indices[i]].x << " "
<<cloud.points[inliers->indices[i]].y << " "
<<cloud.points[inliers->indices[i]].z << std::endl;
return (0);

segmentation 圆柱体分割

任务：
1、估计每个点的表面法线
2、分割出平面模型，存储
3、分割出圆柱体模型，存储

  // All the objects needed
  pcl::PCDReader reader;//PCD文件读取对象
  pcl::PassThrough<PointT> pass;//直通滤波对象
  pcl::NormalEstimation<PointT, pcl::Normal> ne;//法线估计对象
  pcl::SACSegmentationFromNormals<PointT, pcl::Normal> seg;//分割对象 
  pcl::PCDWriter writer;//PCD文件写入对象
  pcl::ExtractIndices<PointT> extract;//点提取对象
  pcl::ExtractIndices<pcl::Normal> extract_normals;
  pcl::search::KdTree<PointT>::Ptr tree (new pcl::search::KdTree<PointT> ());

  // Datasets
  pcl::PointCloud<PointT>::Ptr cloud (new pcl::PointCloud<PointT>);
  pcl::PointCloud<PointT>::Ptr cloud_filtered (new pcl::PointCloud<PointT>);
  pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);
  pcl::PointCloud<PointT>::Ptr cloud_filtered2 (new pcl::PointCloud<PointT>);
  pcl::PointCloud<pcl::Normal>::Ptr cloud_normals2 (new pcl::PointCloud<pcl::Normal>);
  pcl::ModelCoefficients::Ptr coefficients_plane (new pcl::ModelCoefficients), coefficients_cylinder (new pcl::ModelCoefficients);
  pcl::PointIndices::Ptr inliers_plane (new pcl::PointIndices), inliers_cylinder (new pcl::PointIndices);

  // Read in the cloud data
  reader.read ("table_scene_mug_stereo_textured.pcd", *cloud);
  std::cerr << "PointCloud has: " << cloud->points.size () << " data points." << std::endl;

  // Build a passthrough filter to remove spurious NaNs
  pass.setInputCloud (cloud);
  pass.setFilterFieldName ("z");
  pass.setFilterLimits (0, 1.5);
  pass.filter (*cloud_filtered);
  std::cerr << "PointCloud after filtering has: " << cloud_filtered->points.size () << " data points." << std::endl;

  // Estimate point normals
  ne.setSearchMethod (tree);
  ne.setInputCloud (cloud_filtered);
  ne.setKSearch (50);
  ne.compute (*cloud_normals);
 // Create the segmentation object for the planar model and set all the parameters
  seg.setOptimizeCoefficients (true);
  seg.setModelType (pcl::SACMODEL_NORMAL_PLANE);//设置分割模型为平面
  seg.setNormalDistanceWeight (0.1);
  seg.setMethodType (pcl::SAC_RANSAC);
  seg.setMaxIterations (100);
  seg.setDistanceThreshold (0.03);
  seg.setInputCloud (cloud_filtered);
  seg.setInputNormals (cloud_normals);
  // Obtain the plane inliers and coefficients
  seg.segment (*inliers_plane, *coefficients_plane);
  std::cerr << "Plane coefficients: " << *coefficients_plane << std::endl;  

抽取分割点云并保存

  // Extract the planar inliers from the input cloud
  extract.setInputCloud (cloud_filtered);
  extract.setIndices (inliers_plane);
  extract.setNegative (false);
 // Write the planar inliers to disk
  pcl::PointCloud<PointT>::Ptr cloud_plane (new pcl::PointCloud<PointT> ());
  extract.filter (*cloud_plane);
  std::cerr << "PointCloud representing the planar component: " << cloud_plane->points.size () << " data points." << std::endl;
  writer.write ("table_scene_mug_stereo_textured_plane.pcd", *cloud_plane, false);

  // Remove the planar inliers, extract the rest
  extract.setNegative (true);
  extract.filter (*cloud_filtered2);
  extract_normals.setNegative (true);
  extract_normals.setInputCloud (cloud_normals);
  extract_normals.setIndices (inliers_plane);
  extract_normals.filter (*cloud_normals2);

ExtractIndices基于某一分割算法提取点云中的一个子集

class pcl::ExtractIndices< PointT >
例程详见：
https://pointclouds.org/documentation/classpcl_1_1_extract_indices.html#details

pcl::ExtractIndices<PointType> eifilter (true); // Initializing with true will allow us to extract the removed indices
eifilter.setInputCloud (cloud_in);
eifilter.setIndices (indices_in);
eifilter.filter (*cloud_out);
// The resulting cloud_out contains all points of cloud_in that are indexed by indices_in
indices_rem = eifilter.getRemovedIndices ();
// The indices_rem array indexes all points of cloud_in that are not indexed by indices_in
eifilter.setNegative (true);
eifilter.filter (*indices_out);
// Alternatively: the indices_out array is identical to indices_rem
eifilter.setNegative (false);
eifilter.setUserFilterValue (1337.0);
eifilter.filterDirectly (cloud_in);
// This will directly modify cloud_in instead of creating a copy of the cloud
// It will overwrite all fields of the filtered points by the user value: 1337

欧式聚类

概念

聚类是指把不同物体的点云分别组合聚集起来, 从而能让你跟踪汽车, 行人等多个目标. 其中一种对点云数据进行分组和聚类的方法称为欧氏聚类.

欧式聚类是指将距离紧密度高的点云聚合起来. 为了有效地进行最近邻搜索, 可以使用 KD-Tree 数据结构, 这种结构平均可以加快从 o (n)到 o (log (n))的查找时间. 这是因为Kd-Tree允许你更好地分割你的搜索空间. 通过将点分组到 KD-Tree 中的区域中, 您可以避免计算可能有数千个点的距离, 因为你知道它们不会被考虑在一个紧邻的区域中.

欧氏聚类对象 ec 具有距离容忍度. 在这个距离之内的任何点都将被组合在一起. 它还有用于表示集群的点数的 min 和 max 参数. 如果一个集群真的很小, 它可能只是噪音, min参数会限制使用这个集群. 而max参数允许我们更好地分解非常大的集群, 如果一个集群非常大, 可能只是许多其他集群重叠, 最大容差可以帮助我们更好地解决对象检测. 欧式聚类的最后一个参数是 Kd-Tree. Kd-Tree是使用输入点云创建和构建的, 这些输入云点是初始点云过滤分割后得到障碍物点云.

实例

pcl::EuclideanClusterExtractionpcl::PointXYZ类
采用欧式聚类对三维点云进行分割
预先的步骤是：读入一个PCD文件，进行滤波处理，然后RANSAC平面模型分割处理，提取出点云所有在平面上的点集，存储

  // Read in the cloud data
  pcl::PCDReader reader;
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>), cloud_f (new pcl::PointCloud<pcl::PointXYZ>);
  reader.read ("table_scene_lms400.pcd", *cloud);
  std::cout << "PointCloud before filtering has: " << cloud->points.size () << " data points." << std::endl; //*

  // Create the filtering object: downsample the dataset using a leaf size of 1cm
  pcl::VoxelGrid<pcl::PointXYZ> vg;
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered (new pcl::PointCloud<pcl::PointXYZ>);
  vg.setInputCloud (cloud);
  vg.setLeafSize (0.01f, 0.01f, 0.01f);
  vg.filter (*cloud_filtered);
  std::cout << "PointCloud after filtering has: " << cloud_filtered->points.size ()  << " data points." << std::endl; //*

  // Create the segmentation object for the planar model and set all the parameters
  pcl::SACSegmentation<pcl::PointXYZ> seg;
  pcl::PointIndices::Ptr inliers (new pcl::PointIndices);
  pcl::ModelCoefficients::Ptr coefficients (new pcl::ModelCoefficients);
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_plane (new pcl::PointCloud<pcl::PointXYZ> ());
  pcl::PCDWriter writer;
  seg.setOptimizeCoefficients (true);
  seg.setModelType (pcl::SACMODEL_PLANE);
  seg.setMethodType (pcl::SAC_RANSAC);
  seg.setMaxIterations (100);
  seg.setDistanceThreshold (0.02);

  int i=0, nr_points = (int) cloud_filtered->points.size ();
  while (cloud_filtered->points.size () > 0.3 * nr_points)
  {
    // Segment the largest planar component from the remaining cloud
    seg.setInputCloud (cloud_filtered);
    seg.segment (*inliers, *coefficients);
    if (inliers->indices.size () == 0)
    {
      std::cout << "Could not estimate a planar model for the given dataset." << std::endl;
      break;
    }

    // Extract the planar inliers from the input cloud
    pcl::ExtractIndices<pcl::PointXYZ> extract;
    extract.setInputCloud (cloud_filtered);
    extract.setIndices (inliers);
    extract.setNegative (false);

    // Write the planar inliers to disk
    extract.filter (*cloud_plane);
    std::cout << "PointCloud representing the planar component: " << cloud_plane->points.size () << " data points." << std::endl;

    // Remove the planar inliers, extract the rest
    extract.setNegative (true);
    extract.filter (*cloud_f);
    cloud_filtered = cloud_f;
  }

  // Creating the KdTree object for the search method of the extraction创建一个kd树作为提取点云时的搜索方法，再创建一个点云索引向量，用于存储实际的点云索引信息，每个检测到的点云聚类被保存在这里，如cluster_indices[0]包含点云中第一个聚类包含的点集的所有索引
  pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);
  tree->setInputCloud (cloud_filtered);

  std::vector<pcl::PointIndices> cluster_indices;//创建点云索引向量
  pcl::EuclideanClusterExtraction<pcl::PointXYZ> ec;
  ec.setClusterTolerance (0.02); // 设置近邻搜索半径为2cm
  ec.setMinClusterSize (100);//设置一个聚类需要的最少点数目为100
  ec.setMaxClusterSize (25000);//设置一个聚类需要的最大点数目为25000
  ec.setSearchMethod (tree);//设置点云的搜索机制
  ec.setInputCloud (cloud_filtered);
  ec.extract (cluster_indices);//从点云中提取聚类，并将点云索引保存在cluster_indices

ClusterTolerance 搜索半径需要通过测试来设置一个合适的聚类搜索半径

	//迭代访问点云索引cluster_indices,直到分割出所有聚类
	int j=0;
  for (std::vector<pcl::PointIndices>::const_iterator it = cluster_indices.begin (); it != cluster_indices.end (); ++it)
  {
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_cluster (new pcl::PointCloud<pcl::PointXYZ>);//创建新的点云数据集cloud_cluster,将所有当前聚类写入到点云数据集中
    for (std::vector<int>::const_iterator pit = it->indices.begin (); pit != it->indices.end (); pit++)
      cloud_cluster->points.push_back (cloud_filtered->points[*pit]); //*
    cloud_cluster->width = cloud_cluster->points.size ();
    cloud_cluster->height = 1;
    cloud_cluster->is_dense = true;

    std::cout << "PointCloud representing the Cluster: " << cloud_cluster->points.size () << " data points." << std::endl;
    std::stringstream ss;
    ss << "cloud_cluster_" << j << ".pcd";
    writer.write<pcl::PointXYZ> (ss.str (), *cloud_cluster, false); //*
    j++;
  }