python - 如何在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 和 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]

请注意，它优于先前的方法，这是一个可以进一步优化的信号。
方法三
部分 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.]]) , 除了 npi 的情况导致 np.array([[0. 0. 2.], [1. 2. 1.], [1. 1. 3.]]) 的方法
更新
如果您对扁平数组执行相同的算法(这是降维的结果)，则可以进一步优化它。 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 值(大约 20M 点，其中 1M 是不同的点)。最快的版本在 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