点云边界提取方法总结

经纬线扫描法

经纬线扫描法的基本思想是：将经过坐标变换后的二维数据点集按照 x值的大小排序后存储在一个链表中，在二维平面建立点集的最小包围盒并分别计算出 x 和 y 的最大、最小值；令经线从 x 的最小值开始，取步长为dx，在 x 的取值范围内分别计算出每根经线的最大和最小 y 值，并将它们的索引值放在一个新建的链表中，扫描完整个 x 区间；同样的方法扫描纬线，最后将重复的索引值删掉。
下面代码在点云主法线方向和z轴夹角较小时表现较好，其他情况需要修改。

#include <iostream>
#include <algorithm>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>void BoundaryExtraction(pcl::PointCloud<pcl::PointXYZ>::Ptr cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_boundary, int resolution)
{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);std::vector<std::pair<float, float>> minmax_x(resolution, { INT_MAX,-INT_MAX });for (size_t 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] = 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);std::vector<std::pair<float, float>> minmax_y(resolution, { INT_MAX,-INT_MAX });for (size_t 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] = 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;
}int main(int argc, char* argv[])
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>),;pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_boundary(new pcl::PointCloud<pcl::PointXYZ>);pcl::io::loadPCDFile("test.pcd", *cloud);BoundaryExtraction(cloud, cloud_boundary, 200);pcl::io::savePCDFile("boundary.pcd", *cloud_boundary);system("pause");return EXIT_SUCCESS;
}

提取结果：

网格划分法

网格划分法分为三个步骤：（1）网格划分（2）寻找边界网格（3）提取边界线。首先建立数据点集的最小包围盒，并用给定间隔的矩形网格将其分割；然后找出边界网格，将这些边界网格按顺序连接起来形成一条由边界网格组成的“粗边界”；再对每个边界网格按照某种规则判断其内的点是否为边界点，连接初始边界，并对此边界线进一步处理使其光滑。

#include <iostream>
#include <vector>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/common/transforms.h>void BoundaryExtraction( pcl::PointCloud<pcl::PointXYZ>::Ptr cloud,  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_boundary, int resolution)
{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; });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 x_min = px_min.x, x_max = px_max.x, y_min = py_min.y, y_max = py_max.y;float L = sqrt((x_max - x_min)*(y_max - y_min) / resolution);int x_num = (x_max - x_min) / L + 1;int y_num = (y_max - y_min) / L + 1;std::vector<std::vector< pcl::PointCloud<pcl::PointXYZ>::Ptr>> uv(x_num + 1, std::vector< pcl::PointCloud<pcl::PointXYZ>::Ptr>(y_num + 1));for (int i = 0; i <= x_num; ++i){for (int j = 0; j <= y_num; ++j){pcl::PointCloud<pcl::PointXYZ>::Ptr ptcloud(new  pcl::PointCloud<pcl::PointXYZ>);uv[i][j] = ptcloud;}}for (int i = 0; i < cloud->size(); ++i){int x_id = (cloud->points[i].x - x_min) / L;int y_id = (cloud->points[i].y - y_min) / L;uv[x_id][y_id]->push_back(cloud->points[i]);}std::vector<std::vector<bool>> boundary_index(x_num + 1, std::vector<bool>(y_num + 1, false));for (int i = 0; i <= x_num; ++i){if (uv[i][0]->size())boundary_index[i][0] = true;if (uv[i][y_num]->size())boundary_index[i][y_num] = true;}     for (int i = 0; i <= y_num; ++i){if(uv[0][i]->size())boundary_index[0][i] = true;if (uv[x_num][i]->size())boundary_index[x_num][i] = true;}for (int i = 1; i < x_num - 1; ++i){for (int j = 1; j < y_num - 1; ++j){if (uv[i][j]->size()){if (!uv[i - 1][j - 1]->size() || !uv[i][j - 1]->size() || !uv[i + 1][j - 1]->size() || !uv[i][j - 1]->size() ||!uv[i + 1][j - 1]->size() || !uv[i + 1][j]->size() || !uv[i + 1][j + 1]->size() || !uv[i][j + 1]->size())boundary_index[i][j] = true;}}}for (int i = 0; i <= x_num; ++i){for (int j = 0; j <= y_num; ++j){if (boundary_index[i][j]){pcl::PointCloud<pcl::PointXYZ>::Ptr ptcloud(new  pcl::PointCloud<pcl::PointXYZ>);ptcloud = uv[i][j];Eigen::Vector4f cloud_centroid;pcl::compute3DCentroid(*ptcloud, cloud_centroid);cloud_boundary->push_back(pcl::PointXYZ(cloud_centroid(0), cloud_centroid(1), cloud_centroid(2)));}}}
}int main(int argc, char* argv[])
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new  pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_boundary(new  pcl::PointCloud<pcl::PointXYZ>);pcl::io::loadPCDFile("test.pcd", *cloud);BoundaryExtraction(cloud, cloud_boundary, 200*200);pcl::io::savePCDFile("boundary.pcd", *cloud_boundary);system("pause");return 0;
}

提取结果：

法线估计法

可参考PCL：点云数据基于法线的边界提取(从最初的法线估计理论推导到最终的边界提取)

#include <iostream>
#include <vector>
#include <ctime>
#include <boost/thread/thread.hpp>
#include <pcl/io/pcd_io.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/console/parse.h>
#include <pcl/features/eigen.h>
#include <pcl/features/feature.h>
#include <pcl/features/normal_3d.h>
#include <pcl/impl/point_types.hpp>
#include <pcl/features/boundary.h>
#include <pcl/visualization/cloud_viewer.h>int main(int argc, char* argv[])
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);pcl::io::loadPCDFile("test.pcd", *cloud);pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);pcl::PointCloud<pcl::Boundary> boundaries;pcl::BoundaryEstimation<pcl::PointXYZ, pcl::Normal, pcl::Boundary> est;pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normEst; normEst.setInputCloud(cloud);normEst.setSearchMethod(tree);normEst.setRadiusSearch(3);  //法向估计的半径  3//normEst.setKSearch(9);  //法向估计的点数normEst.compute(*normals);est.setInputCloud(cloud);est.setInputNormals(normals);est.setSearchMethod(tree);est.setKSearch(500);  //一般这里的数值越高，最终边界识别的精度越好est.compute(boundaries);pcl::PointCloud<pcl::PointXYZ>::Ptr boundPoints(new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ> noBoundPoints;int countBoundaries = 0;for (int i = 0; i<cloud->size(); i++) {uint8_t x = (boundaries.points[i].boundary_point);int a = static_cast<int>(x); //该函数的功能是强制类型转换if (a == 1){(*boundPoints).push_back(cloud->points[i]);countBoundaries++;}elsenoBoundPoints.push_back(cloud->points[i]);}pcl::io::savePCDFile("boundary.pcd",*boundPoints);system("pause");return EXIT_SUCCESS;
}

提取结果：

alpha shapes算法

原理参考：alpha shapes提取平面点云边界点

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/surface/concave_hull.h>int main(int argc, char* argv[])
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);pcl::io::loadPCDFile("test.pcd", *cloud);pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_hull(new pcl::PointCloud<pcl::PointXYZ>); pcl::ConcaveHull<pcl::PointXYZ> chull;chull.setInputCloud(cloud); //输入点云为投影后的点云chull.setAlpha(10); //设置aLpha值chull.reconstruct(*cloud_hull);pcl::io::savePCDFile("boundary.pcd", *cloud_hull);system("pause");return EXIT_SUCCESS;
}

提取结果：

参考文献：[1]杨雪娇. 点云的边界提取及角点检测算法研究[D].哈尔滨工程大学,2010.