基于Scipy.ndimage对LiDAR高程数组迭代处理的问题

初始数据 #

初始LiDAR高程数组：

import numpy as np

Existing_example_arrayV0 = np.array([
    [0, 0, 1, 0, 0, 0, 0],
    [0, 1, 1, 1, 0, 0, 0],
    [0, 1, 1, 1, 1, 0, 0],
    [1, 1, 1, 1, 1, 0, 0],
    [0, 1, 1, 1, 1, 0, 0],
    [0, 1, 1, 0, 0, 0, 0]
])

第一次迭代实现 #

使用scipy.ndimage.binary_erosion对数组中值≥1的区域做水平3元素结构的腐蚀，为腐蚀区域元素加100，得到V1：

from scipy.ndimage import binary_erosion

# 生成原始有效区域掩码
original_mask = Existing_example_arrayV0 >= 1
# 定义腐蚀结构元素（水平3个相邻像素）
structure = [[1,1,1]]

# 第一次迭代
Existing_example_arrayV1 = Existing_example_arrayV0.copy()
erosion_1 = binary_erosion(original_mask, structure=structure, iterations=1)
Existing_example_arrayV1[erosion_1] += 100

得到结果：

Existing_example_arrayV1 = np.array([
    [0, 0, 1, 0, 0, 0, 0],
    [0, 1, 100, 1, 0, 0, 0],
    [0, 1, 100, 100, 1, 0, 0],
    [1, 100, 100, 100, 1, 0, 0],
    [0, 1, 100, 100, 1, 0, 0],
    [0, 1, 1, 0, 0, 0, 0]
])

目标：第二次迭代结果 #

需要得到V2，仅对第一次腐蚀区域再次腐蚀后的中心像素再加100：

Existing_example_arrayV2 = np.array([
    [0, 0, 1, 0, 0, 0, 0],
    [0, 1, 100, 1, 0, 0, 0],
    [0, 1, 100, 100, 1, 0, 0],
    [1, 100, 200, 100, 1, 0, 0],
    [0, 1, 100, 100, 1, 0, 0],
    [0, 1, 1, 0, 0, 0, 0]
])

错误尝试分析 #

你之前的代码存在三个问题：

1. 导入错误：import scipy.ndimage as binary_erosion 将整个模块重命名为函数名，正确方式是直接导入binary_erosion函数
2. 未修改数组：仅调用binary_erosion但未将结果用于修改数组元素
3. 循环逻辑错误：未基于上一次腐蚀结果迭代，且循环范围仅执行1次却未赋值

解决方案 #

以下是支持任意次数迭代的实现，核心逻辑是每次在上一次腐蚀的结果基础上再次腐蚀，然后为新的腐蚀区域元素加100：

import numpy as np
from scipy.ndimage import binary_erosion

def multi_iteration_elevation(initial_arr, num_iterations):
    current_arr = initial_arr.copy()
    original_mask = current_arr >= 1
    structure = [[1,1,1]]  # 水平3元素腐蚀结构
    prev_erosion = None
    
    for i in range(1, num_iterations + 1):
        if i == 1:
            # 第一次腐蚀基于原始掩码
            current_erosion = binary_erosion(original_mask, structure=structure, iterations=1)
        else:
            # 后续迭代基于上一次的腐蚀结果
            current_erosion = binary_erosion(prev_erosion, structure=structure, iterations=1)
        # 为当前腐蚀区域加100
        current_arr[current_erosion] += 100
        prev_erosion = current_erosion
    
    return current_arr

# 示例：执行2次迭代得到V2
Existing_example_arrayV2 = multi_iteration_elevation(Existing_example_arrayV0, 2)
print(Existing_example_arrayV2)

# 执行N次迭代只需修改num_iterations参数
# 比如5次迭代：
# Existing_example_arrayV5 = multi_iteration_elevation(Existing_example_arrayV0, 5)

逻辑说明 #

• 每次迭代的腐蚀范围会逐步缩小（腐蚀操作的特性），因此每次加100的区域是上一次腐蚀区域的核心部分
• 函数通过prev_erosion保存上一次的腐蚀结果，避免重复计算原始掩码的多次腐蚀，提升效率
• 支持任意次数的迭代，只需修改num_iterations参数即可

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