简介基于MATLAB的非局部均值NLM滤波去噪源码面向图像处理初学者、算法研究人员及需要快速去除高斯噪声的开发者。NLM算法利用图像自相似性通过比较像素块差异并加权平均来还原干净像素对高斯噪声抑制效果良好尤其适用于生物医学图像、自然图像等纹理重复区域较多的去噪场景。资源为单个RAR压缩包仅含1个.m源文件包体大小约910B代码精炼可直接在MATLAB中运行使用标准图像读写函数即可测试无需额外工具箱。程序完整实现了图像预处理、搜索窗口与比较窗口参数设置、块相似度计算、高斯核加权平均及结果输出等核心流程且源码设计模块化参数集中便于修改适合在此基础上扩展其他滤波算法。已有1426人学习下载适合作为数字图像处理课程的实验材料或算法入门参考通过调整窗口尺寸与噪声强度可定量研究去噪质量与边缘保持性能并可与中值滤波、维纳滤波等经典算法对比也可作为向C、Python等语言迁移的原型基础。1. 非局部均值滤波为什么能替代高斯模糊做图像去噪一张暗光环境下拍的夜景ISO推到3200放大看全是颗粒边缘也发虚。常见做法是先上高斯模糊噪声是压住了但文字边缘、头发丝也被一并抹掉图像像蒙了一层纱。非局部均值Non-local Means, NLM走的是另一条路不只看邻域而是到整张图里找与当前像素块相似的块用这些块的均值来还原中心像素。2005年Buades等人提出这个思路后它就成了图像去噪里“用结构换噪声抑制”的经典方案后来很多深度学习降噪网络也拿它做对比基线。这个MATLAB工程里的NLmeansfilt.m就是一个从零实现的NLM滤波器适合做课程设计、算法实验以及作为调参和对比的起点。接下来我会从算法原理一步步拆到能跑的代码再讲清楚参数怎么调。2. 块相似度与权重计算NLM怎么知道哪些块该参与平均2.1 从局部均值到非局部均值我们通常接触的均值滤波、高斯滤波都是用固定窗口对中心像素做加权平均。窗口内的像素因为空间上相邻默认灰度相近所以可以互相平均。这种做法的前提是窗口内像素属于同一个“物体表面”。一旦窗口跨过边缘背景暗像素和前景亮像素被平均在一起边缘自然就糊了。NLM改变的是“参与平均的对象来源”。它认为自然图像中有大量冗余纹理比如一面砖墙、一片草地、一段天空不同位置的小块结构会重复出现。所以给当前像素找相似块时不必局限在它周围几像素而是扩大到整个搜索范围。越相似的块权重越高平均出来的值才能保留原有结构。具体来说对于每个像素 x要估计干净值 f(x)计算公式为f_hat(x) Σ w(x,y) * g(y) / Σ w(x,y)其中 y 是搜索窗口内所有可能的位置g(y) 是带噪图像在 y 处的灰度w(x,y) 是相似度权重。分母是所有权重之和目的是归一化保证输出亮度不变。2.2 块相似度为什么直接比较单个像素不够如果只比较 x 和 y 两个单点灰度很难区分“恰好都被噪声抬高了”和“结构本身相同”。NLM用块对块比较。以 x 为中心取一个大小为 P×P 的块 N_x以 y 为中心取同样大小的块 N_y然后计算两个块的欧氏距离平方d(x,y) ||N_x - N_y||^2d 越小说明两个局部结构越相似。注意带噪图像上的 d 本身也包含噪声干扰所以原始NLM算法在计算距离前会给块内每个像素乘一个高斯核权重。这样中心像素对距离的贡献更大距离结果对高频噪声也不那么敏感。这一点和单纯把两个窗口做差再求平方和不一样高斯核相当于先对距离度量做了一次平滑。权重部分NLM采用指数核映射w(x,y) exp(-d(x,y) / h^2)其中 h 是人为设定的滤波强度参数。用指数形式的原因很直接距离小的时候权重接近1距离大的时候权重指数衰减到接近0比线性阈值更容易区分“相似”和“不相似”。h 选多大决定了算法是更接近均值滤波还是更接近原图保留。工程里通常把 h 和图像噪声标准差 σ 绑定推荐 h c·σc 在0.5到1.0之间第4章会给出一组可用的对照表。2.3 搜索窗口、比较窗口与计算量实现时很少让 x 去匹配全图每一个像素而是限定在一个 S×S 的搜索窗口内。搜索窗口越大能找到的相似块越多去噪效果越强但计算量按 S² 增长。比较窗口 P 也不宜过大P 越大结构判别的语义越宏观同时更难找到完全一致的块。经典论文里常用 S21、P7或者 S11、P5具体效果和图像分辨率有关。我的理解是搜索窗口决定“候选数量”比较窗口决定“筛选严格程度”。这两个概念容易混列个表对比一下对比维度局部均值滤波NLM非局部均值滤波候选来源当前像素周围固定窗口搜索窗口内任意像素相似判定空间距离近即认为相似块间欧氏距离小才相似边缘表现边缘被平均掉近似结构参与平均边缘保留主要参数窗口大小 W搜索窗口 S、比较窗口 P、强度 h这个表解释了为什么NLM在噪声标准差明显时能比高斯模糊少糊边缘。本质上它把一个“空间邻域平均”问题改写成了“结构聚类平均”问题而结构相似度的度量决定了聚类是否准确。在MATLAB里计算两个patch的欧氏距离最直接的方式是% 生成两个patch计算欧氏距离未加高斯权重的版本 diff patchX - patchY; dist sum(diff(:).^2); % 加高斯核后的版本用逐点平方乘高斯核 diffSq diff.^2; distWeighted sum(diffSq .* gaussKernel, all);distWeighted 就是第3章代码里反复用的距离项。gaussKernel 的尺寸与 patch 相同通常用meshgrid生成sigma取1到2。3. NLmeansfilt.m 实现从函数签名到主循环3.1 函数签名与参数设计拿到这个MATLAB工程后核心文件只有一个 NLmeansfilt.m。它把NLM算法封装成一个普通MATLAB函数输入图像和三个参数输出去噪结果。我建议的函数签名如下function denoised NLmeansfilt(img, searchSize, patchSize, h)四个参数的含义如下表参数含义常见取值范围img输入灰度图像支持double或uint8彩色图需逐通道处理searchSize搜索窗口边长奇数11、21、31patchSize比较窗口边长奇数3、5、7h滤波强度与噪声标准差σ相关0.4σ ~ 1.2σ函数内部顺序我按“预处理 - 距离计算 - 权重累加 - 归一化”来组织。编译成常规函数而不是脚本主要为了后面循环调参方便。3.2 预处理与padarray边界处理边界问题是实现里最容易出错的地方。搜索窗口会越过图像边界所以第一步是填充图像。推荐使用replicate填充也就是复制边缘像素它不会像zero padding那样引入黑色人工边界对边缘像素的还原也更真实。function denoised NLmeansfilt(img, searchSize, patchSize, h) if nargin 4 error(请输入图像、搜索窗口大小、比较窗口大小和滤波强度h); end if ~ismatrix(img) error(当前实现仅支持灰度图彩色图请逐通道处理); end img im2double(img); [rows, cols] size(img); searchHalf floor(searchSize / 2); patchHalf floor(patchSize / 2); pad searchHalf patchHalf; padded padarray(img, [pad, pad], replicate); denoised zeros(rows, cols); accumulator zeros(rows, cols); weightSum zeros(rows, cols); % 高斯核对block内像素加权中心像素权重最大 [meshX, meshY] meshgrid(-patchHalf:patchHalf, -patchHalf:patchHalf); sigmaKernel 1; gaussKernel exp(-(meshX.^2 meshY.^2) / (2 * sigmaKernel^2)); gaussKernel gaussKernel / sum(gaussKernel(:));这里把输入转成double避免uint8计算溢出。im2double之后像素范围变成0~1噪声方差也要对应转换。gaussKernel是经典论文里用来加权欧氏距离的高斯核sigma取1是常见设定如果你想让块中心的影响更大可以调小到0.8但视觉差异不大。3.3 主循环搜索窗口内的权重累加接下来是核心循环。图像数据量大完全避免循环不现实但可以把“循环搜索偏移量”和“向量化计算整幅图像距离图”结合起来速度通常比逐像素三重循环快很多。for offsetY -searchHalf : searchHalf for offsetX -searchHalf : searchHalf % 当前偏移下计算整幅图上每个像素对应的块距离 diffSum zeros(rows, cols); for pY -patchHalf : patchHalf for pX -patchHalf : patchHalf baseY (1:rows) pad pY; baseX (1:cols) pad pX; shiftY baseY offsetY; shiftX baseX offsetX; diff padded(baseY, baseX) - padded(shiftY, shiftX); diffSum diffSum diff.^2 * gaussKernel(pY patchHalf 1, pX patchHalf 1); end end % 指数核映射得到权重 weight exp(-diffSum / (h^2)); % 取偏移后的中心像素值做加权平均 centerShiftY (1:rows) pad offsetY; centerShiftX (1:cols) pad offsetX; accumulator accumulator weight .* padded(centerShiftY, centerShiftX); weightSum weightSum weight; end end denoised accumulator ./ weightSum; end这段代码的流程是逐个扫描搜索窗口内的相对位置(offsetY, offsetX)。对每个位置一次性计算整幅图像上所有像素的“当前块与偏移块”的差平方并用高斯核加权得到diffSum。再用指数映射得到权重。accumulator是权重乘以像素值的累加weightSum是权重累加。最后相除就是归一化加权平均。参数说明offsetX和offsetY的取值范围由searchHalf决定代表搜索窗口。例如searchSize21则searchHalf10外层循环共441个偏移位置。每个偏移内部还有patchSize²个内层循环。总循环次数约为441 * patchSize²这就是NLM慢的主要原因也是第5章用积分图加速的出发点。注意这段代码没有跳过中心点自身。当offset都为0时diffSum为0权重exp(0)1自身会参与加权平均符合原始NLM算法的设计。如果你发现输出的噪声残留偏多可以尝试把自身权重去掉但通常保留自身效果更稳。3.4 权重归一化与数值保护最后一步经常被忽略如果图像某些区域的所有patch都和其他patch离得很远权重全部接近0weightSum可能出现NaN。一种稳健做法是给weightSum加保护weightSum(weightSum eps) 1; denoised accumulator ./ weightSum;这一步对纯色大面积区域特别重要。比如天空区域纹理少距离都很小权重不算太小但如果图像本身是黑色背景加少量亮点黑色区域的权重就可能全部为0。加上eps保护后输出至少不会出现NaN白点。4. 参数怎么设搜索窗口、patch大小与噪声强度的调优4.1 一组可以直接上手的参数对照参数选择是NLM被问得最多的部分。我结合工程实践和论文里的推荐值整理了一份对照表。噪声强度用加性高斯白噪声的均方差σ表示注意这里的σ是图像范围在0~1之后的绝对值。噪声水平searchSizepatchSizeh 推荐范围σ0.05轻度1130.5σ ~ 0.6σσ0.1中度2150.6σ ~ 0.8σσ0.2重度3170.8σ ~ 1.0σσ0.3极重度4191.0σ ~ 1.2σ这张表只是起点。h太小权重分布陡峭几乎没有像素参与去噪h太大权重分布趋于平坦NLM退化成均值滤波。我实际调参时会固定patchSize再用PSNR曲线扫描h准备一张基准图加已知噪声σ。固定searchSize遍历h [0.2σ, 0.4σ, …, 1.4σ]。对每个h计算PSNR取峰值对应的h再微调searchSize和patchSize。4.2 噪声方差估计与PSNR验证脚本做算法实验时必须加已知噪声否则没法量化效果。下面是生成带噪图像并调用NLmeansfilt的完整测试脚本% 读入灰度图像并加高斯噪声 img im2double(imread(cameraman.png)); noiseSigma 25 / 255; noisy img noiseSigma * randn(size(img)); % 调用NLMh设为0.7倍噪声标准差 denoised NLmeansfilt(noisy, 21, 5, 0.7 * noiseSigma); % 计算PSNR mse mean((img(:) - denoised(:)).^2); psnrVal 10 * log10(1 / mse); fprintf(PSNR: %.2f dB\n, psnrVal); % 显示结果 subplot(1,3,1); imshow(img); title(原图); subplot(1,3,2); imshow(noisy); title(含噪图); subplot(1,3,3); imshow(im2uint8(denoised)); title(NLM去噪);这段脚本里randn生成标准正态分布噪声乘上noiseSigma再叠加到图上就得到指定方差的高斯噪声。PSNR计算时double图像最大像素值取1所以用1除以MSE。im2uint8是为了让显示结果和imwrite输出的位深正常因为imshow对double图像按0~1范围映射超出范围会显示成全白或全黑。如果真实环境里不知道噪声强度可以用Laplacian高通滤波估计噪声方差对含噪图像做卷积取结果的标准差再乘一个修正系数。这种估算精度够用不必上贝叶斯估计。4.3 一个容易被忽视的坑uint8溢出与h的单位很多人下载源代码后直接用uint8图像跑噪声方差也用整数25而不是25/255。如果函数内部直接计算diff.^2uint8的差值为负时会被截断成0整个距离计算就全错了。所以第3章代码里特别强调im2doubleh也必须和图像范围一致。如果你看到去噪结果出现大块均匀平板十有八九是h的单位不一致。另一个坑是搜索窗口过大图像里唯一的纹理结构会被过度平均。比如拍一张黑色桌面的图搜索窗口覆盖整个桌面不同位置块差异极小权重接近均匀NLM就退化为均值滤波。这种情况下即使算法没问题效果也体现不出来。这也是为什么做对比实验时常用lena、cameraman这类纹理丰富的图像而不会用纯色块拼接图。5. 用积分图把运行时间降一个量级5.1 为什么NLM这么慢第3章的嵌套循环在MATLAB里跑512×512图像searchSize21时通常要几十秒调参时非常耽误时间。慢的根源在于每个offset都要计算整幅图的patch差平方再叠加。一个常见优化是积分图加速一组patch的差平方距离可以拆成当前块平方项、偏移块平方项和交叉项这三项都能用积分图在常数时间内查询把patchSize²的逐点循环变成O(1)查表。5.2 简化版积分图加速demo下面是一段在结构上可运行的积分图核心逻辑只计算一个固定offset的距离图offsetY 2; offsetX -3; % 三项都能用图像平移得到 A padded.^2; B circshift(padded, [offsetY, offsetX]).^2; C padded .* circshift(padded, [offsetY, offsetX]); % 对三者分别做二维累积和得到积分图 intA cumsum(cumsum(A, 1), 2); intB cumsum(cumsum(B, 1), 2); intC cumsum(cumsum(C, 1), 2); % 窗口和可以用积分图O(1)取出 dist (intA - 2 * intC intB) / patchArea;关键点是对于固定offset计算整幅图像的A、B、C三个矩阵后任何patch大小的窗口和都能通过一次加减得到完全不用内层循环遍历patch。缺点是每个offset要保存三个积分图搜索窗口大了以后内存占用明显所以实际工程里会限制搜索范围或者把图像分块处理。5.3 验证你的NLmeansfilt实现是否写对写完NLM后一定要做两个极端自检。把h设成1e-8输入不含噪的原图输出应该几乎等于原图因为任何带噪声的相似块都会被指数映射压成0权重只有完全相同的块能留下。再把h设成10输出应该接近整幅图像均值因为所有权重被压平成接近1退化成全局平均。这两个测试能快速暴露距离计算和归一化的问题。之后再用已知噪声图像做PSNR对比。我一般会把NLM的结果和medfilt2、imgaussfilt做对照NLM的PSNR通常比高斯滤波高1~3dB。如果你复现时发现PSNR提升不明显优先检查h的单位、噪声标准差是否用小数以及搜索窗口是否太小。把你写好的NLmeansfilt.m放进MATLAB路径准备好测试图调用一行代码就能跑起来。这个工程虽然只有一个文件但把非局部均值的距离度量、权重映射、归一化流程完整走了一遍足够当作进一步研究非局部算子的起点。想要更接近论文效果继续加上积分图加速和逐像素噪声方差估计就行。本文还有配套的精品资源点击获取