中值滤波MATLAB代码实战：从原理到自适应去噪

简介这是一份面向MATLAB图像处理初学者的中值滤波示例代码演示如何利用medfilt2函数有效去除椒盐噪声。压缩包内仅含1个m脚本文件整体大小仅640B代码精简、结构清晰适合快速入门。目前已有1997人浏览学习广泛用于图像去噪实战练习。通过这份代码读者可以掌握灰度图像读取、3×3邻域中值滤波调用、滤波前后对比显示等完整流程同时理解注释规范对提升代码可读性的重要作用。这份代码虽小但完整覆盖了从非线性滤波原理到MATLAB编程实现的关键环节可作为课堂演示、课后作业或信号处理交叉应用的起点模板。1. 简单的中值滤波MATLAB代码先把“去椒盐噪声”这件事做对做图像处理的工程师都有过这种经历一张采集卡刚送进来的灰度图满屏黑白噪点均值滤波越抹越糊线性滤波又压不掉那些“孤点”。这时候最省事、最容易解释给别人的办法就是中值滤波。它不依赖线性系统那套频率响应分析逻辑一句话就能讲完把窗口里的像素排序取中间那个数代替当前点。这个非线性操作对脉冲噪声椒盐噪声、传感器坏点、扫描条纹几乎是量身定做而且MATLAB里既有现成函数也适合自己写一版去理解边界效应。这篇文章就从“为什么中值能去噪”讲到“怎么写出一段简单但能改参数的中值滤波MATLAB代码”再补窗口选择、边界填充、彩色图像处理和自适应变体。适合新手照着复现也适合老手在换语言、写嵌入式版本前用MATLAB快速验证算法行为。2. 中值滤波原理与MATLAB内置函数选择2.1 中值为什么能压住椒盐噪声排序统计量的直观解释椒盐噪声的特点是“少数像素值极端偏离邻域”。一个3x3窗口里有8个正常灰度比如100左右和1个盐噪声灰度255均值会把中心点拉高到约117而中值排序后取第5个值仍然是100附近的正常值。关键就在于中值对离群点不敏感——它只关心“比它小的一半”和“比它大的一半”而均值关心“所有值加起来”一个极端值就能把均值拖走。这个性质在数学上对应的是中值是L1损失的极小值点均值是L2损失的极小值点。L1对离群点权重是线性的L2是平方级的所以同样邻域内混入坏点中值估计比均值稳健得多。实际传感器数据里的坏像素、传输误码、数字化时的尖刺本质上都是这种“少数极端值”这正是中值滤波在图像预处理器里长期被放在最前面一级的原因。2.2 MATLAB现成函数medfilt2与medfilt1怎么选MATLAB图像处理工具箱里二维灰度图用medfilt2一维信号用medfilt1。最常见图像去噪场景直接调用medfilt2就够了img imread(cameraman.tif); if size(img,3) 3 img_gray rgb2gray(img); else img_gray img; end img_noise imnoise(img_gray, salt pepper, 0.1); % 3x3窗口默认边界补0 img_med medfilt2(img_noise, [3 3]); % 比较去噪前后信噪比 psnr_after psnr(img_med, img_gray); fprintf(处理后 PSNR %.2f dB\n, psnr_after); imshowpair(img_noise, img_med, montage);medfilt2的第二个参数[3 3]指窗口尺寸必须是奇数单数写法3等价于[3 3]。它默认使用边界补零zeros在图像边缘处窗口覆盖不到的区域补零会引入一圈黑色假边缘。如果图像四边内容重要可以改成symmetric让边缘镜像扩展。内置函数适用维度常用调用形式关键参数medfilt2二维灰度/二值图像medfilt2(I, [m n])窗口尺寸边界选项 padoptmedfilt1一维时序信号medfilt1(x, k)窗口长度 k填充方式 truncatemedfilt3三维体数据/彩色堆叠medfilt3(V, [m n p])窗口尺寸输入需为 double 或 single一维信号处理场景比如传感器时序数据去毛刺用medfilt1(x, 5)窗口长度5返回与x等长的滤波结果。它的边界处理由truncate、zeropad等选项控制默认truncate直接截掉两端的非完整窗口输出长度不变但两端有轻微误差。%% 一维信号的脉冲噪声去除示例 fs 1000; t (0:0.001:0.1); clean sin(2*pi*50*t); noisy clean; idx randperm(length(noisy), 15); noisy(idx) randn(15,1)*3; % 模拟尖刺 med1 medfilt1(noisy, 5); figure; plot(t, clean, k-, t, noisy, r., t, med1, b-); legend(原始, 带刺信号, 中值滤波);medfilt1的长度参数同样要求奇数。选5表示取当前点前后各2个点排序取中值。点数取太大会把真实的窄脉冲一起抹掉取太小又滤不净刺。下面第4章会专门谈窗口参数的定法。3. 手写一个简单的中值滤波MATLAB代码3.1 从三重循环到向量化三种可抄写法内置函数能满足大部分需求但理解中值滤波的本质还是得手写一版。第一种写法最直白对每个像素抠出邻域排序取中值。代码如下function out my_medfilter_loop(I, win) [H, W] size(I); r floor(win/2); out zeros(H, W, like, I); % 边界补零后的扩展图像 padded zeros(H 2*r, W 2*r); padded(r1:Hr, r1:Wr) I; for i 1:H for j 1:W % 取当前邻域并排序中心位置取中值 block padded(i:iwin-1, j:jwin-1); sorted sort(block(:)); out(i,j) sorted(ceil(win*win/2)); end end end调用out my_medfilter_loop(noisy_img, 3);。这套代码胜在逻辑与公式一一对应适合做算法说明、写论文伪代码对照、以及移植到C语言。缺点是三重循环在1MP图像上用3x3窗口大约要0.5秒以上大图会觉得卡MATLAB里逐像素循环不是性能预期内的用法。如果想提速又可以不改逻辑用im2col把每个像素的邻域向量化成矩阵一次sort矩阵消除内层循环function out my_medfilter_im2col(I, win) [H, W] size(I); r floor(win/2); padded padarray(I, [r r], replicate); % 边界复制扩展 B im2col(padded, [win win], sliding); % 每列是一个邻域 sortedB sort(B, 1); medIndex ceil(win*win/2); out reshape(sortedB(medIndex, :), H, W); out cast(out, like, I); endim2col把滑动窗口变成矩阵列B的列数等于输出像素个数。对sort(B,1)按列排序后取中位数所在行即得到整幅图的中值结果。这个版本速度远超三重循环适合课堂演示“向量化给你带来什么”。3.2 边界三种处理法zeros、replicate、symmetric怎么切上面两段代码里边界策略分别是补零和replicate。实际写代码时边界策略会显著影响边缘像素的滤波结果。见表边界模式实现方式典型场景边缘效果zeros补零padarray(I,[r r],0)暗背景显微图、文档扫描图像四周变暗PSNR下降replicate复制边缘padarray(I,[r r],replicate)自然图像、通用场景边缘有轻微延长感更常用symmetric镜像padarray(I,[r r],symmetric)周期性纹理、航拍图边界衔接最自然计算略多我一般建议默认用symmetric或replicate不要用零填充。原因很实际零填充相当于在图像四周人造了一圈黑色像素排序结果会被拉低边缘细节被破坏。但如果是显微镜深色背景且噪声集中在亮目标内部零填充与重建背景更接近反而选零也可接受。所以边界策略没有绝对优劣要在具体图像上跑一遍对比。第二版手写代码用padarray配replicate是最省心的。若想比较三种边界模式对同一图的差异可以写一个循环modes {zeros, replicate, symmetric}; for k 1:3 p padarray(noisy, [1 1], modes{k}); B im2col(p, [3 3], sliding); res sort(B,1); im_f reshape(res(5,:), H, W); fprintf(%s: PSNR %.2f dB\n, modes{k}, psnr(im_f, clean)); end这段代码会把三种边界策略的输出信噪比打印出来根据数值而不是经验选策略。参数说明第5行res(5,:)表示3x3窗口共9个值排序后取第5个即中值padarray的第二个参数[1 1]表示上下左右各扩展1个像素。3.3 运行一段手写代码验证正确性写完代码第一件事是先构造一个已知结果做断言而不是立刻上带噪声的大图。例如%% 手工构造一幅受脉冲噪声污染的横条纹图 I repmat(20:20:200, 50, 1); % 灰度渐变图 J I; J(10:15, 10:15) 255; % 手工撒一块盐噪声 Iref medfilt2(I, [3 3]); Icor my_medfilter_im2col(J, 3); assert(isequal(Iref, Icor), 手写结果与内置函数不一致);这个断言能同时验证索引计算、边界处理和排序位置三个环节。如果isequal不通过先在3x3窗口下比较局部结果找索引偏差再检查sort取的是ceil(win*win/2)而不是floor。小结构验证通过后再上1MP真实图用imshowpair看细节问题基本只在边界处。4. 中值滤波参数设置与三个高频坑4.1 窗口大小怎么选3x3、5x5、7x7的效果边界窗口是中值滤波唯一的超参数直接影响细节保留和去噪强度。下表给出同幅图像在不同噪声密度下的经验范围噪声密度推荐窗口说明 10%3x3保留细节已有明显视觉效果10% ~ 30%5x5在细节损失与去噪间折中30% ~ 60%7x7 或 9x9能压住颗粒但细纹理基本消失 60%二次迭代或自适应普通中值滤波本身难以胜任窗口过小的表现是噪声没滤干净窗口过大则图像变成“漫画感”。验证方法是画一条水平剖面线比较滤波前后同一行的灰度跳变figure; plot(clean(50,:), k); hold on; plot(noisy(50,:), r.); plot(med3(50,:), b); legend(原图, 带噪, 5x5中值);横坐标是列号纵坐标是灰度值。观察蓝线相对红线是否平整以及边沿处从平到凸是否滞后。滞后越明显说明窗口吞掉了原本的锐利边沿应缩小窗口或改用后文的自适应方案。4.2 常见坑uint8溢出、别用均值、迭代要克制第一个坑是数据类型。imread读进来是uint8直接对uint8做加法或减法会产生溢出回绕比如150 - 80在uint8下不是70而是接近255。中值滤波本质只做排序和索引uint8本身不会出错但一旦在滤波前做了图像差分或添加噪声的操作务必先转double再计算最后再转回显示I_d double(I); noise_mask rand(size(I_d)) 0.1; sp I_d; sp(noise_mask) 255; out_d my_medfilter_im2col(sp, 3); out uint8(out_d); % 显示前转换上面的代码先转double再撒盐点避免了uint8下的回绕算术第4行把10%的像素直接置成255模拟盐噪声滤波完成后才转回uint8做显示或imwrite。第二个坑是对椒盐噪声用均值滤波。前面原理部分已经说明均值对极端值敏感3x3下一个255盐点就能把邻域均值拉高近20个灰度级。如果只有工具箱没有图像处理模块、只能filter2做通用卷积也不要拿卷积模板去近似中值。中值没有等价的线性卷积核。第三个坑是迭代滥用。有些人看到中值滤波效果不够彻底就连续跑三次5x5。第一次中值能把孤立脉冲点去掉第二次后边缘持续被“修剪”细结构逐步平滑成块。更合理的做法是先跑一次检测残余噪声再只对噪声像素做第二次滤波见第5章。迭代中值滤波在数学上会收敛到“根信号”但根信号往往严重丢失细节和原图已经不像了。4.3 彩色图与多通道从灰度中值扩展到RGB彩色中值不能简单对三个通道分别滤波然后拼接因为这样会引入原本不存在的颜色。举例红色通道的噪声使某像素R值跳成255单独滤波后R被压回正常但G、B通道没被压到同样程度结果是原像素颜色已经改变还出现了新的色偏。常见做法是把RGB三个分量排序后按矢量中值或联合中值来处理或者退一步实用工程里先转YCbCr对亮度Y通道做中值滤波色度通道只对脉冲噪声严重时才滤波这样既控制成本又避免色偏img_rgb imread(peppers.png); ycbcr rgb2ycbcr(img_rgb); Y ycbcr(:,:,1); Y_f medfilt2(Y, [3 3]); % 对亮度做中值 ycbcr(:,:,1) Y_f; out_rgb ycbcr2rgb(ycbcr);如果确需对彩色图三通道同等地滤除彩色椒盐噪声至少用medfilt2对每个通道分别处理后检查通道间是否出现“边缘错位”并在显示时对比原图颜色直方图是否有通道偏移。严格意义上矢量中值滤波按欧氏距离选邻域内到其他点距离最小的像素作为输出能保持颜色相关性但运算量明显上升除非是离线数据处理一般优先转亮度域处理。5. 用中值滤波做图像增强前处理一个自适应窗口技巧中值滤波在完整图像处理流程里通常只是第一步后续常接对比度拉伸、锐化或结构分析。固定窗口面临两难窗口小了高密度噪声压不干净窗口大了薄纹理被磨平。一个可复现的折中思路是先检测脉冲噪声位置只对可疑噪声点做自适应窗口滤波。检测思路不依赖训练直接用中值的定义脉冲噪声在局部邻域里往往对应极大值或极小值。若当前像素灰度比邻域中值高或低出阈值判定为疑似噪声否则直接保留原值。这样正常像素不过滤细节几乎原样保留只有坏点被替换。function out adaptive_medfilter(I, win_max, thresh) [H, W] size(I); I double(I); out I; r_max floor(win_max/2); for i 1r_max : H-r_max for j 1r_max : W-r_max for w 3:2:win_max r floor(w/2); block I(i-r:ir, j-r:jr); med median(block(:)); cur I(i,j); % 超过阈值才认定是噪声否则保留原值 if abs(cur - med) thresh if w win_max out(i,j) med; % 已到最大窗口直接用中值 end continue; else out(i,j) cur; break; end end end end out uint8(out); end调试时先固定窗口上限为7阈值取25。把输出与原图相减差值图的亮斑就是未滤净的噪声区若亮斑沿边缘连成线说明阈值太小把正常边缘误判成噪声应该上调到3040。若噪声点密集成片则把win_max提到9并保留两层循环。该函数只适合灰度图。需要做三通道时先转YCbCr只对亮度分量调用它色度分量不处理需要在意边界时先padarray做symmetric扩展再进入循环。比较自适应与固定窗口用psnr与ssim两个指标同时看前者反映像素误差后者反映结构保持。若自适应版像素误差略高但SSIM明显更高说明它在压制噪声的同时保留下了更多边缘对后续分割、配准这类任务通常更有利。本文还有配套的精品资源点击获取