基于Numpy实现自定义向量化非线性滤波器的优化问询

在数字图像处理中，许多滤波器是非线性的，例如谐波均值滤波器（Harmonic Mean Filter）。已知Numpy提供大量向量化函数可大幅提升计算速度，但未发现适用于非线性掩码的相关函数，希望优化现有谐波均值滤波器实现，移除耗时的两层Python for循环，现有代码如下：

import math as m
import numpy as np
import cv2 as cv  # 补充原代码缺失的导入

def harmonic(im, ksize):
    # 复制原图像并转换格式
    result = im.copy().astype(np.float32)

    # 计算填充尺寸并对图像填充
    psize = m.floor(ksize/2) # 填充大小
    im = cv.copyMakeBorder(im, psize, psize, psize, psize, cv.BORDER_REFLECT)

    # 执行非线性操作
    for i in range(0, result.shape[0]):
        for j in range(0, result.shape[1]):
            # 获取与核尺寸一致的邻域
            neighbor = im[(i):(i+2*psize+1),(j):(j+2*psize+1)].astype(np.float32)
                        
            # 计算倒数和
            recp_sum = np.sum(np.reciprocal(neighbor,where= neighbor != 0).astype(np.float32))

            # 计算该邻域的谐波均值
            if (recp_sum != 0):
                result[i][j] = (float((ksize*ksize)/(recp_sum)))
                        
    return result.astype(np.uint8)

疑问列表 #

• 能否利用Numpy对数组执行任意自定义向量化操作？还是仅支持有限类型的操作？
• 若可以，具体该如何优化上述谐波均值滤波器代码？
• 关注到np.vectorize，但文档示例与问题关联度低，且不理解pyfunc参数：它是否真的能消除pyfunc中的Python循环？还是仅定义单个像素的映射规则？

解答 #

1. Numpy自定义向量化操作的可能性

Numpy本身不支持完全任意的自定义向量化操作，它的向量化优势来自底层C实现的内置函数（比如np.sum、np.reciprocal这类）。但可以通过两种方式实现自定义逻辑的高效计算：

• 利用Numpy的广播机制和现有函数组合，将自定义逻辑拆解为可向量化的步骤；
• 借助numba这类JIT编译工具，将自定义循环编译为高效机器码（不属于纯Numpy范畴）。
  纯Numpy环境下，只能基于现有向量化函数组合出所需逻辑，无法直接定义任意的向量化操作。

2. 谐波均值滤波器的Numpy优化方案

要移除两层Python循环，核心是用滑动窗口的向量化计算替代逐像素遍历。可以使用np.lib.stride_tricks.sliding_window_view（Numpy 1.20及以上版本支持）生成所有像素的邻域窗口，再批量计算：

优化后的代码：

import numpy as np
import cv2 as cv

def harmonic_vectorized(im, ksize):
    im_float = im.astype(np.float32)
    psize = ksize // 2
    # 对图像进行填充
    padded_im = cv.copyMakeBorder(im_float, psize, psize, psize, psize, cv.BORDER_REFLECT)
    
    # 生成所有滑动窗口，形状为 (图像高度, 图像宽度, 核高度, 核宽度)
    windows = np.lib.stride_tricks.sliding_window_view(padded_im, window_shape=(ksize, ksize))
    
    # 计算每个窗口的倒数和，忽略0值，axis=(2,3)对应窗口的两个维度
    recp_sum = np.sum(np.reciprocal(windows, where=windows != 0), axis=(2, 3))
    
    # 计算谐波均值，用np.where处理倒数和为0的情况
    result = np.where(recp_sum != 0, (ksize * ksize) / recp_sum, im_float)
    
    return result.astype(np.uint8)

说明：

• sliding_window_view通过内存视图生成所有邻域窗口，无Python循环，也不额外占用大量内存；
• 直接对整个窗口数组批量计算倒数和，利用Numpy的广播和轴求和实现高效计算；
• 用np.where替代逐像素条件判断，实现向量化的分支处理。

3. 关于np.vectorize的作用

np.vectorize本质是Python循环的封装，它并没有真正消除pyfunc中的Python循环，只是让函数可以接受数组作为输入，自动对每个元素调用pyfunc：

• pyfunc参数是你定义的、处理单个元素的函数；
• np.vectorize仅能让函数调用方式符合Numpy数组的使用习惯，性能提升非常有限，甚至和手动写Python循环差不多，不适合用来优化这类需要处理邻域的滤波器。

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