python實現CSF地面點濾波算法原理解析

Posted on by

一、算法原理

佈料模擬濾波處理流程：
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、最後得到的結果：

可以看出，非地面點中不能提取到路緣石。

到此這篇關於python實現CSF地面點濾波的文章就介紹到這瞭,更多相關python地面點濾波內容請搜索WalkonNet以前的文章或繼續瀏覽下面的相關文章希望大傢以後多多支持WalkonNet！

推薦閱讀：