NumPy是否存在顺序判断的类where方法 实现点云高程图格网最大值更新

问题背景 #

我正在基于点云生成高程图，将点云投影到高程图时，可能存在多个点投影至同一格网单元的情况，需求是每个格网单元仅保留对应最高高度的点值。
当前涉及的变量定义如下：

• 高程图image：尺寸为(J,K)的float类型NumPy数组
• 待更新的格网单元坐标pix_points：尺寸为(2,N)的数组
• 对应点的高度值heights：尺寸为(1,N)的数组

初始实现的问题 #

我最初尝试的实现代码如下：

image[pix_points[1,:], pix_points[0,:]] = (
np.where(
    image[pix_points[1,:], pix_points[0,:]] <= heights,
    heights, 
    image[pix_points[1,:], pix_points[0,:]]))

该实现无法满足需求：np.where会先整体计算所有条件生成布尔映射，再基于映射选择heights或对应位置的原有值，所有对比都基于更新前的原始值计算，无法处理同一图像位置对应多个不同高度的场景。需求要求在同一次处理流程中，先写入某位置的高度值，后续遇到同一位置的其他点时，需要对比新高度与已写入的值再判断是否更新。
我不希望使用Python循环实现该逻辑，希望找到向量化的实现方案。

补充说明 #

对应操作的Python循环实现代码如下：

for p, h in zip(pix_points.T, heights):
    if image[p[1],p[0]] <= h:
        image[p[1],p[0]] = h

该循环的实际运行速度比预期更快，当N约为50万时运行时长约1.75s，点云规模较小时使用该循环版本完全可行。但由于该模块属于处理相机实时流点云的在线系统，需要进一步提升运行速度。我知道如果要获得更高性能可能需要使用Cython或C重写，但我对Cython编写及C的Python绑定开发不太熟悉，因此想先确认是否存在基于NumPy或其他Python库的实现方案。

解决方案 #

推荐方案：使用NumPy ufunc的at方法 #

这是纯NumPy生态下最匹配需求、性能最高的实现，无额外依赖，代码简洁：

np.maximum.at(image, (pix_points[1], pix_points[0]), heights.reshape(-1))

• 核心原理：NumPy常规花式索引为了提升性能会使用缓冲，同索引多次写入时只会保留最后一次操作的结果；而ufunc.at是无缓冲的原地操作，会按顺序对每个索引位置执行最大值计算，同位置多次命中时会基于已更新的值做比较，逻辑和手写Python循环完全一致。
• 性能表现：50万点规模下运行耗时通常在10~20ms，比纯Python循环快近百倍，完全满足实时点云处理的性能要求。
• 注意事项：使用前需要先过滤掉坐标超出[0,J)、[0,K)范围的点，避免触发索引越界错误；heights需要展平为和坐标数组长度一致的一维数组，避免形状不匹配报错。

备选方案：分箱聚合取最大值 #

由于格网最终只保留最大值，和点的处理顺序无关，也可以先对所有点按格网坐标分组求最大值，再一次性更新到高程图。点规模超过千万级时，可以用scipy.stats.binned_statistic_2d实现，内存占用更低：

from scipy.stats import binned_statistic_2d

# 按格网分箱统计每个格网内的点最大高度
res = binned_statistic_2d(
    pix_points[1], pix_points[0], heights.reshape(-1),
    statistic='max',
    bins=(J, K)
)
grid_max = res.statistic
# 仅更新有点落入的格网，和原有值取最大
valid = ~np.isnan(grid_max)
image[valid] = np.maximum(image[valid], grid_max[valid])

内容的提问来源于stack exchange，提问作者martinako