隧道点云处理（一）：投影边界提取 （Alpha Shapes）

1 算法溯源

Alpha Shapes算法 最早由 Edelsbrunner${}^{[1]}$(1983) 提出。沈蔚${}^{[2,3]}$(2008)首次将该算法应用于LIDAR数据处理（建筑物轮廓提取），并证明了该算法是一种简单高效、稳定的边界点提取算法。

如 图1 所示，用一个半径为 $\alpha$ 的圆在点集 P 外滚动，当 $\alpha$ 足够大时，这个圆就不会滚到点集内部，其滚动的痕迹就是点集 P 的边界线。故理论上存在：

• 当 $\alpha$ 值足够小，趋向于 $0$ 时，点集 P 中的每一个点都是边界点；
• 如果 $\alpha$ 值足够大，趋向于无穷时，Alpha Shapes( $\alpha \to \infty$) 是点集 P 的凸包。当 P 中点密度比较均匀，且 $\alpha$ 取适当的值时，Alpha Shapes 可同时提取点集 P 的内外边界。

图 1 Alpha Shapes 边界提取示意图

2 算法原理

2.1 判断条件

对于有限点集 P，由 n 个点组成，那么这 n 个点可以构成 $n\times (n-1)$ 条有向线段。我们的任务判断哪些线段为边界线段。

Alpha Shapes 模型的判断条件：

在点集 P 内，过任意两点 $p_1$，$p_2$ 绘制半径为 $\alpha$ 的圆，如果这个圆内没有其他点，则认为 $p_1$，$p_2$ 为边界点，其连线 $p_1{p}_2$ 为边界线段。

如 图 2 所示，已知点 $p_1(x_1,y_1)$，$p_2(x_2,y_2)$，求过这两点圆的圆心 $p_3(x_3,y_3)$，可由 后方距离交会法 求得。
$$\{\begin{array}{l}{x}_3=x_1+\frac{1}{2}(x_2-x_1)+H(y_2-y_1)\\ y_3=y_1+\frac{1}{2}(y_2-y_1)+H(x_1-x_2)\end{array}$$
其中，$$H=\sqrt{\frac{{\alpha }^2}{(x_1-x_2{)}^2-(y_1-y_2{)}^2}-\frac{1}{4}}$$

图 2 距离交会法求圆心

获得圆心后，通过判断其他点到圆心的距离与半径 $\alpha$ 的关系即可判断圆内是否存在其他点。

2.2 算法步骤

（1）从点集 P 中任意点 $p_i$ 开始，在半径为 $2\alpha$ 的邻域点集 R 内，取任意点 $p_r$，计算由点 $p_i$ 与 点 $p_r$ 确定的圆心 $p_0$；
（2）计算邻域内其他点到圆心 $p_0$ 的距离 $d_i$。若所有 $d_i$ 均大于 $\alpha$ ，说明圆内没有其他点，则 $p_i{p}_r$ 为边界线段；否则为非边界线段；
（3）对邻域点集 R 中的剩余点重复前两步，直到 R 内的点全部判断完成；
（4）对点集 P 中的剩余点重复上述步骤，直到 P 中的点全部判断完成。

图 3 边界点判断示意图

图 4 边界点提取结果示意图

3 代码实现

步骤： 读取隧道点云→水平投影→投影下采样→边界点提取→保存边界点

由于水平投影后，边界点密度要大于内部点密度，因此要进行下采样，使投影点均匀分布

代码：

#include <pcl/io/pcd_io.h>
#include <pcl/filters/voxel_grid.h>
#include <pcl/surface/concave_hull.h>
#include <pcl/console/time.h>

using namespace std;

typedef pcl::PointXYZ PointT;
typedef pcl::PointCloud<PointT> PointCloudT;

int main()
{
	//----------------------------- 加载点云 ----------------------------
	cout << "->正在加载点云..." << endl;
	PointCloudT::Ptr cloud(new PointCloudT);
	if (pcl::io::loadPCDFile("test.pcd", *cloud) < 0)
	{
		PCL_ERROR("->点云文件不存在！\a\n");
		return -1;
	}
	cout << "->共加载 " << cloud->points.size() << " 个数据点" << endl;
	//==================================================================


	//----------------------------- 点云投影 ----------------------------
	cout << "->正在水平面投影..." << endl;
	for (size_t i = 0; i < cloud->size(); ++i)
	{
		cloud->points[i].z = -1.0f;
	}
	//==================================================================

	//-------------------------- 投影边界下采样 --------------------------
	PointCloudT::Ptr cloud_sub(new PointCloudT);	//投影下采样点云
	pcl::VoxelGrid<PointT> vg;			//创建体素下采样对象
	vg.setInputCloud(cloud);			//设置下采样输入点云
	vg.setLeafSize(0.05f, 0.05f, 0.05f);//设置体素栅格边长
	vg.filter(*cloud_sub);				//执行体素下采样
	//==================================================================

	//-------------------- Alpha Shapes 提取投影边界 --------------------
	pcl::console::TicToc time;
	time.tic();
	cout << "->正在提取边界..." << endl;
	PointCloudT::Ptr cloud_boundary(new PointCloudT);
	pcl::ConcaveHull<PointT> ch;
	ch.setInputCloud(cloud_sub);
	ch.setAlpha(0.15);
	ch.reconstruct(*cloud_boundary);
	cout << "->边界点提取用时：" << time.toc() / 1000 << " s" << endl;
	//==================================================================


	//---------------------------- 保存边界点云 --------------------------
	cout << "->正在保存边界点云..." << endl;
	if (!cloud_boundary->empty())
	{
		pcl::io::savePCDFileBinary("Alpha_Shapes.pcd", *cloud_boundary);
		cout << "->共提取到 " << cloud_boundary->points.size() << " 个边界点" << endl;
	}
	else
	{
		PCL_ERROR("边界点云为空！\a\n");
		system("pause");
		return -1;
	}
	//==================================================================

	return 0;
}

输出结果：

->正在加载点云...
->共加载 744705 个数据点
->正在水平面投影...
->正在提取边界...
->边界点提取用时：9.63 s
->正在保存边界点云...
->共提取到 3782 个边界点

3 结果展示

隧道投影点云： 为便于展示，进行间隔为0.6m的空间距离下采样

边界点云：

$LeafSize$ 为体素下采样半径

投影点云 & 边界点云

4 半径 $\alpha$与边界提取结果的关系

$\alpha =k\cdot LeafSize$	边界点数	耗时(s)	结果展示
k=1	58096	10.975
k=1.5	4091	9.701
k=2	3962	9.647
k=3	3782	9.841
k=4	3654	9.893
k=6	3390	9.677
k=8	3141	9.665
k=10	2945	9.732

结论：

半径 $\alpha$ 应大于体素下采样半径， $k=2$ 时，提取效果较好。继续增大倍数，提取的边界点数越少，但变化不大，且提取效率基本相同。

5 参考文献

[1] on the shape of a set of points in the plane
[2] Building boundary extraction based on LIDAR point clouds data
[3] 沈蔚,李京,陈云浩,等.基于LIDAR数据的建筑轮廓线提取及规则化算法研究[J].遥感学报,2008(05):692-698.