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

方法 2

如果您想要更简单的方法，您还可以使用写在numpy 之上的numpy_indexed，并允许用更少的脚本解决分组问题：

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.]])，除了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 值给出（大约 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