先把结论放在前面标题里写的 NML其实就是圈子里常说的 NLMNon-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 个像素搜索窗口边长 Spatch 边长 P那么总计算量大约是O(N × S² × P²)举个例子一张 256×256 的图搜索半径 7窗口 15×15patch 半径 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 遍历当前 patchq 遍历偏移后的 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偏移数量就是 (2r1)²。r7 时有 225 个偏移每张积分图是一张 (H1)×(W1) 的 float 矩阵。如果图像是 512×512单张积分图大约 2MB225 张就是 450MB 左右。跑到这个规模内存会开始紧张。实际工程里一般有两种处理方式。第一种是分块处理把偏移量分成几组每组算完权重后释放积分图再算下一组。第二种是降低积分图精度用 float32 替代 float64内存直接减半。如果只是学习原理小图上把全部积分图放在内存里完全没问题但心里要清楚这个方案的空间复杂度是 O((2r1)²·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, h10.0, patch_radius2, search_radius7): rows, cols src.shape pad search_radius patch_radius padded np.pad(src, pad, modereflect).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权重会失真。第二偏移量是相对于中心像素而言的所以候选像素坐标是 cidx、cjdy。第三自身位置的 patch 距离是 0权重是 exp(0)1也该参与归一化否则输出会偏暗。你拿小图跑一次就会发现这套写法正确性没问题但慢得没法接受。256×256 图上跑 search_radius7我当年第一次跑等到怀疑代码死循环了。所以在工程里需要一个更快的版本。3.3 积分图加速版完整实现积分图加速版的核心是提前对每个偏移量构建一张“平方差图”的积分图之后每个 patch 距离都是一次查表操作def build_integral(img): H, W img.shape sat np.zeros((H 1, W 1), dtypenp.float64) np.cumsum(np.cumsum(img, axis0), axis1, outsat[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, h10.0, patch_radius2, search_radius7): rows, cols src.shape pad search_radius patch_radius padded np.pad(src, pad, modereflect).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, dtypebool) 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 表示有效范围。例如 dx2 时目标区域是从第 0 行开始的 H−2 行源区域是从第 2 行开始的 H−2 行这样平移后每个有效位置的值就是原图像在该位置偏移之后的值。积分图加速版的主循环仍然是 Python 逐像素遍历但 patch 距离已经从 P² 次运算降到了 4 次数组索引和加减加速比大约等于 patch 面积。在我的环境里256×256、patch_radius2、search_radius7 的测试积分图版比朴素版快了二十倍以上肉眼从“完全等不下去”变成“能等十几秒”但离实时还差很远。3.4 两种实现性能对比与效果验证把三个版本放一起对比out_naive nlm_naive(noisy, h30, patch_radius2, search_radius7) out_integral nlm_integral(noisy, h30, patch_radius2, search_radius7) 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 版本的输出像素值应该高度一致积分图版略有一点点浮点精度差异属于正常现象。性能差异可以整理成一张参考表实现版本256×256 耗时参考计算复杂度关键点朴素版约 7 分钟量级每个 offset 都要做 P² 内层循环积分图版约 12 秒量级每个 offset 一次积分图查表OpenCV fastNlMeansDenoising约 0.1 秒量级C 内部多线程加精心内存布局这个耗时数字不是精确基准不同机器差异很大但数量级差距是真实的。它说明一个道理算法原理相同工程实现和底层优化对性能的影响可能比你想象得还要大。4. 工程落地中的常见坑与经验4.1 边界处理不能一刀切NLM 的边界处理比普通滤波更容易踩坑。原因是它有两种窗口比较相似度用的 patch 窗口和收集候选像素用的搜索窗口。如果只用零填充边缘区域会匹配到大量“黑色背景”的候选 patch这些 patch 和真实的边缘纹理完全不相似权重会畸变。我实测下来reflect 填充的效果最稳它能把边缘附近的图像结构“镜像翻转”出来让边缘像素在搜索范围内至少能找到和自己方向一致的候选块。但要注意这相当于人为引入了不存在的像素参与比较在图幅很小、边缘纹理很强的时候会出现边缘附近纹理被“拉花”的现象。更稳妥的做法是在积分图版本里只计算距离边界至少 search_radiuspatch_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 并不是所有场景的最优解。我常用的选型经验是场景建议算法原因实时视频预览OpenCV fastNlMeansDenoising底层优化成熟速度快强纹理图像BM3D 或 NLM 大 patch更好保留重复纹理高保真科研图像自研 NLM 并精细调参可控性强可复现低照度彩色图像在 Lab 空间对亮度通道做 NLM避免彩色通道互相污染NLM 最大的优点是原理简单、效果好作为图像去噪的入门算法它能让你理解“非局部自相似性”这个在 BM3D、深度学习去噪里仍然核心的底层思想。但要注意NLM 对高噪声σ30的效果会明显下降因为 patch 之间的相似度已经被噪声污染找出来的相似块可信度变差。5. 一些参数调优与后续扩展心得参数调优这件事光看公式和理论永远学不会。我自己在模拟项目 X 里测试时用过一组固定参数patch_radius2、search_radius7、h25在中等噪声图上效果稳定。但换了一张低照度图像后同样的参数把暗部噪点磨得干干净净却也把头发细节糊成一片。后来把 h 降到 15细节立刻回来了。我的经验是调 NLM 参数时一定要同时看两个东西一是整图 PSNR二是局部细节的放大对比图。PSNR 提升不代表视觉质量提升尤其对纹理区域过拟合测试指标的参数往往会牺牲主观质感。h 参数可以从大到小双向试探先在平坦区域表现满意的取值基础上再观察纹理区域是否出现过平滑若出现就逐渐减小 h直到纹理边缘开始保留。这个算法后续可以扩展的方向也很多。彩色图像可以把 RGB 三通道拼成一个向量来计算 patch 距离能避免逐通道独立去噪导致的颜色漂移。多尺度 NLM 可以配合图像金字塔先算粗糙层权重再映射回细尺度能缓解大噪声时的相似度判断失效。再往深走BM3D 本质上就是 NLM 加上了变换域协同滤波理解了 NLM 再看 BM3D 会顺畅很多。最后分享一个实际操作中很实用的技巧如果你在项目里只是需要“速度快一点的 NLM”可以考虑把搜索半径从 7 降到 4。测试下来大部分自然图像的相似 patch 集中在中心附近几个像素范围内半径从 7 降到 4视觉效果差距很小但计算量和内存占用直接少了一半以上。先把搜索范围收紧再慢慢调 h往往比一开始就上大搜索窗口更高效。