PCL教程-点云分割之区域生长分割算法

朱雀 2022-09-10 15:09 370阅读 0赞

> 原文链接：[Region growing segmentation][]

本教程使用到的点云数据：[source files][]

在本篇教程中，我们将学习如何使用由**pcl::RegionGrowing类**实现的区域生长算法。该算法的目的是合并在平滑约束条件下足够接近的点。因此，该算法的输出数据结构是由聚类组成的数组，其中每个聚类都是被认为是同一光滑表面的一部分的点的集合。该算法的工作原理（光滑度的计算）是基于两点法线之间的角度比较。

**目录**

基本原理

算法伪代码

输入

初始化

算法实现

程序代码

代码分析

实验结果

其他数据结果

--------------------

# 基本原理 #

首先，它根据点的曲率值对点进行排序。需要这样做是因为区域从曲率最小的点开始增长。这样做的原因是曲率最小的点位于平坦区域(从最平坦区域生长可以减少段的总数)。

我们有了分类过的云。直到云中没有未标记点时，算法选取曲率值最小的点，开始区域的增长。这个过程如下所示：

*  选中的点被添加到名为种子的集合中。
 *  对于每一个种子点，找到它的邻近点：
    
     *  算出每个相邻点的法线和当前种子点的法线之间的角度，如果角度小于阈值，则将当前点添加到当前区域。
     *  然后计算每个邻居点的曲率值，如果曲率小于阈值，那么这个点被添加到种子中。
     *  将当前的种子从种子列表中移除。

如果种子列表变成空的，这意味着该区域生长已完成，继续重复上述过程。

> **区域生长算法： 将具有相似性的点云集合起来构成区域。**  
> 首先对每个需要分割的区域找出一个种子点作为生长的起点，然后将种子点周围邻域中与种子有相同或 相似性质的点合并到种子像素所在的区域中。而新的点继续作为种子向四周生长，直到再没有满足条件 的像素可以包括进来，一个区域就生长而成了。  
> 算法流程：
> 
> 1.  计算法线normal 和曲率curvatures，依据曲率升序 排序；
> 2.  选择曲率最低的为初始种子点，种子周围的临近点和 种子点云相比较；
> 3.  法线的方向是否足够相近（法线夹角足够 r p y）， 法线夹角阈值；
> 4.  曲率是否足够小（表面处在同一个弯曲程度），区域 差值阈值；
> 5.  如果满足2，3则该点可用做种子点; 6. 如果只满足2，则归类而不做种子;

# 算法伪代码 #

## 输入 ##

*  *Point cloud* = ![\\\{P\\\}][P]
 *  *Point normals* = ![\\\{N\\\}][N]
 *  *Points curvatures* = ![\\\{c\\\}][c]
 *  *Neighbour finding function* ![\\Omega(.)][Omega_.]
 *  *Curvature threshold* ![c\_\{th\}][c_th]
 *  *Angle threshold* ![\\theta\_\{th\}][theta_th]

## 初始化 ##

*区域列表置为空： Region list* ![\{R\}\\leftarrow\{\\O\}][R_leftarrow_O]

可用的点云列表：*Available points list* ![\\\{A\\\}\\leftarrow\\\{1,...,|P|\\\}][A_leftarrow_1_..._P]

## 算法实现 ##

![20210909143940903.png][]

# 程序代码 #

#include <iostream>
    #include <vector>
    #include <pcl/point_types.h>
    #include <pcl/io/pcd_io.h>
    #include <pcl/search/search.h>
    #include <pcl/search/kdtree.h>
    #include <pcl/features/normal_3d.h>
    #include <pcl/visualization/cloud_viewer.h>
    #include <pcl/filters/passthrough.h>
    #include <pcl/segmentation/region_growing.h>
    #include <pcl/console/print.h>
    #include <pcl/console/parse.h>
    #include <pcl/console/time.h>
    #include <windows.h>
    #include <stdio.h>
    #include <psapi.h>
    void PrintMemoryInfo( )
    {
    	HANDLE hProcess;
    	PROCESS_MEMORY_COUNTERS pmc;
    
    	hProcess=GetCurrentProcess();
    	printf( "\nProcess ID: %u\n", hProcess );
    
    	// Print information about the memory usage of the process.
    	//输出进程使用的内存信息
       
    	if (NULL == hProcess)
    		return;
    
    	if ( GetProcessMemoryInfo( hProcess, &pmc, sizeof(pmc)) )
    	{
    		printf( "\tPageFaultCount: 0x%08X\n", pmc.PageFaultCount );
    		printf( "\tPeakWorkingSetSize: 0x%08X\n", 
    				  pmc.PeakWorkingSetSize );
    		printf( "\tWorkingSetSize: 0x%08X\n", pmc.WorkingSetSize );
    		printf( "\tQuotaPeakPagedPoolUsage: 0x%08X\n", 
    				  pmc.QuotaPeakPagedPoolUsage );
    		printf( "\tQuotaPagedPoolUsage: 0x%08X\n", 
    				  pmc.QuotaPagedPoolUsage );
    		printf( "\tQuotaPeakNonPagedPoolUsage: 0x%08X\n", 
    				  pmc.QuotaPeakNonPagedPoolUsage );
    		printf( "\tQuotaNonPagedPoolUsage: 0x%08X\n", 
    				  pmc.QuotaNonPagedPoolUsage );
    		printf( "\tPagefileUsage: 0x%08X\n", pmc.PagefileUsage ); 
    		printf( "\tPeakPagefileUsage: 0x%08X\n", 
    				  pmc.PeakPagefileUsage );
    	}
    
    	CloseHandle( hProcess );
    }
    
    using namespace pcl::console;
    int
    main (int argc, char** argv)
    {
    
    	if(argc<2)
    	{
    		std::cout<<".exe xx.pcd -kn 50 -bc 0 -fc 10.0 -nc 0 -st 30 -ct 0.05"<<endl;
    
    		return 0;
    	}//如果输入参数小于1个，输出提示
    	time_t start,end,diff[5],option;
    	start = time(0); 
    	int KN_normal=50; //设置默认输入参数
    	bool Bool_Cuting=false;//设置默认输入参数
    	float far_cuting=10,near_cuting=0,SmoothnessThreshold=30.0,CurvatureThreshold=0.05;//设置默认输入参数
    	parse_argument (argc, argv, "-kn", KN_normal);
    	parse_argument (argc, argv, "-bc", Bool_Cuting);
    	parse_argument (argc, argv, "-fc", far_cuting);
    	parse_argument (argc, argv, "-nc", near_cuting);
    	parse_argument (argc, argv, "-st", SmoothnessThreshold);
    	parse_argument (argc, argv, "-ct", CurvatureThreshold);//设置输入参数方式
    	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
    	if ( pcl::io::loadPCDFile <pcl::PointXYZ> (argv[1], *cloud) == -1)
    	{
    		std::cout << "Cloud reading failed." << std::endl;
    		return (-1);
    	}// 加载输入点云数据
    	end = time(0); 
    	diff[0] = difftime (end, start); 
    	PCL_INFO ("\Loading pcd file takes(seconds): %d\n", diff[0]); 
    	//Noraml estimation step(1 parameter)
    	pcl::search::Search<pcl::PointXYZ>::Ptr tree = boost::shared_ptr<pcl::search::Search<pcl::PointXYZ> > (new pcl::search::KdTree<pcl::PointXYZ>);//创建一个指向kd树搜索对象的共享指针
    	pcl::PointCloud <pcl::Normal>::Ptr normals (new pcl::PointCloud <pcl::Normal>);
    	pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;//创建法线估计对象
    	normal_estimator.setSearchMethod (tree);//设置搜索方法
    	normal_estimator.setInputCloud (cloud);//设置法线估计对象输入点集
    	normal_estimator.setKSearch (KN_normal);// 设置用于法向量估计的k近邻数目
    	normal_estimator.compute (*normals);//计算并输出法向量
    	end = time(0); 
    	diff[1] = difftime (end, start)-diff[0]; 
    	PCL_INFO ("\Estimating normal takes(seconds): %d\n", diff[1]); 
    	//optional step: cutting the part are far from the orignal point in Z direction.2 parameters
    	pcl::IndicesPtr indices (new std::vector <int>);//创建一组索引
    	if(Bool_Cuting)//判断是否需要直通滤波
    	{
    
    		pcl::PassThrough<pcl::PointXYZ> pass;//设置直通滤波器对象
    		pass.setInputCloud (cloud);//设置输入点云
    		pass.setFilterFieldName ("z");//设置指定过滤的维度
    		pass.setFilterLimits (near_cuting, far_cuting);//设置指定纬度过滤的范围
    		pass.filter (*indices);//执行滤波，保存滤波结果在上述索引中
    	}
    
    
    	// 区域生长算法的5个参数
    	pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;//创建区域生长分割对象
    	reg.setMinClusterSize (50);//设置一个聚类需要的最小点数
    	reg.setMaxClusterSize (1000000);//设置一个聚类需要的最大点数
    	reg.setSearchMethod (tree);//设置搜索方法
    	reg.setNumberOfNeighbours (30);//设置搜索的临近点数目
    	reg.setInputCloud (cloud);//设置输入点云
    	if(Bool_Cuting)reg.setIndices (indices);//通过输入参数设置，确定是否输入点云索引
    	reg.setInputNormals (normals);//设置输入点云的法向量
    	reg.setSmoothnessThreshold (SmoothnessThreshold / 180.0 * M_PI);//设置平滑阈值
    	reg.setCurvatureThreshold (CurvatureThreshold);//设置曲率阈值
    
    	std::vector <pcl::PointIndices> clusters;
    	reg.extract (clusters);//获取聚类的结果，分割结果保存在点云索引的向量中。
    	end = time(0); 
    	diff[2] = difftime (end, start)-diff[0]-diff[1]; 
    	PCL_INFO ("\Region growing takes(seconds): %d\n", diff[2]); 
    
    	std::cout << "Number of clusters is equal to " << clusters.size () << std::endl;//输出聚类的数量
    	std::cout << "First cluster has " << clusters[0].indices.size () << " points." << endl;//输出第一个聚类的数量
    	std::cout << "These are the indices of the points of the initial" <<
    		std::endl << "cloud that belong to the first cluster:" << std::endl;
    	/* int counter = 0;
    	while (counter < clusters[0].indices.size ())
    	{
    	std::cout << clusters[0].indices[counter] << ", ";
    	counter++;
    	if (counter % 10 == 0)
    	std::cout << std::endl;
    	}
    	std::cout << std::endl;
    	*/
    	PrintMemoryInfo();
    	pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
    	pcl::visualization::CloudViewer viewer ("区域增长分割方法");
    	viewer.showCloud(colored_cloud);
    	while (!viewer.wasStopped ())
    	{
    	}//进行可视化
    
    	return (0);
    }

# 代码分析 #

然后设置最小和最大集群大小。这意味着在分割完成后，所有点小于最小值(或大于最大值)的聚类将被丢弃。最小值和最大值的默认值分别为1和“尽可能多”。

算法在内部结构中需要K最近邻搜索，所以这里是提供搜索方法并设置邻居数量的地方。之后，它接收到必须分割的点云、点下标和法线。

pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
      reg.setMinClusterSize (50);
      reg.setMaxClusterSize (1000000);
      reg.setSearchMethod (tree);
      reg.setNumberOfNeighbours (30);
      reg.setInputCloud (cloud);
      //reg.setIndices (indices);
      reg.setInputNormals (normals);

这两行是算法初始化中最重要的部分，因为它们负责上述的平滑约束。第一种方法以弧度为单位设置角度，作为法向偏差的允许范围。如果点之间的法线偏差小于平滑阈值，那么建议它们在同一簇中(新的点-被测试的点-将被添加到簇中)。第二个是曲率阈值。如果两点有一个小的法向偏差，那么它们之间的曲率差被测试。

reg.setSmoothnessThreshold (3.0 / 180.0 * M_PI);
      reg.setCurvatureThreshold (1.0);

RegionGrowing类提供了一个返回彩色云的方法，其中每个集群都有自己的颜色。因此，在这部分代码中，实例化pcl::visualization::CloudViewer以查看分割的结果——相同颜色的云

pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
      pcl::visualization::CloudViewer viewer ("Cluster viewer");
      viewer.showCloud(colored_cloud);
      while (!viewer.wasStopped ())
      {
      }
    
      return (0);
    }

# 实验结果 #

原始点云：

![20210909150126810.png][]

使用区域生长算法分割后的点云（每种颜色代表一个聚类）：

![20210909145758437.png][]

在最后一张图片中，你可以看到彩色云有许多**红点**。这意味着这些点属于被拒绝的聚类，因为它们有太多/太少的点。

使用命令行进行输入：

> D:\\PCLProject\\pcl-project\\pcl\_segmentation\\4\_region\_growing\_segmentation\\cmake\_bin\\Release>region\_growing\_segmentation.exe pig1.pcd -kn 50 -bc 0 -fc 10.0 -nc 0 -st 30 -ct 0.05  
> Loading pcd file takes(seconds): 0  
> Estimating normal takes(seconds): 1  
> Region growing takes(seconds): 0  
> Number of clusters is equal to 141  
> First cluster has 136600 points.  
> These are the indices of the points of the initial  
> cloud that belong to the first cluster:
> 
> Process ID: 4294967295  
> PageFaultCount: 0x00004A9E  
> PeakWorkingSetSize: 0x03A69000  
> WorkingSetSize: 0x03A68000  
> QuotaPeakPagedPoolUsage: 0x0002DD40  
> QuotaPagedPoolUsage: 0x0002D5B8  
> QuotaPeakNonPagedPoolUsage: 0x00003508  
> QuotaNonPagedPoolUsage: 0x00003480  
> PagefileUsage: 0x03427000  
> PeakPagefileUsage: 0x03427000

## 其他数据结果 ##

原始点云：

![20210909150603937.png][]

处理之后的点云：

![20210909150544227.png][]

[Region growing segmentation]: https://pcl.readthedocs.io/projects/tutorials/en/latest/region_growing_segmentation.html#region-growing-segmentation
[source files]: https://github.com/luolaihua/pcl-project/tree/master/pcl_segmentation/4_region_growing_segmentation/source
[P]: https://latex.codecogs.com/gif.latex?%5C%7BP%5C%7D
[N]: https://latex.codecogs.com/gif.latex?%5C%7BN%5C%7D
[c]: https://latex.codecogs.com/gif.latex?%5C%7Bc%5C%7D
[Omega_.]: https://latex.codecogs.com/gif.latex?%5COmega%28.%29
[c_th]: https://latex.codecogs.com/gif.latex?c_%7Bth%7D
[theta_th]: https://latex.codecogs.com/gif.latex?%5Ctheta_%7Bth%7D
[R_leftarrow_O]: https://latex.codecogs.com/gif.latex?%7BR%7D%5Cleftarrow%7B%5CO%7D
[A_leftarrow_1_..._P]: https://latex.codecogs.com/gif.latex?%5C%7BA%5C%7D%5Cleftarrow%5C%7B1%2C...%2C%7CP%7C%5C%7D
[20210909143940903.png]: /images/20220829/81889db30e6f411198fbf8e6c80b6599.png
[20210909150126810.png]: /images/20220829/13290f8c27494121a4cdc1ae1cd0a645.png
[20210909145758437.png]: /images/20220829/307f6ad9e4e84763b2f097f263c170ac.png
[20210909150603937.png]: /images/20220829/fb7f3929882f45ac8ae9d17913a99837.png
[20210909150544227.png]: /images/20220829/01493ff1e2a54fbfaaf0d157f62ca5ba.png