NLM非局部平均去噪算法原理、积分图加速与Python实现

先把结论放在前面：标题里写的 NML，其实就是圈子里常说的 NLM（Non-Local Means），中文一般叫非局部平均算法。当年我第一次听到“非局部”三个字，以为是把整张图所有像素全部拿来做平均，后来真去读论文才发现，它指的是“在一张图里寻找和当前小块长得相似的其它小块”，而不是只盯着身边的邻居。这个思路对图像去噪来说挺反直觉，但效果确实好，尤其在纹理重复性强的区域。

这篇文章我打算直接从原理、复杂度、积分图加速、Python 代码、工程踩坑四个部分一次讲透。适合正在做图像预处理的同学，也适合面试前想快速搞懂经典去噪算法的朋友。代码我会给两个版本，一个是很容易理解的朴素实现，一个是能用在实际项目里的积分图加速版，两者对照着看才会真正理解什么叫“用空间换时间”。

1. 非局部平均：把“局部”两个字忘掉

1.1 为什么局部滤波总是差一口气

传统去噪思路都很“邻居党”：均值滤波把一个像素周围一圈像素取平均，高斯滤波是带距离权重地取平均，双边滤波更聪明一点，额外加上亮度差作为权重，保护边缘。这些方法本质是一样的，都假设“像素和它附近的像素有关系，远处没关系”。

问题恰恰出在这个假设上。自然图像里面大量存在远距离相似结构：一张砖墙照片里，横向平移几个像素的砖块纹理几乎一模一样；一片草地照片里，远处某根草周围的 patch 和当前像素周围的 patch 可能非常接近。对于这种图像，局部滤波只能拿到周围几个邻居的信息，一旦邻居都被噪声污染了，平均出来的结果自然也不干净。

双边滤波虽然通过亮度差权重保护了边缘，但它依然只在固定邻域内找相似像素。换句话说，它对“像素值相似的邻居”给予高权重，却对那些“距离远、结构相似”的像素完全视而不见。这就像你在车库找一把丢在客厅的钥匙，找得再认真，搜索范围错了，结果也不会对。

1.2 非局部自相似：在一张图里找“同类块”

非局部平均算法换了一个思路：不再以单个像素为单位，而是以 patch（小图像块）为单位进行比较。如果两个像素各自周围的 patch 非常相似，那么不管这两个像素隔得多远，都应该互相贡献权重。这就是“非局部”的含义——比较的单位从“局部窗口”提升到了“局部块”，但候选范围扩大到整个搜索窗口，不再局限于紧邻区域。

这个思路的来源是自然图像的自相似性。图像里重复出现的纹理、边缘、角点，形成了一种冗余结构。高斯噪声是逐像素独立的，它破坏不了这种大尺度结构间的相似关系。所以当我们把许多个“长相几乎一样”的 patch 中心像素加权平均，噪声会因为随机性被抵消，而真实结构因为相似性被保留下来。

用一句话概括：NLM 不是把一张图里的像素当作孤立的点去平均，而是假设图里存在大量“失散的双胞胎 patch”，把它们找出来互相印证。这也是它和传统局部算法的最大分水岭。

1.3 权重公式和参数直觉

NLM 的数学表达看起来并不复杂。对每个像素 i，输出值是搜索窗口内所有像素 j 的加权平均：

u(i) = Σj w(i,j) * I(j)

权重 w(i,j) 的定义是：

w(i,j) = exp( -||V(i) − V(j)||² / h² ) / Z(i)

其中 V(i) 表示以像素 i 为中心、大小为 patch_size 的邻域块，||V(i) − V(j)||² 是两个 patch 逐像素差值的平方和，也就是两个小块的“相似度距离”。h 是控制衰减速度的滤波强度参数，Z(i) 是归一化因子，保证所有权重加起来等于 1。

h 的物理意义很直观：噪声越小，patch 之间的距离越能反映真实结构差异，h 就不用太大；噪声越大，即使两个 patch 来自同一结构，像素差值也会被噪声拉大，这时候必须增大 h，让权重衰减得慢一些，才能容忍这种噪声带来的额外距离。经验公式一般取 h ≈ 10σ，σ 是噪声标准差，图像动态范围按 0–255 来算。

这里有个新手经常忽略的点：patch 大小会影响 h 的选择。patch 越大，距离值越大，h 也要相应调大。所以你从网上抄代码时，如果别人用的是 7×7 patch，你换成 5×5 patch，还拿同一套 h 参数跑，结果大概率纹理会糊掉。

2. 朴素NLM的复杂度与积分图加速原理

2.1 朴素实现慢在哪里

原理听懂之后，很多人第一反应是直接写三重循环：遍历每个像素，遍历搜索窗口内的每个候选像素，计算两个 patch 的平方差和。听上去很合理，但算一下复杂度就笑不出来了。

假设图像有 N 个像素，搜索窗口边长 S，patch 边长 P，那么总计算量大约是：

O(N × S² × P²)

举个例子：一张 256×256 的图，搜索半径 7（窗口 15×15），patch 半径 2（块大小 5×5），计算量大约是 65536 × 225 × 25，约 3.7 亿次像素级减法。这还只是一张小图。如果是 512×512，直接到 15 亿次量级，纯 Python 写出来跑一次够你去泡杯咖啡再等它慢慢算完。

朴素实现的瓶颈非常集中，patch 距离计算里的内层 P² 循环被反复执行。积分图加速的核心目标就是把这个 P² 干掉，让一次距离计算变成 O(1) 查表。

2.2 积分图加速的关键：平方差展开

积分图（Integral Image）大家应该不陌生，它能在 O(1) 时间内求出任意矩形区域的和。但要直接对“两个 patch 的差值平方和”做积分图，看起来不太可能，因为差值涉及两个不同位置的小块。

关键思路是把平方差展开：

||V(i) − V(j)||² = Σ (I(p) − I(q))² = Σ I(p)² + Σ I(q)² − 2 Σ I(p)·I(q)

这里 p 遍历当前 patch，q 遍历偏移后的 patch。拆开以后，第一项 Σ I(p)² 只依赖当前像素 i 周围 patch 的平方和，第二项 Σ I(q)² 只依赖候选像素 j 周围 patch 的平方和，这两项都各自只要一张积分图就能 O(1) 查出来。

难点在第三项 Σ I(p)·I(q)。它相当于当前 patch 和偏移 patch 的逐像素乘积和。仔细想会发现，对于一个固定的偏移量 (dx, dy)，q 的位置其实由 p 的位置偏移决定。也就是说，可以预先为每个偏移量构造一张“平方差图”：

B(i) = (I(i) − I(i + dx, i + dy))²

然后对 B 再做一次积分图。这样一来，任意位置的双 patch 距离，都能通过查 B 的积分图一次搞定。这是 NLM 积分图加速的核心思想：把“内层 patch 循环”转换成“预计算 + 查表”。

2.3 分块实现与内存取舍

每个偏移量对应一张积分图，搜索窗口越大，需要存储的积分图就越多。假设搜索半径 r，偏移数量就是 (2r+1)²。r=7 时有 225 个偏移，每张积分图是一张 (H+1)×(W+1) 的 float 矩阵。如果图像是 512×512，单张积分图大约 2MB，225 张就是 450MB 左右。跑到这个规模，内存会开始紧张。

实际工程里一般有两种处理方式。第一种是分块处理，把偏移量分成几组，每组算完权重后释放积分图，再算下一组。第二种是降低积分图精度，用 float32 替代 float64，内存直接减半。如果只是学习原理，小图上把全部积分图放在内存里完全没问题，但心里要清楚，这个方案的空间复杂度是 O((2r+1)²·H·W)，不是 O(H·W)。

我的建议是：先在小图上验证正确性，再考虑内存优化。不要一上来就上 512 图、搜索半径 15，那样很可能内存直接爆掉。

2.4 搜索窗口和 patch 半径怎么选

参数选择直接影响去噪效果和速度。patch 半径一般取 2 到 3，也就是 patch 大小 5×5 到 7×7。patch 太小，相似度判断容易受噪声干扰；patch 太大，又会让“相似块”变少，纹理细节容易被抹平。搜索半径一般取 7 到 15，搜索范围越大，找到真相似块的机会越多，但计算量和内存也越大。

h 参数的经验值前面说过，约等于 10σ。但不同图像内容要不同对待。平坦区域可以稍微加大 h，让平滑更彻底；纹理密集区域要减小 h，否则细节会被过度平均掉。我自己的习惯是先估计噪声水平，再用 10σ 作为基准，肉眼对比几张输出图，做 ±20% 微调。

3. Python代码实现与实操验证

3.1 环境、数据和评价指标

实现 NLM 只需要 Python 3、NumPy 和 OpenCV。OpenCV 这一步主要是用来做测试图像读取、加噪声和调用现成去噪库做对比。数据集不麻烦，直接用一张灰度测试图，再叠加上均值为 0、标准差为 σ 的高斯噪声。评价指标用 PSNR 比较直接。

测试代码可以这样准备：

import cv2 import numpy as np src = cv2.imread("test.png", cv2.IMREAD_GRAYSCALE).astype(np.float64) np.random.seed(42) sigma = 15.0 noisy = src + np.random.normal(0, sigma, src.shape)

先把图像转成 float64，后面计算平方差不容易溢出。噪声强度用 σ=15，在 0–255 动态范围里属于中等噪声，肉眼能看到明显颗粒感，但结构还在，很适合算法效果对比。

3.2 朴素版NLM完整实现

先看朴素版，逻辑最直白，方便和原理一一对应：

def nlm_naive(src, h=10.0, patch_radius=2, search_radius=7): rows, cols = src.shape pad = search_radius + patch_radius padded = np.pad(src, pad, mode="reflect").astype(np.float64) h2 = h * h dst = np.zeros_like(src) offsets = [(dx, dy) for dx in range(-search_radius, search_radius + 1) for dy in range(-search_radius, search_radius + 1)] for i in range(rows): for j in range(cols): ci, cj = i + pad, j + pad center_patch = padded[ci - patch_radius: ci + patch_radius + 1, cj - patch_radius: cj + patch_radius + 1] weight_sum = 0.0 value_sum = 0.0 for dx, dy in offsets: qi, qj = ci + dx, cj + dy patch_q = padded[qi - patch_radius: qi + patch_radius + 1, qj - patch_radius: qj + patch_radius + 1] if dx == 0 and dy == 0: dist2 = 0.0 else: dist2 = np.sum((center_patch - patch_q) ** 2) w = np.exp(-dist2 / h2) weight_sum += w value_sum += w * padded[qi, qj] dst[i, j] = value_sum / weight_sum return dst

这段代码有几个值得注意的点。第一，padding 用 reflect 模式而不是常数零填充，不然边缘像素会匹配到一堆全零 patch，权重会失真。第二，偏移量是相对于中心像素而言的，所以候选像素坐标是 ci+dx、cj+dy。第三，自身位置的 patch 距离是 0，权重是 exp(0)=1，也该参与归一化，否则输出会偏暗。

你拿小图跑一次就会发现，这套写法正确性没问题，但慢得没法接受。256×256 图上跑 search_radius=7，我当年第一次跑，等到怀疑代码死循环了。所以在工程里需要一个更快的版本。

3.3 积分图加速版完整实现

积分图加速版的核心是提前对每个偏移量构建一张“平方差图”的积分图，之后每个 patch 距离都是一次查表操作：

def build_integral(img): H, W = img.shape sat = np.zeros((H + 1, W + 1), dtype=np.float64) np.cumsum(np.cumsum(img, axis=0), axis=1, out=sat[1:, 1:]) return sat def query_sum(sat, x1, y1, x2, y2): return sat[x2, y2] - sat[x2, y1] - sat[x1, y2] + sat[x1, y1] def nlm_integral(src, h=10.0, patch_radius=2, search_radius=7): rows, cols = src.shape pad = search_radius + patch_radius padded = np.pad(src, pad, mode="reflect").astype(np.float64) H, W = padded.shape h2 = h * h sat_list = [] offset_list = [] for dx in range(-search_radius, search_radius + 1): for dy in range(-search_radius, search_radius + 1): shifted = np.zeros_like(padded) lo_x, hi_x = max(0, -dx), min(H, H - dx) lo_y, hi_y = max(0, -dy), min(W, W - dy) shifted[lo_x:hi_x, lo_y:hi_y] = padded[lo_x + dx:hi_x + dx, lo_y + dy:hi_y + dy] mask = np.zeros_like(padded, dtype=bool) mask[lo_x:hi_x, lo_y:hi_y] = True diff = np.where(mask, (padded - shifted) ** 2, 0.0) sat_list.append(build_integral(diff)) offset_list.append((dx, dy)) dst = np.zeros_like(src) for i in range(rows): for j in range(cols): ci, cj = i + pad, j + pad x1, y1 = ci - patch_radius, cj - patch_radius x2, y2 = ci + patch_radius + 1, cj + patch_radius + 1 weight_sum = 0.0 value_sum = 0.0 for (dx, dy), sat in zip(offset_list, sat_list): dist2 = query_sum(sat, x1, y1, x2, y2) w = np.exp(-dist2 / h2) weight_sum += w value_sum += w * padded[ci + dx, cj + dy] dst[i, j] = value_sum / weight_sum return dst

这里偏移构建的逻辑要注意：shifted 数组存的是 padded 平移后的结果，lo_x、hi_x 表示有效范围。例如 dx=2 时，目标区域是从第 0 行开始的 H−2 行，源区域是从第 2 行开始的 H−2 行，这样平移后每个有效位置的值就是原图像在该位置偏移之后的值。

积分图加速版的主循环仍然是 Python 逐像素遍历，但 patch 距离已经从 P² 次运算降到了 4 次数组索引和加减，加速比大约等于 patch 面积。在我的环境里，256×256、patch_radius=2、search_radius=7 的测试，积分图版比朴素版快了二十倍以上，肉眼从“完全等不下去”变成“能等十几秒”，但离实时还差很远。

3.4 两种实现性能对比与效果验证

把三个版本放一起对比：

out_naive = nlm_naive(noisy, h=30, patch_radius=2, search_radius=7) out_integral = nlm_integral(noisy, h=30, patch_radius=2, search_radius=7) psnr_noisy = cv2.PSNR(src, noisy) psnr_naive = cv2.PSNR(src, out_naive) psnr_integral = cv2.PSNR(src, out_integral)

h 取 30 的原因是 σ=15，按 10σ 估算。测试结果通常能看到 PSNR 从 24dB 左右提升到 30dB 上下，纹理区域保留得比高斯滤波、双边滤波都要干净。两个 NLM 版本的输出像素值应该高度一致，积分图版略有一点点浮点精度差异，属于正常现象。

性能差异可以整理成一张参考表：

这个耗时数字不是精确基准，不同机器差异很大，但数量级差距是真实的。它说明一个道理：算法原理相同，工程实现和底层优化对性能的影响，可能比你想象得还要大。

4. 工程落地中的常见坑与经验

4.1 边界处理不能一刀切

NLM 的边界处理比普通滤波更容易踩坑。原因是它有两种窗口：比较相似度用的 patch 窗口，和收集候选像素用的搜索窗口。如果只用零填充，边缘区域会匹配到大量“黑色背景”的候选 patch，这些 patch 和真实的边缘纹理完全不相似，权重会畸变。

我实测下来，reflect 填充的效果最稳，它能把边缘附近的图像结构“镜像翻转”出来，让边缘像素在搜索范围内至少能找到和自己方向一致的候选块。但要注意，这相当于人为引入了不存在的像素参与比较，在图幅很小、边缘纹理很强的时候，会出现边缘附近纹理被“拉花”的现象。

更稳妥的做法是在积分图版本里只计算距离边界至少 search_radius+patch_radius 的内部有效区域，边缘区域用更简单的局部滤波兜底。这样做会损失一圈边缘的高质量去噪，但至少不会出现明显伪影。工程上，我倾向于先跑完核心区域，再把边缘保留原图处理，因为边缘只有几个像素宽，对整体 PSNR 影响很小。

4.2 权重归一化和数值稳定性

权重归一化看着简单，但有个隐藏问题：当 h 比较小、patch 距离又普遍偏大时，所有权重都可能小到接近 0，这时候 value_sum 和 weight_sum 都会出现浮点下溢风险。尤其使用 float32 时，结果可能直接变成 NaN 或整片黑色。

解决方式有两层。第一层是尽量用 float64 做内部计算，输出前再转回原数据类型。第二层是减小 exp 里的负值幅度，可以在做归一化之前从所有 dist2 里减去最小距离，这样只会平移权重刻度，不会改变权重之间的比例关系。

代码上可以做一个小优化：

dist_list = [] for (dx, dy), sat in zip(offset_list, sat_list): dist_list.append(query_sum(sat, x1, y1, x2, y2)) min_dist = min(dist_list) for d, (dx, dy) in zip(dist_list, offset_list): w = np.exp(-(d - min_dist) / h2)

这招在强噪声、小 h 的情况下特别有用，能让权重整体保持在可表示的数级内，但不会让权重相对关系跑偏。

4.3 为什么我的积分图版没有快很多

经常有人跑完积分图版本后抱怨：怎么还是很慢？这个问题的根源在于 Python 最外层还挂着双重 for 循环。积分图只是把 patch 距离的内层循环去掉了，但每个像素仍然要遍历几百个偏移量，每次迭代都伴随 Python 对象开销和 zip 解包。图像越大，这部分开销越不可忽视。

想要进一步提速，最简单的方式是把主循环丢给 Numba JIT。代码结构基本不用改，把双循环包进 njit 函数里即可。另一个方向是用 NumPy 的滑动窗口和整块偏移批量计算，但实现复杂度明显上升。如果你只是要一个能用的去噪模块，我建议直接用 OpenCV 的现成接口，自研版本留着做学习和验证。

4.4 和其他去噪算法怎么选

NLM 并不是所有场景的最优解。我常用的选型经验是：

NLM 最大的优点是原理简单、效果好，作为图像去噪的入门算法，它能让你理解“非局部自相似性”这个在 BM3D、深度学习去噪里仍然核心的底层思想。但要注意，NLM 对高噪声（σ>30）的效果会明显下降，因为 patch 之间的相似度已经被噪声污染，找出来的相似块可信度变差。

5. 一些参数调优与后续扩展心得

参数调优这件事，光看公式和理论永远学不会。我自己在模拟项目 X 里测试时，用过一组固定参数：patch_radius=2、search_radius=7、h=25，在中等噪声图上效果稳定。但换了一张低照度图像后，同样的参数把暗部噪点磨得干干净净，却也把头发细节糊成一片。后来把 h 降到 15，细节立刻回来了。

我的经验是：调 NLM 参数时一定要同时看两个东西，一是整图 PSNR，二是局部细节的放大对比图。PSNR 提升不代表视觉质量提升，尤其对纹理区域，过拟合测试指标的参数往往会牺牲主观质感。h 参数可以从大到小双向试探，先在平坦区域表现满意的取值基础上，再观察纹理区域是否出现过平滑，若出现就逐渐减小 h，直到纹理边缘开始保留。

这个算法后续可以扩展的方向也很多。彩色图像可以把 RGB 三通道拼成一个向量来计算 patch 距离，能避免逐通道独立去噪导致的颜色漂移。多尺度 NLM 可以配合图像金字塔，先算粗糙层权重，再映射回细尺度，能缓解大噪声时的相似度判断失效。再往深走，BM3D 本质上就是 NLM 加上了变换域协同滤波，理解了 NLM 再看 BM3D 会顺畅很多。

最后分享一个实际操作中很实用的技巧：如果你在项目里只是需要“速度快一点的 NLM”，可以考虑把搜索半径从 7 降到 4。测试下来，大部分自然图像的相似 patch 集中在中心附近几个像素范围内，半径从 7 降到 4，视觉效果差距很小，但计算量和内存占用直接少了一半以上。先把搜索范围收紧，再慢慢调 h，往往比一开始就上大搜索窗口更高效。