并行处理激光雷达点云以找到最佳拟合平面

Question

我想计算激光雷达点云的表面粗糙度。为此，我需要找到每个点的邻居的最佳拟合平面并计算从平面到该点的距离。当半径增加时，这个过程需要很长时间。我尝试使用并发 future，但这使得它比没有并行处理的代码更慢。

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)

还有其他方法可以让它更快吗？

Answer 1

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

def process_batch(batch, lidar_data, neighborhood_radius, tree):
    results = []
    for point_idx, point in batch:
        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 = calculate_distance_to_plane(point, plane)
            results.append((point_idx, distance))
        else:
            results.append((point_idx, np.nan))
    return results

def calculate_roughness_parallel(lidar_data, neighborhood_radius, tree, batch_size=100):
    roughness_values = np.zeros(len(lidar_data))
    batches = [(lidar_data[i:i + batch_size], neighborhood_radius, tree) 
               for i in range(0, len(lidar_data), batch_size)]
    
    with ProcessPoolExecutor() as executor:
        futures = [executor.submit(process_batch, [(idx, point) for idx, point in enumerate(batch)], 
                                   lidar_data, neighborhood_radius, tree) 
                   for batch in batches]
                   
        for future in futures:
            for point_idx, roughness in future.result():
                roughness_values[point_idx] = roughness
    
    return roughness_values

if __name__ == "__main__":
    las = laspy.read("path/to/your/lasfile.las")
    lidar_data = np.vstack((las.x, las.y, las.z)).transpose()
    tree = KDTree(lidar_data)
    roughness_result = calculate_roughness_parallel(lidar_data, 1, tree)