体素相关的快速计算

“体素”通常是指在三维空间中具有固定尺寸和位置的小立方体单元。

体素的优点包括：

• 易于处理和计算：在计算机图形学和三维建模中，体素的结构相对简单，计算和操作较为方便。
• 能精确表示物体的内部结构：对于一些需要了解物体内部信息的应用，如医学成像、地质建模等，体素可以提供更详细和准确的内部描述。
• 适合并行计算：因为体素之间相对独立，便于在多核或分布式计算环境中进行并行处理，提高计算效率。
• 数据结构简单：存储和管理体素数据的结构相对简单，降低了数据处理的复杂性。

本文主要介绍在使用体素过程中的一些快速计算的方法

1. 求重心

一般情况下并不需要求体素的重心，可以用中心或随机采样来代替，但是在一些特殊的应用场景，还是需要重心的，而这个重心却很难计算，此处介绍一种快速方法。

通过cumsum函数高效计算voxel中所有Point的X，Y，Z坐标之和/平均值，为了方便理解，下面对该trick的计算过程进行简单演示

# 假设我们有一组Point，以及该组Point中每个点对应的voxel的index
points = np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4], [5, 5, 5]])  # 5x3
indexes = np.array([0, 1, 2, 2, 3])
# 现在我们希望计算属于同一个voxel中的所有点的x,y,z坐标的平均值，暴力方法是通过for循环来实现
# 下面给出如何利用cumsum函数来非常简单的实现上述功能# step1: 对点云的第0维进行累加求和
points_cumsum = points.cumsum(0)
# point_cumsum = [[1, 1, 1], [3, 3, 3], [6, 6, 6], [10, 10, 10], [15, 15, 15]]# step2: 计算每个点与其前一个点是否落在同一个voxel中
kept = np.ones(points.shape[0], dtype=bool)
kept[:-1] = (indexes[1:] != indexes[:-1])
# kept: [True, True, False, True, True]# step3: 计算并输出同一个voxel中所有Point的xyz坐标之和
points_voxel = points_cumsum[kept]
points_voxel = np.concatenate((points_voxel[:1], points_voxel[1:] - points_voxel[:-1]))
# points_voxel: [[1, 1, 1], [2, 2, 2], [7, 7, 7], [5, 5, 5]]

2. 语义投票

对体素中每个点的语义求均值是错误的方法。随机采样并不鲁棒，投票是比较好的方法，但投票的开销很大。此处介绍一个快速计算方法。

2.1. 开销大的方法

当体素较多时，该方法速度极慢

min_values = np.min(points, axis=0)
points_voxel_coordinates = np.round((points - min_values) / voxel_size).astype(int)
sample_voxel_coordinates, sampled_indices = np.unique(points_voxel_coordinates, axis=0, return_inverse=True)
sampled_points = sample_voxel_coordinates * voxel_size + voxel_size * 0.5sampled_semantics = []
for i in range(len(sample_voxel_coordinates)):print(f"Processing voxel {i}/{len(sample_voxel_coordinates)}")voxel_point_indices = np.where(sampled_indices == i)[0]voxel_semantics = semantics[voxel_point_indices]voxel_semantic = scipy.stats.mode(voxel_semantics)[0][0]sampled_semantics.append(voxel_semantic)sampled_semantics = np.array(sampled_semantics)

2.2. 内存大的方法

该方法速度快，但当语义的类别较多时，其占用的内存较大。

num_points = points.shape[0]min_values = np.min(points, axis=0)
points_voxel_coordinates = np.round((points - min_values) / voxel_size).astype(int)
_, sampled_indices, counts = np.unique(points_voxel_coordinates, axis=0, return_inverse=True, return_counts=True)arg_indices = np.argsort(sampled_indices)
sampled_indices = sampled_indices[arg_indices]
points = points[arg_indices, :]
points = points.cumsum(0)kept = np.ones(num_points, dtype=bool)
kept[:-1] = (sampled_indices[1:] != sampled_indices[:-1])
counts_3d = np.tile(counts, (3, 1)).Tpoints = points[kept]
points = np.concatenate((points[:1], points[1:] - points[:-1]))
points = (points / counts_3d).astype(np.float32)semantics = pointcloud.point.semantic.numpy().flatten()
semantic_matrix = np.zeros((num_points, num_semantics), dtype=bool)
semantic_matrix[np.arange(num_points), semantics] = 1semantic_matrix = semantic_matrix[arg_indices, :]
semantic_cols = []for i in range(num_semantics):semantic_col = np.cumsum(semantic_matrix[:, i].astype(np.uint32))semantic_col = semantic_col[kept]semantic_col = np.concatenate((semantic_col[:1], semantic_col[1:] - semantic_col[:-1])).astype(np.uint8)semantic_cols.append(semantic_col)semantic_cols = np.stack(semantic_cols, axis=1)
semantic_matrix = np.argmax(semantic_cols, axis=1).astype(np.uint8)

2.3. 两者兼顾的方法

该方法虽然速度上比上一方法略慢，但解决了内存问题

num_points = points.shape[0]min_values = np.min(points, axis=0)
points_voxel_coordinates = np.round((points - min_values) / voxel_size).astype(int)
_, sampled_indices, counts = np.unique(points_voxel_coordinates, axis=0, return_inverse=True, return_counts=True)arg_indices = np.argsort(sampled_indices)
sampled_indices = sampled_indices[arg_indices]
points = points[arg_indices, :]
points = points.cumsum(0)kept = np.ones(num_points, dtype=bool)
kept[:-1] = (sampled_indices[1:] != sampled_indices[:-1])
counts_3d = np.tile(counts, (3, 1)).Tpoints = points[kept]
points = np.concatenate((points[:1], points[1:] - points[:-1]))
points = (points / counts_3d).astype(np.float32)semantics = pointcloud.point.semantic.numpy().flatten().astype(np.uint8)
semantics = semantics[arg_indices]num_kept = np.sum(kept)
max_values = np.zeros(num_kept, dtype=np.uint64)
max_indices = np.zeros(num_kept, dtype=np.uint8)for i in range(num_semantics):semantic_col = np.where(semantics == i, 1, 0).astype(np.uint64)semantic_col = np.cumsum(semantic_col)semantic_col = semantic_col[kept]semantic_col = np.concatenate((semantic_col[:1], semantic_col[1:] - semantic_col[:-1])).astype(np.uint8)sign = semantic_col > max_valuesmax_values[sign] = semantic_col[sign]max_indices[sign] = i