简介面向图像处理科研与工程应用的MATLAB滤波去噪资料包整合了双线性滤波、Kirsch滤波、逆滤波、双边滤波、同态滤波、小波滤波、约束最小平方滤波、非线性复扩散滤波、Lee滤波、Gabor滤波、Wiener滤波、Kuwahara滤波、Beltrami流滤波、Lucy-Richardson滤波、NonLocalMeans滤波等主流算法适合正在学习图像去噪、复原或开展算法对比的本科生、研究生及相关开发者使用。压缩包内共217个文件以215个m脚本为主另含1个txt说明与1个fig交互界面文件txt提供使用指引fig可直接打开交互演示界面能从算法调用、参数调试到GUI演示形成完整流程整体仅77KB轻量易部署。平台显示已有558人浏览学习。借助这些实现读者能快速获得各类滤波器的可运行代码、统一测试入口及可视化框架便于逐一复现算法效果、对比去噪性能并迁移到自己的图像处理课题中。1. 滤波这行当MATLAB 代码包里到底藏了多少种玩法把这份代码包里十几个滤波脚本全部改造成同一个测试接口之后我的第一个结论可能和大多数人想的不一样在低噪声密度、单一高斯噪声的场景下花力气实现的 Richardson-Lucy、Beltrami 流这类高级滤波和最简单的均值滤波拉不开本质差距真正让它们分出胜负的是混合噪声、纹理区域保真和边缘过冲这三件事。代码包里的Draw_Function_GUI.m、Run_Draw.m和一堆编号脚本3-2.m、15-1.m、17-4.m覆盖了从空间域到频域、从去卷积到变分法去噪的完整链路适合做课程设计、论文复现也适合在 MATLAB 图像处理大作业里直接抽某个模块改参数跑实验。下面按滤波器族分组拆解重点说清每个方法的适用边界和 MATLAB 实现时的关键参数最后给出一套横向评估脚本。2. 空间域滤波怎么选均值、中值到超限邻域与 Kuwahara 的取舍逻辑2.1 邻域统计类滤波的统一实现框架空间域滤波的本质是对每个像素的邻域做加权统计差别只在权重和统计量上。代码包里3-4.m、10-2.m这类脚本属于同一族可以用一个统一框架来跑实验避免每个文件复制粘贴改参数% filter_survey.m % 统一测试入口同一张图、同一噪声跑多个空间域滤波器 img im2double(imread(lena.png)); if size(img,3) 3, img rgb2gray(img); end rng(2024); noisy imnoise(img, gaussian, 0, 0.01); % 高斯噪声方差 0.01 % 1) 均值滤波3x3 邻域 h fspecial(average, [3 3]); out_avg imfilter(noisy, h, replicate); % 2) 中值滤波3x3 邻域 out_med medfilt2(noisy, [3 3], symmetric); % 3) 超限邻域滤波中心与邻域均值差超过阈值才替换 th 0.05; mu imfilter(noisy, h, replicate); diffmap abs(noisy - mu); out_lim noisy; out_lim(diffmap th) mu(diffmap th);参数说明imfilter的边界选项replicate对图像边缘的处理方式是复制边缘像素比默认补零少引入黑边medfilt2的symmetric则是对边界做镜像扩展在边缘纹理场景下更稳。超限邻域滤波的阈值th是关键它决定「多大差异算噪点」阈值太小会把真实边缘也抹掉太大则退化成不滤波。代码包里对这组实验的处理方式是固定噪声方差再对th做 0.01 到 0.1 的扫描生成一组对比图直接在Run_Draw.m里以子图网格展示。2.2 Kirsch 方向模板在滤波里的角色Kirsch 滤波经常被误认为是纯边缘检测算子实际上代码包把它归入滤波模块是有道理的它对八个方向分别做卷积取响应最大的方向作为该像素的「主方向」然后用方向自适应权重做平滑。这属于边缘保持滤波的一个朴素实现比后来基于偏微分方程的方法计算量小一个量级。% kirsch_dir.m % 8 个方向模板只列前两个示例其余按 45 度旋转生成 k1 [5 5 5; -3 0 -3; -3 -3 -3]; % 北方向 k2 [5 5 -3; 5 0 -3; -3 -3 -3]; % 东北方向 kernels zeros(3, 3, 8); for d 1:8 kernels(:,:,d) imrotate(k1, (d-1)*45, crop, bilinear); kernels(:,:,d) kernels(:,:,d) / sum(abs(kernels(:,:,d)(:))); % 归一化 end resp zeros(size(noisy,1), size(noisy,2), 8); for d 1:8 resp(:,:,d) conv2(noisy, kernels(:,:,d), same); end [~, dirIdx] max(abs(resp), [], 3); % 每个像素的主方向编号这里imrotate用来生成旋转后的方向模板但插值会引入轻微误差我通常直接用预定义数组存 8 个方向。方向编号dirIdx的意义在于后续平滑时只沿主方向及其垂直方向做加权平均而不是全向平滑这样能保留线条类纹理。代码包里ssbm.m应该是「双边滤波 方向模板」的变体区别于 2.4 节的标准双边滤波。2.3 Kuwahara 滤波与局部统计量Kuwahara 滤波的核心思路把像素邻域切成四个重叠的子块分别计算均值和方差输出方差最小子块的均值。原理是方差小的区域更均匀噪声被抑制的概率更大同时边缘两侧的方差差异会让输出偏向某一侧从而保住边缘。% kuwahara_impl.m % 窗口尺寸 5x5分成 4 个 3x3 子块 imgP padarray(noisy, [2 2], replicate); out_kuwa zeros(size(noisy)); for i 3:size(imgP,1)-2 for j 3:size(imgP,2)-2 block imgP(i-2:i2, j-2:j2); b1 block(1:3, 1:3); b2 block(1:3, 3:5); b3 block(3:5, 1:3); b4 block(3:5, 3:5); means [mean(b1(:)) mean(b2(:)) mean(b3(:)) mean(b4(:))]; vars [var(b1(:)) var(b2(:)) var(b3(:)) var(b4(:))]; [~, idx] min(vars); out_kuwa(i-2, j-2) means(idx); end end这个双重循环在 MATLAB 里效率低我一般会在代码包里把它替换成nlfilter或blockproc来加速但循环版本的逻辑最直观。注意padarray的replicate选项和imfilter一致避免边界出现异常子块。方差选择还有一个扩展版本就是在方差里加一个正则项vars lambdalambda越大输出越接近全均值越小越倾向于保留边缘纹理代码包15-1.m对lambda做了 0、0.01、0.1 三组对比。2.4 双边滤波与 Lee 滤波权重函数的设计差异双边滤波的空间域权重和灰度域权重相乘灰度差越大权重越小这使它在平坦区域等效于均值滤波在边缘处等效于「边缘另一侧不参与平均」。Lee 滤波则基于局部均值和方差做线性最小均方误差估计它对乘性噪声如 SAR 图像有理论最优性但对高斯噪声的表现不如双边滤波直观。滤波器核心统计量适用噪声边缘保持主要缺陷均值滤波邻域算术均值高斯弱边缘模糊中值滤波邻域中位数脉冲/椒盐中细线纹理丢失超限邻域均值差值阈值高斯中阈值敏感Kuwahara子块最小方差均值高斯中强块效应双边滤波空间灰度联合权重高斯强梯度反转伪影Lee 滤波局部方差加权估计乘性中强噪声模型不符时失效这里把六种空间域滤波器放在同一张表里对比代码包里的3-2.m应当还包含六抽头滤波即六系数 FIR 水平/垂直分离滤波常用于视频去隔行前后处理在静态图上可以看作一个 6×1 的平滑核与均值滤波效果接近但具有更好的频率选择性。空间域这一组跑完后基本能看出规律噪声若是脉冲型中值滤波优先若是高斯型且边缘保护要求高直接上双边或 Kuwahara若不想调参又要稳定输出超限邻域滤波是性价比最高的选择。3. 频域与逆问题Wiener 滤波、约束最小平方和同态滤波的实现边界3.1 逆滤波为什么会在 MATLAB 里直接崩掉频域滤波的基本操作是G F .* H N其中H是退化传递函数N是噪声频谱。最简单的恢复方式是把观测频谱除以H即逆滤波% inverse_filter.m F fft2(noisy); H fspecial(motion, 21, 11); % 运动模糊核 Hf psf2otf(H, size(noisy)); G F ./ Hf; % 直接相除 out_inv real(ifft2(G));这段代码跑完输出往往是一幅布满亮暗斑点的图。原因是Hf在高频区域接近零噪声被放大到无穷人类视觉对高频噪声又极敏感。我在代码包里看到17-4.m处理逆滤波时采取的补偿措施是对Hf做截断把绝对值小于某个阈值的频点置为一个固定小量而不是让它参与除法。这一步牺牲了高频细节但保住了整体灰度动态范围。3.2 Wiener 滤波与约束最小平方滤波的 MATLAB 参数Wiener 滤波在频域的形式是G conj(H) .* F ./ (abs(H).^2 1/SNR)MATLAB 的deconvwnr把它封装成了三种调用方式传噪声信号比、传自相关函数、传噪声功率谱。代码包里的标准用法是第三种% wiener_deconv.m PSF fspecial(gaussian, [7 7], 2); wnr1 deconvwnr(noisy, PSF, 0.01); % 固定 NSR noise_var 0.01; sig_var var(noisy(:)); wnr2 deconvwnr(noisy, PSF, noise_var / sig_var); % 比值形式deconvwnr的第一个参数是退化图像第二个是点扩散函数第三个是信噪比参数。比值形式比固定 NSR 更稳因为var(noisy(:))包含了信号和噪声的联合方差当噪声方差超过信号方差时比值会远大于 1等效于更强的高频衰减防止振铃。约束最小平方滤波deconvreg的用法类似但多了拉格朗日乘子alpha% reg_deconv.m [out_reg, reg_arg] deconvreg(noisy, PSF, 0.4); % alpha0.4reg_arg返回算法实际使用的正则化参数如果输出图像仍有振铃就把alpha调大一个数量级如果图像发糊就把alpha调小。相比 Wiener 滤波约束最小平方对噪声模型不敏感适用面更广代价是参数调节没有解析解只能靠观察。3.3 同态滤波照射反射模型下的动态范围压缩同态滤波把图像建模为照射分量乘反射分量取对数后变成加法再用高通滤波压缩低频照射、保留高频反射% homomorphic.m img_log log(noisy eps); F fft2(img_log); [rows, cols] size(F); u linspace(-0.5, 0.5, rows); v linspace(-0.5, 0.5, cols); [Dx, Dy] meshgrid(v, u); D sqrt(Dx.^2 Dy.^2); gh 1.2; gl 0.4; cutoff 0.3; H (gh - gl) .* (1 - exp(-D.^2 / (2*cutoff^2))) gl; G F .* ifftshift(H); out_homo exp(real(ifft2(G)));gh是高频增益gl是低频增益cutoff是高斯高通截止频率。ifftshift这一步容易漏meshgrid生成的频率原点在左上角而H的构造假设原点在中心必须用ifftshift对齐。同态滤波适合光照不均匀的图像比如一张半边暗半边亮的照片用它压缩亮度差异后后续阈值分割会稳定很多。但它对加性噪声没有建模能力噪声强度高时会把噪声当作反射分量放大所以代码包里的顺序是先把噪声做一次中值预滤波再进同态滤波。4. 反卷积与保边去噪Richardson-Lucy 迭代、NoLocalMeans 与非线性复扩散4.1 Richardson-Lucy 的迭代语义和阻尼参数Richardson-Lucy 最早用于天文图像恢复它假设噪声服从泊松分布用期望最大化推导出迭代格式。MATLAB 中的函数名是deconvlucy它接受的参数远比常见教程里写得多% rl_deconv.m PSF fspecial(gaussian, [9 9], 3); DAMPAR 0.02; % 阻尼阈值控制噪声放大 WEIGHT ones(size(noisy)); % 坏点权重如坏像素置 0 out_rl deconvlucy(noisy, PSF, 10, DAMPAR, WEIGHT);迭代次数默认是 10这个值对结果影响极大。迭代少则模糊残留迭代多则出现「振铃 斑点噪声」因为 RL 对高频噪声的放大是随迭代单调增加的。DAMPAR的作用是限制每一步修正量修正值低于该阈值就不更新WEIGHT则能屏蔽掉坏像素对周围区域的干扰。代码包14-1.m对迭代次数做了 5、10、20、50 四组对照结论是泊松噪声场景下 10 次左右性价比最高高斯噪声下 RL 并不占优应该考虑 Wiener。4.2 NoLocalMeans 的搜索窗口与平滑强度NoLocalMeansNLM和双边滤波的区别在于相似度计算的对象双边滤波逐像素比灰度NLM 逐块比局部邻域结构。MATLAB 在较新版本里直接提供imnlmfilt% nlm_filter.m out_nlm imnlmfilt(noisy, DegreeOfSmoothing, 0.08, ... SearchWindowSize, 21, ... ComparisonWindowSize, 7);DegreeOfSmoothing控制指数核的衰减速度值越大平滑越强一般取噪声标准差的 0.5 到 1.5 倍SearchWindowSize是搜索窗越大越能找到更多相似块但计算量平方增长ComparisonWindowSize是块大小。NLM 对纹理区域的保护能力明显强于双边滤波主要代价是速度21×21 的搜索窗在大图上跑几十秒很正常。代码包给出的优化思路是先降采样加速最后对结果做一次边缘保持上采样。4.3 非线性复扩散滤波的迭代格式非线性复扩散滤波是 Beltrami 流之外的另一个变分思路把扩散系数从实数扩展到复数虚部提供边缘增强。常见半隐式离散格式如下% complex_diffusion.m phi noisy 1i * 0.1 * noisy; % 初值加入虚部 lambda 0.15; dt 0.2; alpha 0.2; % kappa 控制边缘敏感度 kappa 2.0; for iter 1:15 g 1 ./ (1 abs(imag(phi)) / kappa); % 边缘停止函数 lap del2(phi); % 拉普拉斯算子 phi phi dt * (lambda * del2(g .* real(phi)) 1i * alpha * lap); end out_cdiff real(phi);这里的del2在 MATLAB 中计算拉普拉斯时已包含 0.25 的归一化系数所以dt的取值范围与显式差分不同一般 0.1 到 0.3 之间稳定。虚部imag(phi)在迭代中积累的是二阶导数信息边缘处虚部较大对应的g变小扩散被抑制平坦区域虚部接近零扩散正常进行最终效果是在去噪的同时增强边缘。5. Beltrami 流与小波阈值去噪两类非线性方法的 MATLAB 实现路径5.1 Beltrami 流把它理解为「图像即曲面」的几何流Beltrami 流把灰度图像看作嵌入高维空间中的二维流形去噪过程就是让这个流形按面积最小化方向演化。相比 Perona-Malik 各向异性扩散Beltrami 流的扩散张量由图像梯度构造保边能力更好且对参数不敏感。代码包里没有现成的函数通常自己写显式迭代% beltrami_flow.m u noisy; dt 0.05; iter 20; sigma 1.0; gauss fspecial(gaussian, [5 5], sigma); for t 1:iter ux imfilter(u, [-1 0 1], replicate) / 2; uy imfilter(u, [-1; 0; 1], replicate) / 2; uxx imfilter(u, [1 -2 1], replicate); uyy imfilter(u, [1; -2; 1], replicate); uxy imfilter(ux, [-1; 0; 1], replicate) / 2; W2 1 ux.^2 uy.^2; % Beltrami 流离散格式1/sqrt(W2) * (拉普拉斯 - 二次项) lap uxx .* (1 uy.^2) - 2 .* ux .* uy .* uxy uyy .* (1 ux.^2); u u dt .* lap ./ W2.^2; end out_beltrami u;W2是度量张量的行列式相关量lap ./ W2.^2实现了非线性扩散。dt超过 0.1 时数值不稳定表现为迭代若干步后图像出现棋盘格检测方法很简单打印每次迭代前后的max(abs(u(:) - u_old(:)))这个值如果跳变到 1e-2 量级就说明步长过大。Beltrami 流的优点是迭代次数增加时图像细节不会像均值滤波那样持续模糊而是在保边和平滑之间收敛到一个平衡点。5.2 小波阈值去噪的三层参数小波去噪是分解、阈值、重构三步。MATLAB 的wavedec2负责分解wthresh负责阈值% wavelet_denoise.m [coefs, books] wavedec2(noisy, 3, db4); % 三层分解db4 小波 [thr, sorh, keepapp] ddencmp(den, wv, noisy); denoised wdencmp(gbl, coefs, books, db4, 3, thr, sorh, keepapp); % 自定义硬阈值写法 c_hard coefs; c_hard(abs(c_hard) thr) 0; out_wav waverec2(c_hard, books, db4);ddencmp自动估计全局阈值sorh为s表示软阈值h表示硬阈值。软阈值把所有系数往零收缩视觉更平滑但没有硬阈值锐利硬阈值保留大于阈值的系数细节好但可能出现小幅震荡。第三层分解后近似分量keepapp不要做阈值处理否则图像整体亮度会被压低。小波去噪在混合噪声场景下的表现比空间域滤波稳定因为它天然把信号和噪声按频带分开高频层阈值处理对低频分量的扰动极小。5.3 小波与 Beltrami 在代码包里的搭配顺序代码包的习惯是先用小波去掉高频高斯噪声再用 Beltrami 流处理剩余的结构性伪影。这个顺序有实际依据小波阈值对孤立噪声点有效但对纹理边缘处的噪声块处理不干净而 Beltrami 流恰好擅长修复这类边缘附近的不规则噪声。反过来顺序不行先跑 Beltrami 会把高频细节磨掉一部分小波再去噪时可能把真正的纹理误判为噪声。6. 用一份测试脚本快速横向评估滤波结果代码包里Run_Draw.m的职责是统一出图但只出图很难量化对比。我在工程实践中的做法是补一个eval_filters.m用 PSNR 和 SSIM 自动打分% eval_filters.m results struct(); results.avg filter_eval(noisy, out_avg, img); results.med filter_eval(noisy, out_med, img); results.nlm filter_eval(noisy, out_nlm, img); results.rl filter_eval(noisy, out_rl, img); ssim_table struct2table(results, RowNames, {PSNR, SSIM}); function scores filter_eval(original, noisy, clean) scores [psnr(noisy, clean); ssim(noisy, clean)]; end真正调参时只看两个数不够我一般会在15-1.m的基础上加一个残差图imshowpair(noisy, clean, diff)红色区域代表滤波过度绿色区域代表残留噪声。若红色集中在边缘说明滤波器保边太强适当调大 NLM 的DegreeOfSmoothing或小波的软阈值若绿色均匀散布说明阈值偏低把ddencmp返回的thr乘以 1.2 再跑一轮。最后把每组实验的 PSNR、SSIM、运行时间和参数量填进表格哪个滤波器值得继续调参哪组参数已经逼近上限直接看哪一行震荡最大。本文还有配套的精品资源点击获取