3D点空间网格化及网格编号标注的Python实现

现有存储3D点坐标的矩阵D，需将所有3D点量化到n×n×n规格的空间网格中，最终用每个点所属的网格编号表征原始点矩阵。

基础准备 #

首先导入numpy库，定义3D点矩阵：

import numpy as np
D = np.array([[-45.08341177,  34.40457052,   7.63253164],
       [-46.81391587,  34.35034554,   8.08080687],
       [-48.25988676,  34.29036266,   8.11271667],
       [-42.78217419,  24.51951629,   5.04130893],
       [-44.50295096,  25.34408096,   5.26748238],
       [-46.18559432,  25.60061165,   5.19218065]])

拆分三个维度的坐标数组：

x = D[:,0]
y = D[:,1]
z = D[:,2]

初始化网格相关参数，生成三个维度的网格边界：

r = np.zeros(x.shape)
n = 25 # 单维度网格数量，总网格数为n*n*n
xE = np.linspace(np.min(x), np.max(x), num = n+1 )
yE = np.linspace(np.min(y), np.max(y), num = n+1 )
zE = np.linspace(np.min(z), np.max(z), num = n+1 )

网格编号标注实现 #

原逻辑为按x→y→z的顺序遍历所有网格单元，对落在当前单元坐标区间内的点赋从1开始递增的编号，区间判断规则为左开右闭：即大于网格左边界、小于等于网格右边界。

直译版本（完全对齐原有逻辑） #

Python索引从0开始，对应调整循环范围即可直接复刻逻辑：

idx = 1
for ii in range(n):
    for jj in range(n):
        for kk in range(n):
            mask = (x > xE[ii]) & (x <= xE[ii+1]) & (y > yE[jj]) & (y <= yE[jj+1]) & (z > zE[kk]) & (z <= zE[kk+1])
            r[mask] = idx
            idx += 1

边界提示：x/y/z维度的最小值点不满足> 维度最小值（即第一个网格左边界）的判断条件，会保留初始值0，若需要覆盖最小值点，可将第一个网格的左边界判断改为>=。

高效向量化版本（无三重循环） #

点数量较大时三重循环运行效率极低，可通过numpy内置函数直接计算每个点的网格索引，速度提升明显：

# 计算每个点在三个维度上的网格位置（从0开始计数）
ix = np.searchsorted(xE, x, side='right') - 1
iy = np.searchsorted(yE, y, side='right') - 1
iz = np.searchsorted(zE, z, side='right') - 1
# 修正最小值点的索引越界问题
ix[ix < 0] = 0
iy[iy < 0] = 0
iz[iz < 0] = 0
# 按原有编号顺序计算全局网格编号（从1开始）
r = ix * n * n + iy * n + iz + 1

两个版本输出的网格编号结果完全一致，可根据数据规模选择使用。

内容的提问来源于stack exchange，提问作者The Peace of Programming