对大型numpy.memmap执行ndimage.convolve报56GiB内存分配失败如何解决

对大型numpy.memmap对象执行scipy.ndimage.convolve卷积操作时触发内存异常，具体报错信息：

异常类型：_ArrayMemoryError
异常详情：无法为形状为(3710, 1056, 3838)、数据类型为float32的数组分配56.0 GiB内存。
问题根因：ndimage.convolve运算过程中会创建常规内存驻留的numpy数组，全量数组大小超出系统可用内存上限。

• 手动分块卷积
  沿数组最长维度将大体积数据切分为内存可承载的重叠子块，逐块执行卷积后写回预创建的磁盘memmap结果文件。必须在块边缘预留与卷积核半径等长的重叠区域，卷积完成后裁掉重叠部分再写回，否则块拼接位置会出现明显边界伪影。单块大小建议控制在系统可用内存的70%以内，避免运算过程中临时变量占满内存。参考实现逻辑如下：

  import numpy as np
  from scipy import ndimage
  
  # 输入为只读磁盘映射数组，输出为提前创建的写模式磁盘映射数组，全程不占全量内存
  input_arr = np.memmap("input_vol.dat", dtype=np.float32, mode="r", shape=(3710, 1056, 3838))
  output_arr = np.memmap("conv_result.dat", dtype=np.float32, mode="w+", shape=input_arr.shape)
  
  kernel = np.ones((3,3,3), dtype=np.float32) / 27  # 替换为实际使用的卷积核
  pad_radius = [s//2 for s in kernel.shape]  # 各维度边缘补长
  tile_size = 256  # 根据实际可用内存调整，示例值下单块内存占用约3.7GiB
  
  # 沿最长的第三维度分块
  for z_start in range(0, input_arr.shape[2], tile_size):
      z_end = min(z_start + tile_size, input_arr.shape[2])
      # 取带重叠边的子块
      z_pad_start = max(0, z_start - pad_radius[2])
      z_pad_end = min(input_arr.shape[2], z_end + pad_radius[2])
      tile = input_arr[:, :, z_pad_start:z_pad_end]
      # 子块卷积
      tile_conv = ndimage.convolve(tile, kernel, mode="reflect")
      # 裁掉重叠边缘，写回输出数组对应位置
      write_start = z_start - z_pad_start
      write_end = write_start + (z_end - z_start)
      output_arr[:, :, z_start:z_end] = tile_conv[:, :, write_start:write_end]
      # 主动释放临时变量内存
      del tile, tile_conv
  

• 用支持外存惰性计算的框架自动分块
  不需要手动实现分块逻辑，使用dask.array加载磁盘映射数组，通过map_overlap接口封装卷积操作，框架会自动处理分块、边缘重叠、内存调度，计算结果可直接写回磁盘，全程不会将全量数组加载到内存。如果使用的卷积核是可分离类型，拆分为多个一维卷积逐维度执行，可将内存占用降低一个数量级。
• 优化运算参数减少内存开销
  若业务精度允许，将数组dtype从float32转换为float16可直接将内存需求减半；调用卷积函数时主动传入output参数指向预先创建好的memmap输出数组，避免函数内部自动创建全量大小的内存驻留数组。注意部分scipy版本对非内存驻留的output数组兼容性较差，正式使用前先用小尺寸样本测试验证。
• 替换为大体积影像专用处理库
  针对TB级三维体数据、遥感影像设计的处理库内置了分块IO、卷积优化逻辑，不需要手动实现分块，直接传入磁盘文件路径或memmap对象即可完成运算，内存控制表现远好于通用scipy接口。

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