python - 如何有效地将点云数据的大型 numpy 数组转换为下采样的二维数组？

Question

我有一个大型无序激光雷达点云数据的 numpy 数组，形状为 [num_points, 3]，它们是每个点的 XYZ 坐标。我想将其下采样为平均高度值的 2D 网格 - 为此，我想将数据拆分为 5x5 X-Y 箱并计算每个箱中的平均高度值(Z 坐标)。

有人知道有什么快速/有效的方法吗？

当前代码:

import numpy as np
from open3d import read_point_cloud

resolution = 5

# Code to load point cloud and get points as numpy array
pcloud = read_point_cloud(params.POINT_CLOUD_DIR + "Part001.pcd")
pcloud_np = np.asarray(pcloud.points)

# Code to generate example dataset
pcloud_np = np.random.uniform(0.0, 1000.0, size=(1000,3))

# Current (inefficient) code to quantize into 5x5 XY 'bins' and take mean Z values in each bin
pcloud_np[:, 0:2] = np.round(pcloud_np[:, 0:2]/float(resolution))*float(resolution) # Round XY values to nearest 5

num_x = int(np.max(pcloud_np[:, 0])/resolution)
num_y = int(np.max(pcloud_np[:, 1])/resolution)

mean_height = np.zeros((num_x, num_y))

# Loop over each x-y bin and calculate mean z value 
x_val = 0
for x in range(num_x):
    y_val = 0
    for y in range(num_y):
        height_vals = pcloud_np[(pcloud_np[:,0] == float(x_val)) & (pcloud_np[:,1] == float(y_val))]
        if height_vals.size != 0:
            mean_height[x, y] = np.mean(height_vals)
        y_val += resolution
    x_val += resolution

Answer 1

这里是在扁平二维网格上使用 np.bincount 习惯用法的建议。我还冒昧地对原始代码添加了一些小修复:

import numpy as np
#from open3d import read_point_cloud

resolution = 5

# Code to load point cloud and get points as numpy array
#pcloud = read_point_cloud(params.POINT_CLOUD_DIR + "Part001.pcd")
#pcloud_np = np.asarray(pcloud.points)

# Code to generate example dataset
pcloud_np = np.random.uniform(0.0, 1000.0, size=(1000,3))

def f_op(pcloud_np, resolution):
    # Current (inefficient) code to quantize into 5x5 XY 'bins' and take mean Z values in each bin
    pcloud_np[:, 0:2] = np.round(pcloud_np[:, 0:2]/float(resolution))*float(resolution) # Round XY values to nearest 5

    num_x = int(np.max(pcloud_np[:, 0])/resolution) + 1
    num_y = int(np.max(pcloud_np[:, 1])/resolution) + 1

    mean_height = np.zeros((num_x, num_y))

    # Loop over each x-y bin and calculate mean z value 
    x_val = 0
    for x in range(num_x):
        y_val = 0
        for y in range(num_y):
            height_vals = pcloud_np[(pcloud_np[:,0] == float(x_val)) & (pcloud_np[:,1] == float(y_val)), 2]
            if height_vals.size != 0:
                mean_height[x, y] = np.mean(height_vals)
            y_val += resolution
        x_val += resolution

    return mean_height

def f_pp(pcloud_np, resolution):
    xy = pcloud_np.T[:2]
    xy = ((xy + resolution / 2) // resolution).astype(int)
    mn, mx = xy.min(axis=1), xy.max(axis=1)
    sz = mx + 1 - mn
    flatidx = np.ravel_multi_index(xy-mn[:, None], sz)
    histo = np.bincount(flatidx, pcloud_np[:, 2], sz.prod()) / np.maximum(1, np.bincount(flatidx, None, sz.prod()))
    return (histo.reshape(sz), *(xy * resolution))

res_op = f_op(pcloud_np, resolution)
res_pp, x, y = f_pp(pcloud_np, resolution)

from timeit import timeit

t_op = timeit(lambda:f_op(pcloud_np, resolution), number=10)*100
t_pp = timeit(lambda:f_pp(pcloud_np, resolution), number=10)*100

print("results equal:", np.allclose(res_op, res_pp))
print(f"timings (ms) op: {t_op:.3f} pp: {t_pp:.3f}")

示例输出:

results equal: True
timings (ms) op: 359.162 pp: 0.427

加速近 1000 倍。