当前位置 博文首页 > python实现CSF地面点滤波算法原理解析

    python实现CSF地面点滤波算法原理解析

    作者:点云侠 时间:2021-09-19 17:45

    目录
    • 一、算法原理
    • 二、读取las点云
    • 三、算法源码
    • 四、结果展示
    • 五、CloudCompare实现

    一、算法原理

    布料模拟滤波处理流程:
    1)利用点云滤波算法或者点云处理软件滤除异常点;
    2)将激光雷达点云倒置;
    3)设置模拟布料,设置布料网格分辨率 G R GR GR,确定模拟粒子数。布料的位置设置在点云最高点以上;
    4)将布料模拟点和雷达点投影到水平面,为每个布料模拟点找到最相邻的激光点的高度值,将高度值设置为 I H V IHV IHV;
    5)布料例子设置为可移动,布料粒子首先受到重力作用,当粒子高度 C H V CHV CHV小于 I H V IHV IHV时,将粒子高度设置为 I H V IHV IHV;粒子设置为不可移动;
    6)计算布料粒子之间的内力作用,根据设置的布料刚性参数,调整布料粒子之间的相对位置;
    7)重复进行5)和6)计算,迭代次数达到设置的最大迭代次数;
    8)计算激光雷达点与对应布料模拟点的距离,距离小于阈值标记为地面点,距离大于阈值标记为非地面点。

    点云地面点滤波(Cloth Simulation Filter, CSF)“布料”滤波算法介绍

    二、读取las点云

    参考链接: python读取las
    1、GitHub: laspy
    2、基础教程:Laspy: Documentation
    3、安装:pip install laspy
    4、使用example:

    import laspy
    #============读取las格式的点云===========
    inFile = laspy.file.File(r"40m1.las", mode='r') # 读取点云
    print('X,Y,Z',inFile.x,inFile.y,inFile.z) # 输出点云坐标
    print('点云个数:',len(inFile)) #读取点云个数
    #============保存点云为las文件===========
    h = inFile.header
    outFile = laspy.file.File('666.las', mode = "w", header=h)
    points = inFile #对点云进行的相关操作
    outFile.points = points
    outFile.close() #关闭文件完成保存

    三、算法源码

    1、算法细节:CSF
    2、源码获取:https://github.com/jianboqi/CSF
    3、源码编译:下载源代码。在python文件夹下:
    python setup.py build
    python setup.py install
    4、读取las并可视化算法结果

    import laspy
    import CSF
    import numpy as np
    import open3d as o3d
    #============读取las文件=============
    inFile = laspy.file.File(r"40m1.las", mode='r') # read a las file
    points = inFile.points
    xyz = np.vstack((inFile.x, inFile.y, inFile.z)).transpose() # extract x, y, z and put into a list
    #============布料模拟滤波============
    csf = CSF.CSF()
    # 参数设置
    csf.params.bSloopSmooth = False    #粒子设置为不可移动
    csf.params.cloth_resolution = 0.1  #布料网格分辨率
    csf.params.rigidness = 3  #布料刚性参数
    csf.params.time_step = 0.65
    csf.params.class_threshold = 0.03 #点云与布料模拟点的距离阈值
    csf.params.interations = 500      #最大迭代次数
    # more details about parameter: http://ramm.bnu.edu.cn/projects/CSF/download/
    csf.setPointCloud(xyz)
    ground = CSF.VecInt()  # 地面点索引列表
    non_ground = CSF.VecInt() # 非地面点索引列表
    csf.do_filtering(ground, non_ground) # 执行滤波
    #============保存为las文件==========
    outFile = laspy.file.File(r"non_ground.las",
                              mode='w', header=inFile.header)
    outFile.points = points[non_ground] # 提取非地面点保存到las
    outFile.close() # 关闭文件夹
    
    a=xyz[ground]
    b=xyz[non_ground]
    #=============可视化===============
    def view_cloud(a, b):
        pcd = o3d.geometry.PointCloud()
        # =====numpy转point=======
        pcd.points = o3d.utility.Vector3dVector(a)
    
        pcd1 = o3d.geometry.PointCloud()
    
        pcd1.points = o3d.utility.Vector3dVector(b)
        #=======自定义颜色========
        pcd.paint_uniform_color([0, 1, 0])
        pcd1.paint_uniform_color([1, 0, 0])
        o3d.visualization.draw_geometries([pcd, pcd1],window_name='提取结果')
        o3d.visualization.draw_geometries([pcd1],window_name='非地面点')
        o3d.visualization.draw_geometries([pcd],window_name='地面点')
    view_cloud(a,b)
    

    5、读取pcd文件并可视化结果

    import open3d as o3d
    import CSF
    import numpy as np
    
    pc = o3d.io.read_point_cloud("数据//100m1.pcd")
    xyz = np.asarray(pc.points)
    csf = CSF.CSF()
    # prameter settings
    csf.params.bSloopSmooth = False
    csf.params.cloth_resolution = 0.1
    csf.params.rigidness = 3
    csf.params.time_step = 0.65
    csf.params.class_threshold = 0.03
    csf.params.interations = 500
    # more details about parameter: http://ramm.bnu.edu.cn/projects/CSF/download/
    csf.setPointCloud(xyz)
    ground = CSF.VecInt()  # a list to indicate the index of ground points after calculation
    non_ground = CSF.VecInt() # a list to indicate the index of non-ground points after calculation
    csf.do_filtering(ground, non_ground) # do actual filtering.
    
    # o3d.io.write_point_cloud("trans_of_source.pcd", non_ground)#保存点云
    a=xyz[ground]
    b=xyz[non_ground]
    def view_cloud(a, b):
        pcd = o3d.geometry.PointCloud()
        # From numpy to Open3D
        pcd.points = o3d.utility.Vector3dVector(a)
    
        pcd1 = o3d.geometry.PointCloud()
        # From numpy to Open3D
        pcd1.points = o3d.utility.Vector3dVector(b)
    
        pcd.paint_uniform_color([0, 1, 0])
        pcd1.paint_uniform_color([1, 0, 0])
        o3d.visualization.draw_geometries([pcd, pcd1],window_name='提取结果')
        o3d.visualization.draw_geometries([pcd1],window_name='非地面点')
        o3d.visualization.draw_geometries([pcd],window_name='地面点')
    view_cloud(a,b)
    

    四、结果展示

    五、CloudCompare实现

    1、加载点云数据,点击Plugins中的CSF Filter功能

    2、弹出如下窗口:



     图中:Cloth resolution:是指用于覆盖地形的布的网格大小(单位与点云的单位相同)。你设置的布分辨率越大,你得到的DTM就越粗糙;Max iterations:是指地形仿真的最大迭代次数。500对大多数场景来说都足够了。Classification threshold:是指根据点与模拟地形之间的距离,将点云划分为地面和非地面部分的阈值。0.5适用于大多数场景
      这里的网格分辨率和距离阈值最小只能设置为10cm,地面10cm的范围默认是地面点,精确度不如自己代码实现中的高。
    3、最后得到的结果:

    可以看出,非地面点中不能提取到路缘石。

    jsjbwy
    下一篇:没有了