非局部均值滤波(Non-Local Means)算法解析：从2005年原始论文到Python实现

1次阅读

没有评论

共计 2157 个字符，预计需要花费 6 分钟才能阅读完成。

在图像处理中，噪声是常见的问题。传统线性滤波器如高斯模糊通过加权平均来平滑噪声，但这种方法有个致命缺点——它会同时模糊图像的边缘和细节。中值滤波在保留边缘方面表现稍好，但对于复杂纹理区域仍然会造成信息丢失。

2005 年 Buades 等人提出的非局部均值 (NLM) 算法突破了这一局限。其核心思想是：图像中可能存在多个相似区域，可以利用这些非局部信息来进行去噪。

权重的计算公式为：

$$w(p,q)=e^{-\frac{|v(p)-v(q)|^2}{h^2}}$$

其中：
– $p$ 是当前像素
– $q$ 是搜索窗口内的参考像素
– $v(p)$ 和 $v(q)$ 是像素邻域的特征向量
– $h$ 是控制衰减的参数

以下是完整的 Python 实现，使用 OpenCV 和 NumPy 进行优化：

import numpy as np
import cv2

def nlm_denoise(img, h=10, search_window=21, patch_size=7):
    """
    非局部均值去噪实现
    :param img: 输入图像(灰度)
    :param h: 滤波参数，控制权重衰减
    :param search_window: 搜索窗口大小
    :param patch_size: 相似块大小
    """
    pad = patch_size // 2
    img_pad = cv2.copyMakeBorder(img, pad, pad, pad, pad, cv2.BORDER_REFLECT)

    # NOTE: 使用积分图加速距离计算
    integral_sq = cv2.integral(np.square(img_pad).astype(np.float32))
    integral = cv2.integral(img_pad.astype(np.float32))

    H, W = img.shape
    output = np.zeros_like(img)

    for i in range(H):
        for j in range(W):
            # 当前像素在 padding 后图像中的位置
            p_i, p_j = i + pad, j + pad

            # 搜索窗口范围
            i_min = max(p_i - search_window//2, pad)
            i_max = min(p_i + search_window//2, H + pad - 1)
            j_min = max(p_j - search_window//2, pad)
            j_max = min(p_j + search_window//2, W + pad - 1)

            # 计算当前块特征
            v_p = integral[p_i+pad:p_i+pad+1, p_j+pad:p_j+pad+1] - \
                  integral[p_i-pad:p_i-pad+1, p_j+pad:p_j+pad+1] - \
                  integral[p_i+pad:p_i+pad+1, p_j-pad:p_j-pad+1] + \
                  integral[p_i-pad:p_i-pad+1, p_j-pad:p_j-pad+1]

            weights = []
            values = []

            # 遍历搜索窗口
            for m in range(i_min, i_max+1):
                for n in range(j_min, j_max+1):
                    # 计算块距离
                    v_q = integral[m+pad:m+pad+1, n+pad:n+pad+1] - \
                          integral[m-pad:m-pad+1, n+pad:n+pad+1] - \
                          integral[m+pad:m+pad+1, n-pad:n-pad+1] + \
                          integral[m-pad:m-pad+1, n-pad:n-pad+1]

                    dist = np.sum((v_p - v_q)**2) / (patch_size**2)
                    weight = np.exp(-dist/(h**2))

                    weights.append(weight)
                    values.append(img_pad[m, n])

            # 归一化权重
            weights = np.array(weights)
            weights = weights / np.sum(weights)

            # 加权平均
            output[i, j] = np.sum(weights * np.array(values))

    return output

if __name__ == "__main__":
    # 测试用例
    img = cv2.imread("noisy_image.png", 0)

    # 参数设置
    h = 15  # 经验值，噪声越大 h 应该越大
    search_window = 21  # 搜索范围
    patch_size = 5  # 相似块大小

    # 执行去噪
    denoised = nlm_denoise(img, h, search_window, patch_size)

    # 保存结果
    cv2.imwrite("denoised.png", denoised)