点云下采样/抽稀python-pcl：pcl::VoxelGrid::applyFilter

Posted 2022-11-23 程序媛一枚~

这篇博客将介绍如何使用python-pcl对点云las/laz进行下采样/抽稀，可以根据设置的参数决定下采样到元数据的90%，80%，60%，或者40%，不会完整的按这个梯度递减，但参数rate顺序下降，基本能大致达到这个概率。

下采样设置的voxelGridFilter.set_leaf_size(rate,rate,rate) 值越大，最后保留的点云越少。

虽然用的是python-pcl的api调用下采样算法，实质上调用的仍然是C++的VoxelGridFilter算法。

import pcl
import math
import time
from laspy.file import File
import numpy as np


# 初始文件路径  输出文件路径  抽稀参数（单位m）
def cx(filePath, outputPath, rate):
    end1 = time.time()
    f = File(filePath, mode='r')
    total = len(f.points)
    print('points： ', total)

    # 获取偏移量
    offset = f.header.offset
    x0 = offset[0]
    y0 = offset[1]
    z0 = offset[2]
    inFile = np.vstack((f.x - x0, f.y - y0, f.z - z0, f.intensity)).transpose()
    sp = math.ceil(len(inFile) / 1000000)

    cloud = pcl.PointCloud_PointXYZI()
    m2 = []
    for i in range(0, sp):
        end = (i + 1) * 1000000
        start = i * 1000000
        clPoints = inFile[start:end]
        cloud.from_array(np.array(clPoints, dtype=np.float32))

        # 抽稀
        sor = cloud.make_voxel_grid_filter()
        sor.set_leaf_size(rate, rate, rate)
        cloud_filtered = sor.filter()
        for i in range(0, cloud_filtered.size):
            m2.append(
                (cloud_filtered[i][0] + x0, cloud_filtered[i][1] + y0, cloud_filtered[i][2] + z0, cloud_filtered[i][3]))

    x = [k[0] for k in m2]
    y = [k[1] for k in m2]
    z = [k[2] for k in m2]
    c = [k[3] for k in m2]

    print('【outputPath】', outputPath)
    cxlen = len(x)
    outFile = File(outputPath, mode='w', header=f.header)
    outFile.x = np.array(x)
    outFile.y = np.array(y)
    outFile.z = np.array(z)
    outFile.intensity = np.array(c)
    outFile.header.set_pointrecordscount(cxlen)
    outFile.close()

    print('cx ' + str(rate) + ' success...')
    print('total: ', total, 'cx: ', cxlen)
    print('percent: :.2%'.format(cxlen / total))
    end2 = time.time()
    print("cx 耗时:%.2f秒" % (end2 - end1))
    return cxlen


def main():
    cx('D:/project/las/1001140020191217.las',
       'D:/project/las/1001140020191217_cx.las', 0.03)
    cx('D:/project/las/1001140020191217.las',
       'D:/project/las/1001140020191217_cx_005.las', 0.05)



if __name__ == "__main__":
    main()