【PCL自学：Recognition 2】隐式形状模型识别原理及示例代码

  在本章节中，我们将学习如何使用在pcl::ism::ImplicitShapeModel类中实现的隐式形状模型算法。Jan Knopp、Mukta Prasad、Geert Willems、Radu Timofte和Luc Van Gool在文章“ Hough Transforms and 3D SURF for robust three dimensional classification” F中描述了这种算法，用于鲁棒的三维分类。该算法是广义Hough变换和特征集方法的结合，其目的如下。假设有一堆已知来自不同物体的点云训练集，算法会为这些点云计算一个特定的模型，然后这些模型将被用于在一些没有训练过的点云中进行物体识别并预测物体中心。

理论基础

  该算法由两个步骤组成，第一步是训练，第二步是识别训练集之外的点云中物体。让我们来看看训练是如何完成的。它包括六个步骤：

【1】首先进行了关键点检测。在给定的实现中，它只是训练云的简化。在这一步中，所有的点云都通过体素网格方法进行简化;其余的点被声明为关键点。
【2】对每个关键点进行特征估计。例如FPFH估计。
【3】利用k-means算法对所有特征进行聚类，构建视觉词(或几何词)典。获得的聚类表示某物体的视觉词。聚类中的每个特性都是这个可视化词的实例。
【4】对于每一个单独的实例，计算到中心的方向——从关键点(获得特征的地方)到给定云的质心的方向。
【5】对于每个视觉单词，由公式计算统计权重

  统计权重$W_{st}(c_i,v_j)$对视觉单词$v_j$为类$c_i$投出的所有选票进行加权。这里$n_{vot}(v_j)$是视觉单词$v_j$的投票总数，$n_{vot}(c_i,v_j)$是$v_j$对类$c_i$的投票数，$n_{vw}(c_i)$是对类$c_i$投票的视觉单词的数量，$n_{ftr}(c_i)$是$c_i$从中学习到的特征的数量。$C$是所有类的集合。
【6】对于每个关键点(估计特征的点)，根据公式计算学习到的权值

  论文作者将${\lambda }_{ij}$定义为视觉单词$v_j$的一个特定实例对类$c_i$的一个特定训练形状的投票;也就是说，${\lambda }_{ij}$记录了视觉单词$v_j$的特定实例到其所在的训练形状中心的距离。这里$A$是与单词$v_j$有关的所有特征的集合，集合的形状为类$c_i$。$\sigma$的推荐值是形状尺寸的10%。函数f只是一个中值。$d_a({\lambda }_{ij})$是投票中心与实际中心之间的欧氏距离。

  训练过程完成，得到训练的模型(权值、方向等)后，进行目标搜索(或识别)过程。它由接下来的四个步骤组成:
【1】特征点提取
【2】点云中每个关键点的特征估计。
【3】对于每个特征，在字典中搜索最近的视觉词(即一个聚类)。
【4】对每个特征做出如下操作：
  [4.1]为训练过的模型中每个视觉单词的每个实例投票(为感兴趣的类投票)
  [4.2]将根据公式计算出的对应方向的选票和票数相加

【5】前面的步骤为我们计算了预测的中心和每一票的权重。为了得到与中心相对应的单点，需要对这些选票进行分析。为此，算法采用非极大值抑制方法。用户只需要传递感兴趣的对象的半径，其余的将由ISMVoteList::findStrongestPeaks()方法完成。

更为详细的理论过程，请阅读这篇论文。

示例代码及注释

  首先，将本章需要的点云集-训练集和识别云集下载到合适的位置。下面是非常适合本教程的云的列表。

【训练集点云】：Cat (train) Horse (train) Lioness (train) Michael (train) Wolf (train)

【测试集点云】：Cat Horse Lioness Lioness Michael Wolf

  接下来需要做的是在您喜欢的任何编辑器中创建一个文件implicit_shape_model.cpp，并将以下代码复制到其中

 #include <iostream>
 #include <pcl/io/pcd_io.h>
 #include <pcl/features/normal_3d.h>
 #include <pcl/features/feature.h>
 #include <pcl/visualization/cloud_viewer.h>
 #include <pcl/features/fpfh.h>
 #include <pcl/features/impl/fpfh.hpp>
 #include <pcl/recognition/implicit_shape_model.h>
 #include <pcl/recognition/impl/implicit_shape_model.hpp>
 
 int
 main (int argc, char** argv)
 {
   if (argc == 0 || argc % 2 == 0)
     return (-1);
 
   unsigned int number_of_training_clouds = (argc - 3) / 2;
 
   pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;
   normal_estimator.setRadiusSearch (25.0);
 
   std::vector<pcl::PointCloud<pcl::PointXYZ>::Ptr> training_clouds;
   std::vector<pcl::PointCloud<pcl::Normal>::Ptr> training_normals;
   std::vector<unsigned int> training_classes;
 
// =====================================================================
//	以下只是加载将用于训练的云。算法需要法线，因此估计了点云的法向。
// 在循环计算后，所有的云将被插入到training_clouds向量中。T
// raining_normals和training_classes将存储对应对象的法线和类索引。
//======================================================================
   for (unsigned int i_cloud = 0; i_cloud < number_of_training_clouds - 1; i_cloud++)
   {
     pcl::PointCloud<pcl::PointXYZ>::Ptr tr_cloud(new pcl::PointCloud<pcl::PointXYZ> ());
    if ( pcl::io::loadPCDFile <pcl::PointXYZ> (argv[i_cloud * 2 + 1], *tr_cloud) == -1 )
       return (-1);
 
     pcl::PointCloud<pcl::Normal>::Ptr tr_normals = (new pcl::PointCloud<pcl::Normal>)->makeShared ();
     normal_estimator.setInputCloud (tr_cloud);
     normal_estimator.compute (*tr_normals);
 
    unsigned int tr_class = static_cast<unsigned int> (strtol (argv[i_cloud * 2 + 2], 0, 10));

    training_clouds.push_back (tr_cloud);
    training_normals.push_back (tr_normals);
    training_classes.push_back (tr_class);
  }

// =====================================================================
// 以下创建了特征估计器的实例，使用FPFH方法进行的特征估计。
//	在将其传递给ISM算法之前，必须完全设置好它。这就是我们定义所有特征估计设置的地方。
// =====================================================================
   pcl::FPFHEstimation<pcl::PointXYZ, pcl::Normal, pcl::Histogram<153> >::Ptr fpfh
     (new pcl::FPFHEstimation<pcl::PointXYZ, pcl::Normal, pcl::Histogram<153> >);
  fpfh->setRadiusSearch (30.0);
 pcl::Feature< pcl::PointXYZ, pcl::Histogram<153> >::Ptr feature_estimator(fpfh);
 
// 在此实例中提供了训练数据和特征估计器。最后一行提供了前面提到的用于云简化的采样大小值。
  pcl::ism::ImplicitShapeModelEstimation<153, pcl::PointXYZ, pcl::Normal> ism;
  ism.setFeatureEstimator(feature_estimator);
  ism.setTrainingClouds (training_clouds);
  ism.setTrainingNormals (training_normals);
  ism.setTrainingClasses (training_classes);
  ism.setSamplingSize (2.0f);
// ======================================================================
//这一行代码只是创建了一个pcl::ism:: implicitshapemodelestimate的实例。它是一个有三个参数的模板类
// FeatureSize -特征(直方图)的大小，以计算
// PointT -点模板类
// NormalT -法线模板类
// ======================================================================
  pcl::ism::ImplicitShapeModelEstimation<153, pcl::PointXYZ, pcl::Normal>::ISMModelPtr model (new pcl::features::ISMModel);
 	 // 启动训练过程
  ism.trainISM (model);
	// 保存训练模型
   std::string file ("trained_ism_model.txt");
   model->saveModelToFile (file);
 	// 加载训练模型文件(这一步其实没有必要，只为了展示如何从其他地方加载已经训练好的模型)
   model->loadModelFromfile (file);
 

   unsigned int testing_class = static_cast<unsigned int> (strtol (argv[argc - 1], 0, 10));
   pcl::PointCloud<pcl::PointXYZ>::Ptr testing_cloud (new pcl::PointCloud<pcl::PointXYZ> ());
   if ( pcl::io::loadPCDFile <pcl::PointXYZ> (argv[argc - 2], *testing_cloud) == -1 )
     return (-1);
     
  // 分类过程需要点云及其法向，训练过程也需要点云和法向。因此，以下代码加载测试点云并计算法线。
   pcl::PointCloud<pcl::Normal>::Ptr testing_normals = (new pcl::PointCloud<pcl::Normal>)->makeShared ();
   normal_estimator.setInputCloud (testing_cloud);
   normal_estimator.compute (*testing_normals);
 
  // 以下启动分类过程。
  // 它告诉算法在给定的云testing_cloud中查找类型为testing_class的对象。
  // 注意，该算法将使用您通过的任何训练过的模型。分类完成后，将返回中心的投票名单。
  // pcl::ism::ISMVoteList是一个单独的类，它的目的是帮助分析投票。
   pcl::features::ISMVoteList<pcl::PointXYZ>::Ptr vote_list = ism.findObjects (
     model,
     testing_cloud,
     testing_normals,
     testing_class);
 
	// 以下代码负责在选票中找到最强的峰值。该搜索基于非最大抑制思想，这就是为什么非最大半径等于从训练模型中获得的目标半径。
   double radius = model->sigmas_[testing_class] * 10.0;
   double sigma = model->sigmas_[testing_class];
   std::vector<pcl::ISMPeak, Eigen::aligned_allocator<pcl::ISMPeak> > strongest_peaks;
   vote_list->findStrongestPeaks (strongest_peaks, testing_class, radius, sigma);
 
 // 以下是代码的其余部分，非常简单。它负责将云可视化并计算出最强的峰值，这些峰值代表testing_class类型对象的估计中心。
   pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = (new pcl::PointCloud<pcl::PointXYZRGB>)->makeShared ();
   colored_cloud->height = 0;
   colored_cloud->width = 1;
 
   pcl::PointXYZRGB point;
   point.r = 255;
   point.g = 255;
   point.b = 255;
 
   for (std::size_t i_point = 0; i_point < testing_cloud->size (); i_point++)
   {
     point.x = (*testing_cloud)[i_point].x;
     point.y = (*testing_cloud)[i_point].y;
     point.z = (*testing_cloud)[i_point].z;
     colored_cloud->points.push_back (point);
   }
   colored_cloud->height += testing_cloud->size ();

  point.r = 255;
  point.g = 0;
  point.b = 0;
 for (std::size_t i_vote = 0; i_vote < strongest_peaks.size (); i_vote++)
  {
    point.x = strongest_peaks[i_vote].x;
    point.y = strongest_peaks[i_vote].y;
    point.z = strongest_peaks[i_vote].z;
    colored_cloud->points.push_back (point);
  }
  colored_cloud->height += strongest_peaks.size ();

  pcl::visualization::CloudViewer viewer ("Result viewer");
  viewer.showCloud (colored_cloud);
  while (!viewer.wasStopped ())
  {
  }

  return (0);
}

---

【博主简介】
  斯坦福的兔子，男，天津大学机械工程工学硕士。毕业至今从事光学三维成像及点云处理相关工作。因工作中使用的三维处理库为公司内部库，不具有普遍适用性，遂自学开源PCL库及其相关数学知识以备使用。谨此将自学过程与君共享。
博主才疏学浅，尚不具有指导能力，如有问题还请各位在评论处留言供大家共同讨论。
若前辈们有工作机会介绍欢迎私信。