作者 | 老司机带我飞 编辑 | 汽车人
原文链接：https://www.zhihu.com/question/465073203/answer/2925688433
点击下方卡片，关注“自动驾驶之心”公众号
ADAS巨卷干货，即可获取
点击进入→自动驾驶之心【3D目标检测】技术交流群
后台回复【3D检测综述】获取最新基于点云/BEV/图像的3D检测综述！
最近因为要做货车自动装卸货，使用3D视觉获取车厢点云，需要滤除地面点云。目前滤波地面点有多种方式，然后实际项目使用了CSF和PCL Progressive Morphological Filter方式，对比了下，CSF效果好很多。以下关于两种方法的编译使用参考。
一、使用CSF对kitti的点云数据过滤出地面点云，结合PCL使用，C++实现
环境：win11+vs2022
使用CSF进行地面点云滤波，使用了PCL库读取显示点云，CSF算法使用github开源代码，自己编译生成CSF.lib。github地址：
https://github.com/jianboqi/CSF
点云测试数据下载链接：链接：https://pan.baidu.com/s/1mdZ4RBObvVZkyq9GOTv_hg 提取码：6hdx
#include <string>
#include <pcl/point_cloud.h>
#include <pcl/point_types.h>//PCL对各种格式的点的支持头文件
#include <pcl/io/pcd_io.h>//PCL的PCD格式文件的输入输出头文件
#include <pcl/visualization/cloud_viewer.h>//点云查看窗口头文件
#include <pcl/filters/filter.h> //滤波相关头文件
#include <pcl/filters/passthrough.h>
#include <pcl/filters/statistical_outlier_removal.h>
#include "CSF.h"
using namespace std;
void clothSimulationFilter(const vector< csf::Point >& pc,vector<int> &groundIndexes,vector<int> & offGroundIndexes)
{
//step 1 read point cloud
CSF csf;
csf.setPointCloud(pc);// or csf.readPointsFromFile(pointClouds_filepath);
//pc can be vector< csf::Point > or PointCloud defined in point_cloud.h
//step 2 parameter settings
//Among these paramters:
//time_step interations class_threshold can remain as defualt in most cases.
csf.params.bSloopSmooth = false;
csf.params.cloth_resolution = 0.5;
csf.params.rigidness = 3;
csf.params.time_step = 0.65;
csf.params.class_threshold = 0.5;
csf.params.interations = 500;
//step 3 do filtering
//result stores the index of ground points or non-ground points in the original point cloud
csf.do_filtering(groundIndexes, offGroundIndexes);
//csf.do_filtering(groundIndexes, offGroundIndexes,true);
//if the third parameter is set as true, then a file named "cloth_nodes.txt" will be created,
//it respresents the cloth nodes.By default, it is set as false
}
void addPointCloud(const vector<int>& index_vec, const pcl::PointCloud<pcl::PointXYZI>::Ptr cloud, pcl::PointCloud<pcl::PointXYZI>::Ptr cloud_filtered)
{
auto& points = cloud_filtered->points;
const auto& pointclouds = cloud->points;
for_each(index_vec.begin(), index_vec.end(), [&](const auto& index) {
pcl::PointXYZI pc;
pc.x = pointclouds[index].x;
pc.y = pointclouds[index].y;
pc.z = pointclouds[index].z;
pc.intensity = pointclouds[index].intensity;
points.push_back(pc);
});
cloud_filtered->height = 1;
cloud_filtered->width = cloud_filtered->points.size();
}
int main(int argc, char** argv)
{
string pcd_path = "000060.pcd";
// Generate pointcloud data，新建指针cloud存放点云
pcl::PointCloud<pcl::PointXYZI>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZI>);
if (pcl::io::loadPCDFile<pcl::PointXYZI>(pcd_path, *cloud) == -1)//打开点云文件。
{
PCL_ERROR("Couldn't read that pcd file\n");
return(-1);
}
vector<csf::Point> pc;
const auto& pointclouds = cloud->points;
pc.resize(cloud->size());
transform(pointclouds.begin(), pointclouds.end(), pc.begin(), [&](const auto& p)->csf::Point {
csf::Point pp;
pp.x = p.x;
pp.y = p.y;
pp.z = p.z;
return pp;
});
std::vector<int> groundIndexes, offGroundIndexes;
clothSimulationFilter(pc, groundIndexes, offGroundIndexes);
pcl::PointCloud<pcl::PointXYZI>::Ptr cloud_filtered(new pcl::PointCloud<pcl::PointXYZI>);
addPointCloud(groundIndexes, cloud,cloud_filtered);
pcl::PCDWriter writer;
writer.write<pcl::PointXYZI>("groundPointCloud.pcd", *cloud_filtered, false);
cloud_filtered->points.clear();
addPointCloud(offGroundIndexes, cloud,cloud_filtered);
writer.write<pcl::PointXYZI>("nonGroundPointCloud.pcd", *cloud_filtered, false);
pcl::visualization::CloudViewer viewer("this is a point cloud viewer haha!!");
viewer.showCloud(cloud_filtered);
while (!viewer.wasStopped())
{
}
return 0;
}二、相关配置
1.配置CSF工程
在github上下载CSF的源码，在windows下使用cmake-gui进行配置（相关操作略，可自行百度），配置好后，会出现如图所示的build文件夹。
进入所在目录下，右键打开终端，输入上面复制的地址（提前配置好github或者直接下载.zip文件解压）git clone https://github.com/jianboqi/CSF.git
命令行配置工程
cd CSF
mkdir build && cd build
cmake ..
cmake-gui配置
打开cmake-gui，输入文件路径，按下图图示配置，然后点击Tools->Configure->Generate即可.
打开VS工作，开始编译，如下图所示：
2. 新建vs工程测试
为了方便，使用cmake来配置工程，在cmakelist中添加如下，如下图所示：
最后运行结果如下图所示：
二、点云地面点滤波(Progressive Morphological Filter)算法介绍（PCL库）
1. 引言
机载LiDAR可以获取快速、低成本地获取大区域的高精度地形测量值。为了获取高精度DTM/DEM需要区分测量点中的地面点（由地面直接返回）及非地面点（建筑、车、植被）。众多学者采用了各种各样的方法来进行"点云地面点滤波"。（此篇博客中也进行了相关介绍，不再骜述）
2. Morphological Filters（形态学滤波）
2.1 膨胀/腐蚀
膨胀/腐蚀是其中的两个基础操作，通俗的说膨胀/腐蚀可以扩大/减小特征的尺寸，并以此组合为开/闭操作。针对LiDAR测量点p(x, y, z)，高程 z 值在(x, y)处对应的膨胀操作可以定义为：
式中：(xp, yp, zp) 代表p点的相邻点，w为操作窗口(可以为一维“线”也可以为二维“矩形/圆/其他形状等”)。膨胀操作完成后会输出p点在窗口w内具有最大高程值的近邻点。
与之类似的，腐蚀操作为在p点窗口w内找到具有最低高程值的近邻点。可以通过下式进行定义：
了解膨胀/腐蚀这两个基础操作之后，可以通过对其进行简单组合来来形成开/闭操作，其中开操作为先进行腐蚀再进行膨胀操作，而闭操作为先膨胀再进行腐蚀。在LiDAR数据处理中使用了“开”操作，处理效果如下图中所示：
可以从图中得知：“虚线”是先进行“腐蚀”操作所形成的表面，这个表面剔除了“树木”点，但是大型建筑物却变得不完整。“实线”是对“腐蚀”操作结果进行“膨胀”操作所形成的表面，可以看出其又恢复了大型建筑的形状。基于此，我们可以知道，“开操作”具备去除地面上的细小地物，保留大型地物的能力，这种能力对于后续处理是非常重要的。
2.2 形态学滤波
上述的“开操作”只是去除了细小地物，保留了大型地物，并没有去除所有非地面点去除，而且仅仅通过一个“开操作”也不可能实现对复杂地表的提取。因此，怎么利用好“开操作”的能力进行下一步骤的提取是进一步提升的关键。
Kilian等人提出，可以在“开操作”处理后的结果中的每一个“窗口”内找到一个“最低点”，然后此窗口内的其他点若落在“最低点”的一个高程范围内则认为是地面点。这个高程范围通常根据机载LiDAR系统扫描的精度来定义，正常为20-30cm。
此方法中有两个方面对最后的结果好坏非常重要：
滤波窗口的尺寸，如果窗口尺寸设置的比较小，则可以很好的保留地面起伏的细节，但是只能去除像汽车、树木等细小地物，而对建筑物则去除效果较差（屋顶通常被认为是地面）。相反，若窗口尺寸设置的较大，则会过度去除一些“地面点”，例如，一条道路与一条小水沟相邻，若窗口尺寸大于道路的宽度，则道路可能就会被认为是非地面点（因为小水沟中的点高程较低，会被认为是窗口内的最低点，而道路点较高，被判断为非地面点）。此外，一些局部的小山丘、沙丘都极可能被“切除”。
建筑与树木在特定/局部区域的分布。
注：一个最理想的情况是我们可以设置一个“窗口”，这个“窗口”的尺寸可以足够小，能够保留地面细节。同时，还可以足够大，能够去除建筑、汽车、树木等地物。但是，这种理想情况在实际数据集中国并不存在。
为了解决这一问题，Kilian等人继续提出了可以通过改变窗口大小来多次进行滤波。每个点都被赋予一个与窗口大小相关的权重，窗口尺寸越大，点的权重越高。这种方法虽然得到了更好一些的效果，但是没有从"point level"进行区分地面点与非地面点。（"point level"区分的地面点与非地面点之后可以通过插值的方法使得DEM/DTM的生成效果更好。）
3. Progressive Morphological Filters
由上述2.2节中的分析可知，传统的形态学滤波很难通过一个固定大小的窗口去检测出各种尺寸变化的不同地物。这一问题可以通过逐渐改变窗口大小来解决。
如下图中所示，首先使用一个尺寸为l1的窗口来对原始数据进行开操作。由图中的“虚线”可以看出树木等尺寸小于l1的地物被去除，且地形特征中小于l1的部分被“切除”（山丘顶部高程被替换成了l1中的最小值），但是，尺寸大于l1的建筑物被保留了下来。接着，进行下一次迭代，窗口尺寸变为l2，对上一次的处理结果进行开操作处理，结果从“实线”中可以看出，l2大于建筑的尺寸，所以建筑也被去除，但同时山丘顶部被“切除”的范围更大。
需要注意的是：通过逐渐增加窗口尺寸解决了去除不同大小地物的问题，但是又引入了"山丘"顶部等小于窗口尺寸的地形特征部分被“切除”的问题。
为了解决这一新出现的问题，可以通过引入一个高度差阈值来解决。建筑屋顶和地面点之间的高程差是“突变”（abrupt change），而地面高程是“渐变”（gradual change）。因此，二者之间高程变化中的明显区别可以帮助我们进行区分。假设dhp,1代表原始LiDAR测量值与在任意给定p点处第一次迭代表面之间的高程差，dhT,1代表高程差阈值，则如果dhp,1 ≤ dhT,1点p就被认为是地面点，反之如果dhp,1 > dhT,1就认为点p是一个非地面点。此后，令dhmax(t),1为当前迭代中初始地面点与滤波表面之间差值的最大值，则如果选取的dhT,1 > dhmax(t),1则所有的测量值都会保留。
在第二次迭代中假设上一次滤波表面和本次滤波表面的最大高程差为dhmax(t),2，则如果dhmax(t),2 < dhT,2,则高程差值在dhmax(t),2范围内的地面点都会被保留。类似的，假设在上次迭代和本次迭代之间建筑高程差值最小为dhmin(b),2（通常近似为建筑的高度），如果dhmin(b),2 > dhT,2，则建筑就会被移除。
通常设置dhT,k为研究区域第k次迭代中建筑物的最矮高程值。以dhT,k作为阈值，对于第k次迭代中的任意点p如果dhp,k < dhT,k则将其标记为地面点，反之为非地面点。通过这种方式，不同尺寸的建筑物（树）可以随着迭代窗口尺寸的增加逐步被去除。
综上所述，Progressive Morphological Filters的详细流程如下图所示：
我们可以对上图总结以下四个步骤：
加载不规则点云，划分为规则网格，在每个网格中选取高程最低点（如果网格中没有点则根据最近邻点进行插值），构建最小高程表面。
使用输入的初始滤波窗口尺寸、1）中获取的最小高程表面作为第一次迭代的参数进行第一次迭代。随后，在后续的迭代中，以前一次获取的滤波表面及3）中计算的滤波窗口尺寸作为输入。每次迭代的输出有两部分：a) 形态学滤波后得到的更加光滑的表面；b) 基于不同阈值所检测出来的非地面点。
计算新的滤波窗口尺寸及不同的高程插值阈值。重复步骤2)、步骤3)直到窗口尺寸大于预设的最大窗口。
基于去除非地面测量值的数据进行生成DEM/DTM。
注意：每一次迭代中的“开操作”实际都是作用在步骤1）所划分网格中的点，所以Progressive Morphological
Filters是"point level"来对LiDAR测量值进行滤波处理的。
3.1 参数计算(窗口尺寸/高程差阈值)
在上述步骤3）中我们要变化窗口尺寸 wk和高程差阈值 dhT,k两个参数的值，以进行下一次迭代，那么这两个值是怎么计算的呢？
3.1.1 窗口尺寸首先是窗口尺寸 wk有两种计算方式，第一种是：
式中，k为迭代次数，b是初始窗口大小（由用户进行输入），最后+1是为了保证 wk为一个奇数，窗口对称。但是，如果一个研究区域具有非常大的地物，这种增加窗口尺寸速度太慢则会耗费较多时间。因此，可以使用第二种方式，通过指数增长来改变窗口大小，计算如下式：
同样的，式中k为迭代次数，b是初始窗口大小（由用户进行输入），这种方式的增长速度较第一种方式快很多。
3.1.2 高程差阈值
高程差阈值与研究区域的地形坡度密不可分，因此可以通过下式进行计算：
式中，dh0为初始高程差阈值，s为坡度，c为格网大小，dhmax为最大高程差阈值。在城市区域，树木、汽车相对于建筑的尺寸小很多，所以通常是最后滤除建筑，最大高程差阈值dhmax可以设置为一个固定值（如最矮建筑物高度）。而在山区，主要的非地面点为植被，所以并没有必要设置固定的最大高程差阈值dhmax，于是dhmax通常被设置为测区内的最大高程差。
此外，坡度s通过第k次迭代的最大高程差dhmax(t),k，以及窗口尺寸wk进行计算，如下式所示：
3.2 参数输入/输出
3.2.1 参数输入
原始LiDAR点云数据，每个点都由(x, y, z)进行表示
划分格网大小c 参数b（计算窗口尺寸）
最大窗口尺寸（判断是否停止迭代）
地形坡度s
初始高程差阈值dh0
最大高程差值dhmax
3.2.1 参数输出
地面点
非地面点
3.3 PCL官方示例代码
#include <iostream>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/filters/extract_indices.h>
#include <pcl/segmentation/progressive_morphological_filter.h>
int
main (int argc, char** argv)
{
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PointIndicesPtr ground (new pcl::PointIndices);
// Fill in the cloud data
pcl::PCDReader reader;
// Replace the path below with the path where you saved your file
reader.read<pcl::PointXYZ> ("samp11-utm.pcd", *cloud);
std::cerr << "Cloud before filtering: " << std::endl;
std::cerr << *cloud << std::endl;
// Create the filtering object
pcl::ProgressiveMorphologicalFilter<pcl::PointXYZ> pmf;
pmf.setInputCloud (cloud);
pmf.setMaxWindowSize (20);
pmf.setSlope (1.0f);
pmf.setInitialDistance (0.5f);
pmf.setMaxDistance (3.0f);
pmf.extract (ground->indices);
// Create the filtering object
pcl::ExtractIndices<pcl::PointXYZ> extract;
extract.setInputCloud (cloud);
extract.setIndices (ground);
extract.filter (*cloud_filtered);
std::cerr << "Ground cloud after filtering: " << std::endl;
std::cerr << *cloud_filtered << std::endl;
pcl::PCDWriter writer;
writer.write<pcl::PointXYZ> ("samp11-utm_ground.pcd", *cloud_filtered, false);
// Extract non-ground returns
extract.setNegative (true);
extract.filter (*cloud_filtered);
std::cerr << "Object cloud after filtering: " << std::endl;
std::cerr << *cloud_filtered << std::endl;
writer.write<pcl::PointXYZ> ("samp11-utm_object.pcd", *cloud_filtered, false);
return (0);
}接下来的很多项目都需要使用3D视觉，感觉有必要好好学习3D点云相关算法，然后买深蓝学院的三维点云处理课程，还在学习中，新的一年，加油！
（一）视频课程来了！
自动驾驶之心为大家汇集了毫米波雷达视觉融合、高精地图、BEV感知、传感器标定、传感器部署、自动驾驶协同感知、语义分割、自动驾驶仿真、L4感知、决策规划、轨迹预测等多个方向学习视频，欢迎大家自取（扫码进入学习）
（扫码学习最新视频）
视频官网：www.zdjszx.com
（二）国内首个自动驾驶学习社区
近1000人的交流社区，和20+自动驾驶技术栈学习路线，想要了解更多自动驾驶感知（分类、检测、分割、关键点、车道线、3D目标检测、Occpuancy、多传感器融合、目标跟踪、光流估计、轨迹预测）、自动驾驶定位建图（SLAM、高精地图）、自动驾驶规划控制、领域技术方案、AI模型部署落地实战、行业动态、岗位发布，欢迎扫描下方二维码，加入自动驾驶之心知识星球，这是一个真正有干货的地方，与领域大佬交流入门、学习、工作、跳槽上的各类难题，日常分享论文+代码+视频，期待交流！
（三）【自动驾驶之心】全栈技术交流群
自动驾驶之心是首个自动驾驶开发者社区，聚焦目标检测、语义分割、全景分割、实例分割、关键点检测、车道线、目标跟踪、3D目标检测、BEV感知、多传感器融合、SLAM、光流估计、深度估计、轨迹预测、高精地图、NeRF、规划控制、模型部署落地、自动驾驶仿真测试、产品经理、硬件配置、AI求职交流等方向；
添加汽车人助理微信邀请入群
备注：学校/公司+方向+昵称
6538
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
