LiDAR点云粗糙度计算提速问询：并行处理失效后的优化方案

问题背景 #

需要计算LiDAR点云的表面粗糙度，方法为：对每个点的邻域点集拟合最佳拟合平面，再计算该点到平面的距离。但邻域半径增大时计算耗时剧增，尝试用concurrent.futures并行处理后速度反而比串行更慢，寻求可行的提速方法。

现有代码 #

import numpy as np
import laspy
from concurrent.futures import ProcessPoolExecutor
from scipy.spatial import cKDTree as KDTree
from skspatial.objects import Plane

def process_point(point_idx, point, lidar_data, neighborhood_radius, tree):
    neighbor_indices = tree.query_ball_point(point, neighborhood_radius)
    neighbors = lidar_data[neighbor_indices]
    
    if neighbors.shape[0] >= 3:
        plane = Plane.best_fit(neighbors)
        # Distance calculation remains the same
        # Return roughness value and index
        return point_idx, calculate_distance_to_plane(point, plane)
    else:
        return point_idx, np.nan

# Parallel processing wrapper function
def calculate_roughness_parallel(lidar_data, neighborhood_radius, tree):
    roughness_values = np.zeros(len(lidar_data))
    with ProcessPoolExecutor() as executor:
        futures = [executor.submit(process_point, point_idx, point, lidar_data, neighborhood_radius, tree)
                   for point_idx, point in enumerate(lidar_data)]
        for future in futures:
            point_idx, roughness = future.result()
            roughness_values[point_idx] = roughness
    return roughness_values


if __name__ == "__main__":
    las = laspy.read("las0.laz")
    in_point = np.vstack((las.x, las.y, las.z)).transpose()
    tree = KDTree(in_point)
    roughness_result = calculate_roughness_parallel(in_point, 1, tree)

可行提速方法 #

1. 替换并行方式：用线程池替代进程池 #

ProcessPoolExecutor的进程间通信开销极大，尤其是传递大量点云数据时。改用ThreadPoolExecutor（线程池），利用多线程处理IO密集型的KDTree查询，同时结合NumPy向量化运算减少GIL锁的影响：

from concurrent.futures import ThreadPoolExecutor

def calculate_roughness_threadpool(lidar_data, neighborhood_radius, tree):
    roughness_values = np.zeros(len(lidar_data))
    with ThreadPoolExecutor() as executor:
        futures = [executor.submit(process_point, point_idx, point, lidar_data, neighborhood_radius, tree)
                   for point_idx, point in enumerate(lidar_data)]
        for future in futures:
            point_idx, roughness = future.result()
            roughness_values[point_idx] = roughness
    return roughness_values

2. 重写平面拟合逻辑：用NumPy替代skspatial #

skspatial.objects.Plane.best_fit的底层实现存在额外的类实例化开销，改用NumPy直接计算最佳拟合平面，大幅提升计算效率：

def calculate_best_fit_plane(neighbors):
    # 计算邻域点中心
    centroid = np.mean(neighbors, axis=0)
    # 去中心化处理
    centered = neighbors - centroid
    # 协方差矩阵SVD分解求法向量
    _, _, vh = np.linalg.svd(centered.T @ centered)
    normal = vh[-1]
    # 计算平面方程参数d
    d = -np.dot(normal, centroid)
    return normal, d

def calculate_distance_to_plane(point, normal, d):
    # 点到平面距离公式
    return abs(np.dot(normal, point) + d) / np.linalg.norm(normal)

# 修改process_point函数适配新逻辑
def process_point(point_idx, point, lidar_data, neighborhood_radius, tree):
    neighbor_indices = tree.query_ball_point(point, neighborhood_radius)
    neighbors = lidar_data[neighbor_indices]
    
    if neighbors.shape[0] >= 3:
        normal, d = calculate_best_fit_plane(neighbors)
        return point_idx, calculate_distance_to_plane(point, normal, d)
    else:
        return point_idx, np.nan

3. 批量邻域查询：减少循环调用开销 #

利用KDTree的query_ball_point支持批量输入的特性，一次性获取所有点的邻域索引，避免循环调用的额外开销：

def calculate_roughness_batch(lidar_data, neighborhood_radius, tree):
    roughness_values = np.full(len(lidar_data), np.nan)
    # 批量获取所有点的邻域索引
    all_neighbor_indices = tree.query_ball_point(lidar_data, neighborhood_radius)
    
    for point_idx, neighbor_indices in enumerate(all_neighbor_indices):
        neighbors = lidar_data[neighbor_indices]
        if len(neighbors) >=3:
            normal, d = calculate_best_fit_plane(neighbors)
            roughness_values[point_idx] = calculate_distance_to_plane(lidar_data[point_idx], normal, d)
    return roughness_values

4. 降采样预处理：减少计算量 #

若精度要求允许，先对原始点云进行体素降采样，在保留整体特征的前提下减少需要处理的点数量：

def voxel_downsample(points, voxel_size):
    # 计算每个点所属的体素索引
    voxel_indices = np.floor(points / voxel_size).astype(int)
    # 按体素分组，取均值作为降采样点
    unique_indices, indices = np.unique(voxel_indices, axis=0, return_inverse=True)
    downsampled = np.zeros((len(unique_indices), 3))
    np.add.at(downsampled, indices, points)
    counts = np.bincount(indices)
    downsampled /= counts[:, np.newaxis]
    return downsampled

# 主函数中使用降采样
if __name__ == "__main__":
    las = laspy.read("las0.laz")
    in_point = np.vstack((las.x, las.y, las.z)).transpose()
    # 体素大小可根据需求调整
    downsampled_point = voxel_downsample(in_point, voxel_size=0.5)
    tree = KDTree(downsampled_point)
    roughness_result = calculate_roughness_batch(downsampled_point, 1, tree)

5. Numba JIT编译：核心函数加速 #

用numba对平面拟合、距离计算等核心函数进行JIT编译，将Python代码转为机器码执行，大幅提升计算速度：

from numba import jit

@jit(nopython=True)
def calculate_best_fit_plane_numba(neighbors):
    centroid = np.mean(neighbors, axis=0)
    centered = neighbors - centroid
    cov_matrix = np.dot(centered.T, centered)
    _, _, vh = np.linalg.svd(cov_matrix)
    normal = vh[-1]
    d = -np.dot(normal, centroid)
    return normal, d

@jit(nopython=True)
def calculate_distance_to_plane_numba(point, normal, d):
    return abs(np.dot(normal, point) + d) / np.linalg.norm(normal)

内容的提问来源于stack exchange，提问作者Purple_Ad