如何在 numpy 中有效地从 Nx3 点云阵列获取表面点？

Question

假设我有一个点云的 numpy 数组（形状：Nx3，每行是 (x, y, z)）。我想从点云生成高度图（投影到 xy 平面）。对于 xy 平面上的每个网格，我只想保留投影到该网格的所有点的最大 z 值。

例如，

A=np.array([[0, 0, 1],
            [0, 0, 1.5],
            [0, 0, 2.0],
            [1, 1, 1],
            [1, 2, 1],
            [1, 1, 3]])

那么我期望的输出是

B=np.array([[0, 0, 2.0],
             [1, 1, 3],
             [1, 2, 1]])

如何在 numpy 中高效地执行此操作？谢谢。

Answer 1

方法 1a

这可以用np.lexsort, np.unique and np.ufunc.reduceat来完成：

A = A[np.lexsort((A[:, 0], A[:, 1]))]
_, idx = np.unique(A[:,:2], return_index = True, axis=0)
output = np.maximum.reduceat(A, idx)

方法 1b

我们还可以以一种稍微有效的方式提高 reduceat 的速度：

A = A[np.lexsort((A[:, 0], A[:, 1]))]
u, idx = np.unique(A[:, :2], return_index = True, axis=0)
return np.c_[u, np.maximum.reduceat(A[:,2], idx)]

方法二

如果你想要一个更简单的方法，你也可以使用 numpy_indexed，它写在 numpy 之上，可以用更少的脚本解决分组问题：

import numpy_indexed as npi
_, idx = npi.group_by(A[:, :2]).argmax(A[:, 2])
output = A[idx]

请注意，它优于以前的方法，这表明它可以进一步优化。

方法 3

某些 pandas 方法比 numpy 工作得更快。看来您的方法很幸运：

import pandas as pd
df = pd.DataFrame(A)
return df.loc[df.groupby([0,1])[2].idxmax()].values

输出

所有输出都是：np.array([[0. 0. 2.], [1. 1. 3.], [1. 2. 1.]])，除了导致 np.array([[0. 0. 2.], [1. 2. 1.], [1. 1. 3.]])

的 npi 方法的情况

更新

如果对降维结果的展平数组执行相同的算法，则可以进一步优化它。 numexpr 包在这里可以达到极速。新方法的名称是：approach1a_on1D、approach1b_on1D、approach2_on1D.

import numexpr as ne

def reduct_dims(cubes):
    m0, m1 = np.min(cubes[:,:2], axis=0)
    M0 = np.max(cubes[:,0], axis=0)
    s0  = M0 - m0 + 1
    d = {'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2],
         's0':s0,'m0':m0, 'm1':m1}
    return ne.evaluate('c0-m0+(c1-m1)*s0', d)

def approach1a(A):
    A = A[np.lexsort((A[:, 0], A[:, 1]))]
    _, idx = np.unique(A[:, :2], return_index = True, axis=0)
    return np.maximum.reduceat(A, idx)

def approach1b(A):
    A = A[np.lexsort((A[:, 0], A[:, 1]))]
    u, idx = np.unique(A[:, :2], return_index = True, axis=0)
    return np.c_[u, np.maximum.reduceat(A[:,2], idx)]

def approach2(A):
    _, idx = npi.group_by(A[:, :2]).argmax(A[:, 2])
    return A[idx]

def approach3(A):
    df = pd.DataFrame(A)
    return df.loc[df.groupby([0,1])[2].idxmax()].values

def approach1a_on1D(A):
    A_as_1D = reduct_dims(A)
    argidx = np.argsort(A_as_1D)
    A, A_as_1D = A[argidx], A_as_1D[argidx] #sort both arrays
    _, idx = np.unique(A_as_1D, return_index = True)
    return np.maximum.reduceat(A, idx)

def approach1b_on1D(A):
    A_as_1D = reduct_dims(A)
    argidx = np.argsort(A_as_1D)
    A, A_as_1D = A[argidx], A_as_1D[argidx]
    _, idx = np.unique(A_as_1D, return_index = True)
    return np.c_[A[:,:2][idx], np.maximum.reduceat(A[:,2], idx)]

def approach2_on1D(A):
    A_as_1D = reduct_dims(A)
    _, idx = npi.group_by(A_as_1D).argmax(A[:,2])
    return A[idx]

%timeit reduct_dims(cubes)

160 ms ± 7.44 ms per loop (mean ± std. dev. of 7 runs, 10 
loops each)

使用perfplot分析性能

我已经测试了每种方法在来自激光雷达的真实点云数据上的效率，这些数据以厘米为单位给出 3D 值（大约 2000 万个点，其中 100 万个是不同的点）。最快版本运行 2 秒。让我们看看 perfplot:

上的结果

import tensorflow as tf
import perfplot
import matplotlib.pyplot as plt
from time import time
path = tf.keras.utils.get_file('cubes.npz', 'https://github.com/loijord/lidar_home/raw/master/cubes.npz')
cubes = np.load(path)['array'].astype(np.int64) // 50
t = time()
fig = plt.figure(figsize=(15, 10))
plt.grid(True, which="both")
out = perfplot.bench(
        setup = lambda x: cubes[:x],
        kernels = [approach1a, approach1b, approach2, approach3, approach1a_on1D, approach1b_on1D, approach2_on1D],
        n_range = [2 ** k for k in range(22)],
        xlabel = 'cubes[:n]',
        title = 'Testing groupby max on cubes',
        show_progress = False,
        equality_check = False)
out.show()
print('Overall testing time:', time() -t)
# Overall testing time: 129.78826427459717

如何在 numpy 中有效地从 Nx3 点云阵列获取表面点？

How to get surface point from Nx3 point cloud array efficiently in numpy?

python

arrays

numpy

point-clouds

更新

使用perfplot分析性能