尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

自适应中值滤波原理详解与Matlab去噪仿真实现

自适应中值滤波原理详解与Matlab去噪仿真实现 简介面向图像处理初学者、课程设计学生以及需要快速复现去噪算法的研究人员这份资源提供了基于MATLAB的自适应中值滤波图像去噪完整仿真方案可有效应对椒盐噪声污染并在去噪同时保留更多边缘细节是理解自适应滤波机制的良好起点。压缩包共3个文件包含核心算法脚本.m、演示操作录像.avi与测试图像.jpg整体仅2.14MB轻量便携。脚本实现主控制逻辑与可视化输出录像逐步展示从环境配置到运行出图的完整过程测试图则方便直接验证算法效果。运行前需将MATLAB当前文件夹切换至工程路径建议使用2021a或更高版本配套视频可辅助新手避开常见环境问题。已有712人学习下载适合课程设计、论文实验或算法对比等场景。借助包内完整的源码、测试图与操作演示可深入理解自适应窗口尺寸随噪声密度动态变化的思想并灵活修改参数对比去噪效果。1. 自适应中值滤波为什么比普通中值滤波更值得仿真在图像去噪这个老话题里中值滤波一直是最容易上手、最不挑硬件的一类方法。但真把一张叠加了高密度椒盐噪声的图丢给固定窗口的中值滤波结果往往不是去噪而是把图像弄成一块块模糊的色斑。原因不复杂噪声密度高时固定窗口里的中值本身可能就是被污染像素排序取中值失去了意义。更麻烦的是中值滤波对边缘和细纹理的处理是“一视同仁”的——它分不清某像素是噪声还是本来就在细节区域。自适应中值滤波Adaptive Median FilterAMF解决的就是这个矛盾窗口大小和是否滤波都由局部统计特性决定。能确认是脉冲噪声的才动这个像素确认是边缘细节的主动放行。这个“先在噪声检测上做文章、再决定滤波强度”的思路在图像去噪这个方向上几乎是所有现代去噪算法从加权中值到非局部均值的雏形。用 Matlab 做这个算法的仿真核心不是调一个medfilt2而是自己把窗口扩展、噪声检测、双阈值判断这些逻辑一步步写出来。这篇帖子会从算法判定逻辑讲到 Matlab 代码实现再对比不同参数下的去噪效果最后聊怎么把仿真过程做成能放进 PPT 或演示视频的动态展示。适合做图像处理课程实验、毕设仿真以及刚开始接触自适应滤波但不想只看公式的工程师。2. 自适应中值滤波的两级判定逻辑从噪声检测到窗口扩展2.1 标准中值滤波的失效边界先明确中值滤波到底在算什么。对图像中的某个像素点 ((x,y))取一个大小为 (W \times W) 的邻域窗口窗口内所有像素按灰度排序取中间值作为当前像素的新灰度。Matlab 里一行g medfilt2(f, [3 3])就能跑但这一行隐藏了两个默认行为窗口大小固定且对窗口内的每个像素都做替换。固定窗口的隐患在噪声密度上来之后非常明显。假设椒盐噪声密度是 0.3一个 3×3 窗口内平均会有约 2.7 个噪声点。排序后的中值落在真实像素上的概率尚可但如果噪声密度到 0.53×3 窗口里一半左右是噪声排序取到的中值大概率本身就是噪声或受噪声影响的灰度滤波后图像会残留大量黑白点。更隐蔽的情况是边缘像素一个干净的边缘像素和一个被噪声污染的非边缘像素在固定窗口中各占 9 个灰度值排序后的输出几乎总会压平边缘造成细节丢失而这种丢失在视觉上表现得比噪声本身更难接受。自适应策略的核心就是不把每个像素都交给同一个窗口。它对每个像素做先判断、后处理判断不通过时自动扩大窗口直到满足条件或到达窗口上限。这套机制业界通常称为“双级判定”英文资料里常写为 two-level decision它也是后面代码实现的骨架。2.2 第一级判定区分噪声点和信号点先看窗口内部的灰度分布。设当前窗口为 (S_{xy})窗口半径为 (r)(z_{min})、(z_{max})、(z_{med}) 分别是窗口内灰度的最小值、最大值和中值。第一级判定规则如下若 (z_{min} z_{med} z_{max})说明中值本身不是极值即中值不是被椒盐噪声污染的亮度值。这种情况下窗口中心像素的原始灰度 (z_{xy}) 如果满足 (z_{min} z_{xy} z_{max})那就认为它是正常信号直接保留原始灰度不做滤波。若 (z_{xy}) 本身落在极值上即等于 (z_{min}) 或 (z_{max})判定为疑似噪声进入第二级判定。第一级判定的意义很多人会看漏它真正在做的是“保护非噪声极值”。图像中天然存在接近全黑或全白的像素如夜晚的路灯、白墙上的黑点这些像素的灰度在邻域中可能恰好是极小值或极大值。如果只按单一阈值判断这类像素会被误伤。所以第一级先看中值是否可靠再看中心像素是否真的是极值双条件同时满足才认为是噪声候选。2.3 第二级判定中值本身也是极值时怎么办当第一级判定发现 (z_{med}) 等于 (z_{min}) 或 (z_{max})说明当前窗口内噪声密度高排序后的中值已经不可靠。此时的做法不是直接滤波而是扩大窗口半径从 (r1) 重新执行第一级判定。具体流程从最小窗口半径 (r_{min})通常 1即 3×3 窗口开始。执行第一级判定若不通过则 (r \leftarrow r1)重新计算窗口灰度统计。若 (r) 已经达到最大窗口半径 (r_{max})常取 3 或 4即 7×7 或 9×9 窗口仍然不满足第一级条件说明局部区域噪声非常密集。此时不再无限扩大直接输出当前窗口的中值 (z_{med}) 作为替代值。这个“扩窗上限”的设计是自适应中值滤波的关键参数。(r_{max}) 设小了高密度噪声区域救不回来设大了虽然能覆盖更宽的空间但窗口跨越物体边界的概率增加容易把两个不同区域的灰度混在一起取中值产生伪纹理。第二级判定的另一条重要规则是当中心像素被判定为噪声、但窗口扩展后的中值仍落在极值上时也有资料主张用“相邻窗口中值”或“最近有效像素”来替代但这属于变体做法。课程和常见的毕业设计仿真里标准 AMF 多采用“扩窗后输出中值”的方式实现简单、效果可预期。2.4 AMF 与传统中值滤波的对比表维度传统中值滤波自适应中值滤波窗口大小固定动态调整从 3×3 自动扩至上限是否检测噪声不做区分全图统一处理两级判定只处理疑似噪声点边缘保护差边缘像素易被排序后的中值改写较好极值判定的目的是先识别再处理高密度噪声0.5残留大量斑点能通过扩窗提高恢复能力但也受 (r_{max}) 限制Matlab 内置函数medfilt2无内置需手动实现适用场景中低密度高斯/椒盐噪声做预处理高密度椒盐噪声、对细节有要求的去噪任务AMF 的代价是计算量显著上升每个像素都要做窗口统计和条件判断噪声密度越高窗口扩张越频繁排序次数也就越多。实际落地时通常先用快速排序算法如sort的底层实现来减少排序开销后面代码实现里会用这个思路。3. Matlab 实现从零手写自适应中值滤波函数3.1 灰度图预处理与噪声类型选择写滤波函数前先把测试图准备好。Matlab 里最常用的实验图是自带的camerman.tif也可以用imread读入自己的图像。关键是第一步要把图像转成double类型否则灰度的加减比较在uint8下会发生溢出问题。clear; clc; close all; % 读取灰度图转为 double归一化到 [0,1] I imread(cameraman.tif); if size(I, 3) 3 I rgb2gray(I); end I im2double(I); % 关键后续所有操作都在 double 下做 % 加入椒盐噪声密度 0.3 In imnoise(I, salt pepper, 0.3);参数说明imnoise的第三个参数是噪声密度取值范围 0 到 1表示总像素中多少比例会被替换为 0 或 255归一化后为 0 或 1。密度 0.3 意味着三成像素被污染这个密度对普通 3×3 中值滤波已经有明显压力但自适应滤波大多能恢复得不错。im2double把灰度从 0-255 的整数区间映射到 0-1 的浮点区间窗口内排序和比较都更稳。注意因为imnoise用的是随机数生成器每次运行输出的噪声位置不同如果需要可复现结果建议先执行一次rng(0)固定随机种子。3.2 自写 min 和 max 的朴素自适应窗口版本下面是一个原理最直观的实现每个像素独立处理代码量短逻辑和上面讲的两级判定一一对应。function out adaptive_median_filter(I, r_min, r_max) % 输入 I 为 double 灰度图值域 [0,1] % r_min 最小窗口半径通常为 13×3 % r_max 最大窗口半径建议 3 或 4 [M, N] size(I); % 边界填充用边缘镜像避免窗口越界 pad r_max; I_pad padarray(I, [pad pad], symmetric); out zeros(M, N); for i pad1 : Mpad for j pad1 : Npad for r r_min : r_max % 当前窗口半径下的局部块 block I_pad(i-r:ir, j-r:jr); zmin min(block(:)); zmax max(block(:)); zmed median(block(:)); zxy I_pad(i, j); % 第一级判定中值是否为极值 if zmin zmed zmed zmax % 第二级判定当前像素是否是疑似噪声 if zmin zxy zxy zmax out(i-pad, j-pad) zxy; % 是信号直接保留 else out(i-pad, j-pad) zmed; % 是噪声输出中值 end break; else % 中值不可靠若已到最大半径只能用当前中值兜底 if r r_max out(i-pad, j-pad) zmed; end end end end end end逻辑说明外层r循环从最小半径到最大半径递增内层对每个半径做窗口统计。break语句保证一旦满足第一级判定条件就不再扩大窗口。这里容易踩坑的是padarray的symmetric选项——它把边缘像素做镜像填充比补零好得多。如果补零图像四周边界会出现一圈由 0 引发的伪边缘滤波后视觉上会有一圈黑框。还有一点zmin zmed zmed zmax是严格小于不能取。在后文 3.3 会解释为什么会让平坦区域的像素完全走错误分支。3.3 优化版一维排序、复用统计量、避免窗口全量排序朴素版本的问题是慢。对 256×256 图像在密度 0.5 噪声下每个像素可能要重复计算多个窗口的median。如果实验图是 1024×1024朴素版本跑一次要几十秒根本不适合做参数扫描。优化思路分三层利用sort一次得到排序后的整行数据同时提取zmin、zmax、zmed避免min、max、median三次扫描。外层扩窗时上一轮窗口包含在下一轮窗口内。可以维护窗口内灰度的直方图扩窗时只增量更新边界新增像素而不是重新对整块排序。窗口半径递增时复用上一轮的排序结果理论上扩窗代价从 O(W^2 log W) 降到 O(W)但在 Matlab 里向量化的收益往往比手动维护数据结构的收益更明显所以折中做法是保留全窗口排序只是用sort(block(:))一次性取值。function out adaptive_median_filter_fast(I, r_min, r_max) % 优化版利用 sort 一次取三值减少重复扫描 [M, N] size(I); % 对中心像素做扩展时先预留输出矩阵 out zeros(M, N); pad r_max; I_pad padarray(I, [pad pad], symmetric); for i pad1 : Mpad for j pad1 : Npad for r r_min : r_max block I_pad(i-r:ir, j-r:jr); sorted sort(block(:)); % 升序 zmin sorted(1); zmax sorted(end); zmed sorted(ceil(length(sorted)/2)); zxy I_pad(i, j); if zmin zmed zmed zmax if zmin zxy zxy zmax out(i-pad, j-pad) zxy; else out(i-pad, j-pad) zmed; end break; else if r r_max out(i-pad, j-pad) zmed; end end end end end end核心改动是sorted sort(block(:))这一行block(:)将二维窗口拉成一维列向量sort升序排列后最小值是第一个元素最大值是最后一个中位数是正中间那个。注意在元素个数为偶数时中值定义有两种这里统一用ceil效果偏差在图像上肉眼基本看不出来。3.4 与传统 medfilt2 的对比脚本实现完滤波器立刻就能做对比实验。下面脚本同时跑出三种结果普通 3×3 中值、自适应中值r_min1, r_max3、自适应中值r_min1, r_max4并计算 PSNR 和 SSIM。% 直接调用上面写的函数 I_med3 medfilt2(In, [3 3]); % 传统中值滤波 I_amf3 adaptive_median_filter_fast(In, 1, 3); % 自适应最大 7×7 I_amf4 adaptive_median_filter_fast(In, 1, 4); % 自适应最大 9×9 % 计算峰值信噪比 PSNR psnr_med3 psnr(I_med3, I); psnr_amf3 psnr(I_amf3, I); psnr_amf4 psnr(I_amf4, I); % 计算结构相似性 SSIM ssim_med3 ssim(I_med3, I); ssim_amf3 ssim(I_amf3, I); ssim_amf4 ssim(I_amf4, I); fprintf(PSNR: medfilt2%.2f dB, AMF(r_max3)%.2f dB, AMF(r_max4)%.2f dB\n, ... psnr_med3, psnr_amf3, psnr_amf4); fprintf(SSIM: medfilt2%.4f, AMF(r_max3)%.4f, AMF(r_max4)%.4f\n, ... ssim_med3, ssim_amf3, ssim_amf4);说明psnr和ssim都是 Matlab 图像处理工具箱Image Processing Toolbox里的函数。psnr输入两张图必须是同尺寸同类型这里都是 double 且在 0-1 范围所以不需要额外转换。SSIM 是一个 0 到 1 之间的值比 PSNR 更能反映结构相似性。对于椒盐噪声去噪通常 AMF 的 PSNR 会高于普通中值滤波 3-8 dB具体数值受噪声密度和图像内容影响。3.5 参数调试r_max 选择与边界处理r_max窗口最大尺寸适应噪声密度计算耗时256×256主要风险13×3≤0.2快高密度噪声下残留大量椒盐点25×5≤0.4中等窗口跨边缘时细节区出现轻度糊化37×7≤0.6慢宽窗口下平坦区域可能出现块状伪纹理49×9≤0.7很慢边缘过冲明显细线结构丢失边界处理上padarray用symmetric已经规避了大部分越界问题但要留意r_max大于图像最小边长的一半时pad会大于实际图像尺寸这时镜像范围也会出现异常。正常实验图都在 256×256 以上r_max4 完全安全。提示如果输入图像不是方形padarray的pad参数是一个二维向量[pad_row, pad_col]行和列方向的扩展尺寸可以不一样。上面代码统一用了pad r_max对非方形图像误差也不大但严谨的做法是行方向用r_max、列方向也用r_max保持一致。4. 仿真实验噪声密度变化与去噪效果对比分析4.1 噪声密度从 0.1 到 0.7 的梯度测试单张图固定密度的对比只能说明一种情况下的表现。要系统地评价自适应中值滤波得让噪声密度在一个区间内变化观察三个指标PSNR、SSIM、主观视觉残留。下面脚本完成这个任务rng(0); % 固定随机种子 density_list 0.1:0.1:0.7; psnr_med zeros(length(density_list), 1); psnr_amf zeros(length(density_list), 1); for k 1:length(density_list) d density_list(k); In imnoise(I, salt pepper, d); I_med medfilt2(In, [3 3]); I_amf adaptive_median_filter_fast(In, 1, 3); psnr_med(k) psnr(I_med, I); psnr_amf(k) psnr(I_amf, I); end % 打印表格数据 T table(density_list, psnr_med, psnr_amf, ... VariableNames, {NoiseDensity, MedianPSNR, AMF_PSNR}); disp(T);rng(0)让每次循环生成的噪声位置一致这样对比的是两种算法在同一种噪声分布下的差异而不是随机误差。这一步在写论文或报告时尤其重要——可复现性是仿真实验的基本要求。当噪声密度达到 0.7 时3×3 固定窗口中几乎每个窗口都有 5 个左右的噪声像素排序中值大概率落在噪声灰度上medfilt2的结果会呈现明显的水渍状色斑而 AMF 因为可以扩到 7×7 或 9×9取到有效信号的概率就高一些。但注意r_max3 时 0.7 的密度已经逼近理论上限如果密度继续升高效果会快速下降。4.2 边缘像素在滤波前后的灰度剖面中值滤波最被诟病的点是边缘过度平滑而自适应版本在边缘处通常选择保留原始值。为了可视化这一点取图中一行像素画出原始图像、噪声图像、普通中值滤波、自适应滤波四种情况下的灰度剖面曲线。row 100; % 取第 100 行观察 figure; subplot(2,1,1); plot(I(row, :), k-, LineWidth, 1.2); hold on; plot(In(row, :), b., MarkerSize, 4); legend(原始, 带噪); title(第100行灰度剖面原始与噪声); ylim([0 1]); subplot(2,1,2); plot(I_med3(row, :), r--, LineWidth, 1.2); hold on; plot(I_amf3(row, :), g-, LineWidth, 1.2); legend(medfilt2, AMF); title(第100行灰度剖面滤波结果对比); ylim([0 1]);剖面图能非常直观地说明问题在明显的灰度阶跃物体边缘处medfilt2会把跳变沿拉平成一个斜坡而 AMF 的曲线会在阶跃附近保留更陡的过渡。对比原始图像曲线能看出 AMF 在平坦区域有明显去噪能力同时在边缘处更接近原始曲线。这个可视化手法常出现在课程报告和论文的结果分析部分也是答辩时评委爱看的内容。4.3 局部区域放大对比细节保留程度全局 PSNR 无法反映局部细节保留情况尤其是细线、纹理这类高频信息。建议用imcrop截取图像中一块包含细节的区域放大后对比。rect [80 60 60 60]; % 裁剪区域 [x, y, width, height] I_crop_orig imcrop(I, rect); I_crop_noise imcrop(In, rect); I_crop_med imcrop(I_med3, rect); I_crop_amf imcrop(I_amf3, rect); % 拼接成 2×2 montage 显示 figure; subplot(2,2,1); imshow(I_crop_orig); title(原始); subplot(2,2,2); imshow(I_crop_noise); title(带噪); subplot(2,2,3); imshow(I_crop_med); title(medfilt2); subplot(2,2,4); imshow(I_crop_amf); title(AMF);放大区域的视觉差异通常比全图更明显。medfilt2输出中能看到的细线断裂、屋顶纹理变模糊在 AMF 中会清晰不少。如果看不出来可以尝试截取 cameraman 图中的三脚架区域那个位置的密集细线是检验中值滤波类算法的最佳试金石。注意imcrop的坐标原点在图像左上角宽度高度按像素计算数据归一化不影响裁剪结果。4.4 运行时间统计性能这块在论文里是硬指标。用tic/toc记录运行时间注意多次运行取平均值排除抖动。funcs {() medfilt2(In, [3 3]), ... () adaptive_median_filter_fast(In, 1, 3), ... () adaptive_median_filter_fast(In, 1, 4)}; labels {medfilt2-3x3, AMF-rmax3, AMF-rmax4}; for k 1:length(funcs) t zeros(1, 5); for trial 1:5 f funcs{k}; tic; f(); t(trial) toc; end fprintf(%s 平均耗时: %.4f s\n, labels{k}, mean(t)); endmedfilt2底层是高效的 C 实现自适应版本因为是纯 Matlab 循环速度上被拉开差距是预期内的事。对 256×256 图像medfilt2耗时 0.01s 量级而 AMF 的纯循环实现通常要 2s 左右。如果仿真图是 1024×1024纯循环版本可能要到 30 秒以上。实际项目中要提速可以考虑把最内层循环用mex重构或者用parfor并行化——后面第 5 章会讲一个不需要离开 Matlab 环境的加速思路。5. 把仿真做成操作演示视频录屏与动图两种路线5.1 用 Matlab 自带的 VideoWriter 录制滤波过程标题里提到“含代码操作演示视频”说明这个项目的交付物不只是算法本身还包含可演示的视频材料。Matlab 官方推荐的做法是用VideoWriter把多帧图像写成一个 AVI 或 MP4 文件。最常见的演示视频思路是先展示原始图 → 展示加噪图 → 窗口逐渐扩大的动画 → 滤波结果图。下面给出一个录制窗口动画的代码。% 生成演示视频窗口逐帧扩大的滤波过程 v VideoWriter(AMF_demo.avi, Motion JPEG AVI); v.FrameRate 5; % 每秒 5 帧视觉上刚好能看清变化 open(v); % 先写 15 帧原始图停顿感靠多帧重复实现 for f 1:15 imshow(In); title(带噪声图像); frame getframe(gcf); writeVideo(v, frame); end % 窗口从 1 扩到 3每级窗口显示 10 帧 for r 1:3 I_tmp adaptive_median_filter_fast(In, 1, r); for f 1:10 imshow(I_tmp); title(sprintf(自适应中值滤波 (r %d), r)); frame getframe(gcf); writeVideo(v, frame); end end close(v);参数说明VideoWriter第一个参数是输出文件名扩展名决定封装格式。.avi和Motion JPEG AVI是兼容性最好的组合生成的视频体积较大但任何播放器都能打开。.mp4需要系统里有对应的编码器Windows 上通常能用Linux 上可能需要确认是否安装了 H.264 编码支持。getframe(gcf)捕获当前 figure 窗口的内容注意 gcf 窗口大小会直接影响视频分辨率建议在录制前用set(gcf, Position, [100 100 640 480])固定窗口尺寸避免不同帧分辨率抖动。视频里的“代码操作演示”部分如果不需要展示实时代码输入可以用live script.mlx 文件逐段运行代码块用系统录屏软件Windows 的 Xbox Game Bar 或 OBS录屏这比 Matlab 内部录制的画面更丰富。5.2 用 write 命令把关键中间结果输出为 GIF视频适合汇报答辩GIF 则更适合放进课程报告或博客附图中。Matlab 里把多帧图像写进 GIF 的方式是用imwrite的WriteMode, append选项。% 输出一组 GIF 动画原始 - 加噪 - 滤波结果 frames {I, In, I_amf3}; filename amf_steps.gif; for k 1:length(frames) [A, map] gray2ind(frames{k}, 256); if k 1 imwrite(A, map, filename, gif, LoopCount, Inf, DelayTime, 1); else imwrite(A, map, filename, gif, WriteMode, append, DelayTime, 1); end endgray2ind把 double 灰度图转成带 colormap 的索引图像这是 GIF 格式要求的。DelayTime单位是秒1 表示每帧停留 1 秒如果想让加噪过程更快改成 0.5 即可。LoopCount, Inf表示无限循环播放。这套写法在写博客、做实验报告时很实用生成的 GIF 可以直接插入 Markdown 或者 Word不需要任何视频剪辑软件。5.3 演示视频中如何展示代码与参数录屏演示的场景下很多同学会直接把脚本从头到尾运行一遍视频里全是代码滚动观看者根本抓不住重点。更好的顺序是先用一句中文注释交代实验目的例如% 比较固定中值与自适应中值在噪声密度0.5下的去噪效果。运行加噪代码图窗弹出停顿 2 秒让观众看到噪声图。运行medfilt2结果再运行 AMF 结果两张图并排对比。命令行再打印一组 PSNR 数值。最后运行VideoWriter片段把之前录的视频文件播一遍。配合录屏软件Matlab 编辑器里的代码高亮和命令行输出都能被记录下来效果比直接导出 PDF 脚本要生动得多。如果环境支持还可以在.mlx活页脚本里把代码、输出图和文字说明放在同一个文档流中每段代码块执行后立即在右侧显示生成的图窗。这样录屏时就只滚动活页脚本不需要反复切换窗口观感上更接近现代 Notebook 风格。5.4 验证演示视频正确性的自查清单做完视频不要直接交差花两分钟检查这几个点窗口尺寸是否真的在变化看 title 中的 r 值或者画面中噪声残留的变化、不同算法对比时用是否是同一张噪声图由于rng没有固定两次imnoise产生的噪声位置可能不同对比会失真、视频文件是否能被主流播放器解码。最容易翻车的其实是最后一个——Matlab 在 Linux 下生成的某些编码格式 Windows 播放器不支持保险起见统一用 Motion JPEG AVI。注意演示视频中的结果图应当与正文中的 PSNR 数值来源于同一组固定随机种子的运行结果。如果录视频时重新执行了脚本而没有设置rng噪声分布就可能与论文表格中的数据对不上这种不一致在答辩时大概率会被追问。本文还有配套的精品资源点击获取
返回列表