python實(shí)現(xiàn)CSF地面點(diǎn)濾波算法原理解析

import laspy
#============讀取las格式的點(diǎn)云===========
inFile = laspy.file.File(r"40m1.las", mode='r') # 讀取點(diǎn)云
print('X,Y,Z',inFile.x,inFile.y,inFile.z) # 輸出點(diǎn)云坐標(biāo)
print('點(diǎn)云個(gè)數(shù):',len(inFile)) #讀取點(diǎn)云個(gè)數(shù)
#============保存點(diǎn)云為las文件===========
h = inFile.header
outFile = laspy.file.File('666.las', mode = "w", header=h)
points = inFile #對(duì)點(diǎn)云進(jìn)行的相關(guān)操作
outFile.points = points
outFile.close() #關(guān)閉文件完成保存

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()
# 參數(shù)設(shè)置
csf.params.bSloopSmooth = False    #粒子設(shè)置為不可移動(dòng)
csf.params.cloth_resolution = 0.1  #布料網(wǎng)格分辨率
csf.params.rigidness = 3  #布料剛性參數(shù)
csf.params.time_step = 0.65
csf.params.class_threshold = 0.03 #點(diǎn)云與布料模擬點(diǎn)的距離閾值
csf.params.interations = 500      #最大迭代次數(shù)
# more details about parameter: http://ramm.bnu.edu.cn/projects/CSF/download/
csf.setPointCloud(xyz)
ground = CSF.VecInt()  # 地面點(diǎn)索引列表
non_ground = CSF.VecInt() # 非地面點(diǎn)索引列表
csf.do_filtering(ground, non_ground) # 執(zhí)行濾波
#============保存為las文件==========
outFile = laspy.file.File(r"non_ground.las",
                          mode='w', header=inFile.header)
outFile.points = points[non_ground] # 提取非地面點(diǎn)保存到las
outFile.close() # 關(guān)閉文件夾

a=xyz[ground]
b=xyz[non_ground]
#=============可視化===============
def view_cloud(a, b):
    pcd = o3d.geometry.PointCloud()
    # =====numpy轉(zhuǎn)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='提取結(jié)果')
    o3d.visualization.draw_geometries([pcd1],window_name='非地面點(diǎn)')
    o3d.visualization.draw_geometries([pcd],window_name='地面點(diǎn)')
view_cloud(a,b)

import open3d as o3d
import CSF
import numpy as np

pc = o3d.io.read_point_cloud("數(shù)據(jù)//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)#保存點(diǎn)云
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='提取結(jié)果')
    o3d.visualization.draw_geometries([pcd1],window_name='非地面點(diǎn)')
    o3d.visualization.draw_geometries([pcd],window_name='地面點(diǎn)')
view_cloud(a,b)