《点云处理》点云边缘提取

Meditation-

已于 2024-01-04 16:26:01 修改

阅读量1k

点赞数 6

文章标签：算法 c++ 点云点云处理

于 2023-12-18 20:54:24 首次发布

本文链接：https://blog.csdn.net/dyk4ever/article/details/135071536

版权

前言

点云边缘或者说点云边界提取是点云处理过程中经常需要用到的方法。目前，常用的点云边缘提取的方法有凹包检测、基于法向和曲率的边缘提取、经纬线扫描法等。最终，还可能会涉及到使用提取的边界进行点云分割。

环境：
Windows11 + PCL 1.11.1

1.凹包检测

凹包检测是识别点云中凹多边形的过程，PCL中有集成的凹包检测算法。但是使用该方法需要注意的是，如果是多线程运行的软件或者算法，并且在每一个线程中都需要使用凹包检测。那么可能会存在问题。因为凹包检测需要单线程运行，这可能是与其内部实现方式有关系，通过加锁就可以解决掉凹包检测引起的程序崩溃问题。

代码：

/// <summary>
/// 使用凹包检测算法进行点云轮廓提取
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloudout">输出点云</param>
/// <param name="polygons">对应于轮廓的索引点</param>
/// <param name="alpha">系数</param>
/// <returns>return false表示轮廓提取失败</returns>
bool ConcaveHull(const pcl::shared_ptr<pcl::PointCloud<pcl::PointXYZ>>& cloud, pcl::shared_ptr<pcl::PointCloud<pcl::PointXYZ>>& cloudout, std::vector<pcl::Vertices>& polygons, const double& alpha)
{
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 1000) return false;
	if (cloudout == nullptr) cloudout.reset(new pcl::PointCloud<pcl::PointXYZ>);
	pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());
	pcl::ConcaveHull<pcl::PointXYZ> chull;
	chull.setInputCloud(cloud);													// 输入点云为投影后的点云
	chull.setAlpha(alpha);														// 设置alpha值为0.5
	chull.setSearchMethod(tree);
	chull.reconstruct(*cloudout, polygons);

	return true;
}

可以使用OpenMP加锁，比较常用的是critical语句。

#pragma omp critical
{
	ConcaveHull(cloud, cloudout, polygons, alpha);
}

一般alpha都设定为0.5，如果点云比较稀疏，则需要适当增大这个数值。

2.凸包检测

凸包检测与其说是提取轮廓，其实不如说是提取点云凸多边形的顶点。还可以得到凸多边形的面积大小，一般来说都会设定为2维，如果是三维，那么还可以计算得到凸多变体的体积。总的来说，使用凸包检测是很难得到一个稠密整齐的轮廓点云的，而是得到顶点。

/// <summary>
/// 使用凸包检测提取点云轮廓顶点
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloudout">输出点云</param>
/// <param name="polygons">轮廓点云顶点</param>
/// <returns>return false表示轮廓顶点提取失败</returns>
bool ConvexHull(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr& cloudout, std::vector<pcl::Vertices>& polygons)
{

	pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());
	pcl::ConvexHull<pcl::PointXYZ> hull;
	hull.setInputCloud(cloud);
	hull.setDimension(2);              // 设置凸包维度
	hull.setSearchMethod(tree);
	hull.reconstruct(*cloudout, polygons);
	//double convex_volume = hull.getTotalVolume();

	return true;
}

凹包检测点云分割
与凹包检测或者凸包检测可以配合使用的是crophull，使用该方法可以利用凹包检测或者凸包得到的轮廓点云cloudout和polygons对其他点云进行分割。使用这种方法去分割点云只适用于点云数量不大，点云数量越大，速度越慢。

代码如下：

/// <summary>
/// 使用凹包或者凸包检测结果进行点云分割
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="hullcloud">轮廓点云</param>
/// <param name="polygons">多边形顶点</param>
/// <param name="cloudout">输出点云</param>
/// <param name="is_in">true是提取轮廓内的点，false提取轮廓外的点</param>
/// <returns>return true表示分割成功</returns>
bool CropHull(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, const pcl::PointCloud<pcl::PointXYZ>::Ptr& hullcloud, 
	const std::vector<pcl::Vertices>& polygons, pcl::PointCloud<pcl::PointXYZ>::Ptr& cloudout, bool is_in)
{
	if (cloud == nullptr || hullcloud == nullptr || cloudout == nullptr) return false;
	pcl::CropHull<pcl::PointXYZ> SURF_HULL;
	SURF_HULL.setDim(2);
	SURF_HULL.setInputCloud(cloud);           // 输入点云
	SURF_HULL.setHullIndices(polygons);       // 凸包或者凹包的索引值
	SURF_HULL.setHullCloud(hullcloud);        // 凸包或者凹包的点云
	SURF_HULL.setCropOutside(is_in);          // false剔除 true 提取
	SURF_HULL.filter(*cloudout);
	
	return true;
}

3.经纬线扫描法

从这篇帖子看到了通过经纬线扫描法这种轮廓提取方法，这个方法特别适用于使用线扫相机得到的点云。而且只是遍历了点云两遍，速度还比较快。但是如果得到的点云使用体素滤波进行稀疏的话，通常稀疏点云是错乱的，即排列不整齐，在x和y方向上。这种稀释后的点云使用经纬线扫描法得到的轮廓通常也不是一个干净整齐的轮廓。此外，还存在一个问题，如果点云空间外面漂浮了一些噪点，那么这些点可能会被扫描成轮廓。所以就需要扫描前对原始点云进行去噪，也可以扫描之后对轮廓点云进行去噪。但是如果噪点较多，很可能导致得到的轮廓去噪后产生缺口，轮廓不闭合。所以到底是轮廓提取前去噪还是提取之后去噪需要根据点云质量而定。

代码：

/// <summary>
/// 使用经纬线扫描法进行点云轮廓提取
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloud_boundary">输出点云边界</param>
/// <param name="resolution">单方向上扫描滑动次数</param>
/// <returns>return true表示轮廓提取成功</returns>
bool BoundaryExtraction(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_boundary, int resolution)
{
	if (cloud == nullptr) return false;
	if (cloud->points.size() < 100) return false;
	
	pcl::PointXYZ px_min = *std::min_element(cloud->begin(), cloud->end(), [](pcl::PointXYZ pt1, pcl::PointXYZ pt2) {return pt1.x < pt2.x; });
	pcl::PointXYZ px_max = *std::max_element(cloud->begin(), cloud->end(), [](pcl::PointXYZ pt1, pcl::PointXYZ pt2) {return pt1.x < pt2.x; });

	float delta_x = (px_max.x - px_min.x) / resolution;
	float min_y = INT_MAX, max_y = -INT_MAX;
	std::vector<int> indexs_x(2 * resolution + 2);
	std::vector<std::pair<float, float>> minmax_x(resolution + 1, { INT_MAX,-INT_MAX });

#pragma omp parallel for 
	for (int i = 0; i < cloud->size(); ++i)
	{
		int id = (cloud->points[i].x - px_min.x) / delta_x;
		if (cloud->points[i].y < minmax_x[id].first)
		{
			minmax_x[id].first = cloud->points[i].y;
			indexs_x[id] = i;
		}
		else if (cloud->points[i].y > minmax_x[id].second)
		{
			minmax_x[id].second = cloud->points[i].y;
			indexs_x[id + resolution + 1] = i;
		}
	}

	pcl::PointXYZ py_min = *std::min_element(cloud->begin(), cloud->end(), [](pcl::PointXYZ pt1, pcl::PointXYZ pt2) {return pt1.y < pt2.y; });
	pcl::PointXYZ py_max = *std::max_element(cloud->begin(), cloud->end(), [](pcl::PointXYZ pt1, pcl::PointXYZ pt2) {return pt1.y < pt2.y; });

	float delta_y = (py_max.y - py_min.y) / resolution;
	float min_x = INT_MAX, max_x = -INT_MAX;
	std::vector<int> indexs_y(2 * resolution + 2);
	std::vector<std::pair<float, float>> minmax_y(resolution + 1, { INT_MAX,-INT_MAX });

#pragma omp parallel for
	for (int i = 0; i < cloud->size(); ++i)
	{
		int id = (cloud->points[i].y - py_min.y) / delta_y;
		if (cloud->points[i].x < minmax_y[id].first)
		{
			minmax_y[id].first = cloud->points[i].x;
			indexs_y[id] = i;
		}
		else if (cloud->points[i].x > minmax_y[id].second)
		{
			minmax_y[id].second = cloud->points[i].x;
			indexs_y[id + resolution + 1] = i;
		}
	}

	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_xboundary(new pcl::PointCloud<pcl::PointXYZ>);
	pcl::copyPointCloud(*cloud, indexs_x, *cloud_xboundary);
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_yboundary(new pcl::PointCloud<pcl::PointXYZ>);
	pcl::copyPointCloud(*cloud, indexs_y, *cloud_yboundary);
	*cloud_boundary = *cloud_xboundary + *cloud_yboundary;

	return true;
}

4.基于法向轮廓提取

PCL中封装有基于法向的轮廓提取算法，直接称之为轮廓提取。

代码如下：

/// <summary>
/// PCL点云边界提取
/// </summary>
/// <param name="cloud">输入点云</param>
/// <param name="cloudout">输出边界点云</param>
/// <param name="angle">边界角度阈值</param>
/// <param name="radius">搜索半径</param>
/// <returns>return true表示边界提取成功</returns>
bool BoundaryEstimate(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr& cloudout, const float& angle, 
	const double& radius)
{
	if (cloudout == nullptr) cloudout.reset(new pcl::PointCloud<pcl::PointXYZ>);
	// 多线程法线估计
	pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());
	pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);
	pcl::NormalEstimationOMP<pcl::PointXYZ, pcl::Normal> ne;
	ne.setInputCloud(cloud);
	ne.setSearchMethod(tree);
	ne.setRadiusSearch(radius); // 法线估计的搜索半径
	ne.setNumberOfThreads(12);  // 设置线程数
	ne.compute(*normals);
	// 曲率边缘提取
	pcl::PointCloud<pcl::Boundary> boundaries;
	pcl::BoundaryEstimation<pcl::PointXYZ, pcl::Normal, pcl::Boundary> est;
	est.setInputCloud(cloud);
	est.setInputNormals(normals);
	est.setAngleThreshold(angle * M_PI / 180 );
	//est.setRadiusSearch(radius); // 边缘提取的搜索半径
	est.setKSearch(50);
	est.setSearchMethod(tree);
	est.compute(boundaries);
	// 提取边缘点
	
	for (size_t i = 0; i < cloud->size(); ++i)
	{
		if (boundaries[i].boundary_point > 0) 
		{
			cloudout->points.emplace_back(cloud->at(i));
		}
	}

	return true;
}