高效计算点云的体素索引

Question

假设我有一个输入点云X，由一个维度数组N x 3表示。此数组中的每一行对应于XYZ空间中的一个点，介于-1和1之间。现在，让k成为定义体素网格分辨率的参数（例如k=3对应于维度3 x 3 x 3的体素网格）我正在寻找一种有效的方法来计算每个点对应体素的索引。这或多或少是我目前使用NumPy的方式（为了清晰起见，写得更生动）：

    # Generate some random input point cloud, 100 points
    X = np.random.randn(100, 3)
    # Define the resolution of the grid (e.g. say 3 x 3 x 3)
    k = 3
    # Iterate through points of input point cloud
    for j, point in enumerate(X):
        # Placeholder for voxel membership
        partitions = np.zeros(3,)

        for i, p in enumerate(point):
            for d in range(k):
                # Points are between -1 and 1, so the interval for each dimension is [-1, 1]
                # Length 2, "left"/"lower" end is -1
                if p <= (-1 + (d + 1) * 2 / k):
                    partitions[i] = d

        # Compute the index of the voxel this point belongs to
        # Can think of this as the base 10 representation of the base k number given by (x, y, z) in partitions
        # Voxels are indexed such that (0, 0, 0) --> index 0, (0, 0, 1) --> index 1, (0, 0, 2) --> 
        # index 2, (0, 1, 0) --> index 3, etc.
        p_reversed = np.flip(partitions)
        idx= 0
        for d in range(3):
            idx += (k ** d) * p_reversed[d]

       # Now idx corresponds to the index of the voxel to which point j belongs

随着N和k的增加，这一比例明显降低；是否有更有效的实施方案

Answer 1

以下是我的方法：

N, k = 3, 10
np.random.seed(1)
X = np.random.normal(-1,1, (100,3))

# X:
# [[ 0.62434536 -1.61175641 -1.52817175]
# [-2.07296862 -0.13459237 -3.3015387 ]
# [ 0.74481176 -1.7612069  -0.6809609 ]]

# output:
np.argmax(X[:,:,None] < np.linspace(-1,1,k)[None, None, :], axis=-1)

输出：

array([[8, 0, 0],
       [0, 4, 0],
       [8, 0, 2]])

Answer 2

实际上，您正在将点值作为浮点值与一系列介于-1和1之间的其他浮点值进行比较

但是，您要做的是计算（一次）一个生成值的函数。执行简单的计算，而不是迭代

理想情况下，这个简单的函数应该是numpy可以分布在点值的各个列上的函数

更理想的情况是，它可以分布在整个阵列中，允许您在单个操作或一系列操作中使用二维阵列

因为您使用的是固定的体素大小，并且因为您对所有维度使用相同的大小和范围，我认为您可以通过简单的减法和乘法来实现：

1. 首先，减去将量程起点移动到零所需的量。在这种情况下，由于您的范围是[-1,1]，您可以将点值减去-1（或加+1），使其从0开始。

2. 接下来，将点值“缩放”到范围[0，1]中。您可以乘以范围长度的倒数（高-低==1--1==2），因此乘以1/2==0.5。

3. 此时，中间值是该点出现的范围的分数。因此，通过将该分数乘以体素范围的大小（3），将它们映射到体素空间。并将结果值显式（通过函数）或隐式（通过数据类型）转换为整数。Zen of Python说显式比隐式好，所以这是我的偏好。

如何编写此代码

RANGE_START = -1
RANGE_END = 1
VOXEL_K = 3

# Make up some points in the valid range
CLOUD = np.random.random((100, 3)) * (RANGE_END - RANGE_START) + RANGE_START

# shift points to 0-based:
zero_based_points = CLOUD - RANGE_START

# convert points to [0, 1) fraction of range
fractional_points = zero_based_points / (RANGE_END - RANGE_START)

# project points into voxel space: [0, k)
voxelspace_points = fractional_points * VOXEL_K

# convert voxel space to voxel indices (truncate decimals: 0.1 -> 0)
voxel_indices = voxelspace_points.astype(int)

注意：浮点值可以是nan或inf，它们不能很好地转换为整数。因此，您可能应该预处理您的点，以某种方式过滤这些值（将它们替换为sentinel值或从数据集中删除它们或…？）