NAFSM中值滤波器讲解与实现

去年大三时课上，医学图像处理老师讲解了这个算法，当时布置的作业就是实现这个中值滤波器，今天突然从一篇论文中看到过类似这种思想，就又想起来了，发现网上对这个算法的讲解很少，那我就在这分享一下咯！下面是原论文：
K. K. V. Toh and N. A. Mat Isa, “Noise Adaptive Fuzzy Switching Median Filter for Salt-and-Pepper Noise Reduction,” in IEEE Signal Processing Letters, vol. 17, no. 3, pp. 281-284, March 2010.，链接：https://ieeexplore.ieee.org/document/5356178

接下来我将从我对该算法理解、算法实现、算法MATLAB代码、实验报告这个顺序进行我的讲解，希望能帮助到一些人。
一、NAFSM算法理解：
首先发一张我的总结：

第一步：是一个二值化处理，对于图像像素灰度值最大和最小当做椒盐噪声，该像素值变成0。
注意图中公式N（i，j）和X（i，j）的区别，这对于理解后面推导很重要，其中N代表二值化后的矩阵，值为1的像素代表该点不是噪声，为0的像素代表该点为噪声点，而对于X则还是原图的灰度值，后面M代表准恢复图像像素灰度值，Y代表最后处理后的图像像素灰度值大小。
第二步：定义一个矩形框，图中s即代表矩形框的大小参数，如果s=1，代表矩形框大小为33；如果s=2，代表矩形框大小为55，以此类推；
第三步：因为对于处理这类噪声，主要难点就是噪声点灰度大小恢复的问题，作者考虑当在该像素点矩形框内，有存在不是噪声点的点，也就是图片中G的大小不为0，则直接可以提取矩形框内那些不为0的点（非噪声像素点）灰度值的中值作为该点灰度值；
第四步：然而，如果发现矩形框内都是噪声点，也就是矩形框内的N（i，j）都为0，这个时候就需要将矩形框扩大，比如s=1，变成s=2，继续第二步，如果最后发现s=3，也就是77矩形框内结果还是都为噪声点，即G=0，则直接将该点灰度值定义为其33框内，左上角四个像素灰度值的中值（这里解释一下，由于一张图处理一般都是从左上角像素开始，所以这里就假设了正在处理的像素点之前的像素点都已经获得比较靠谱的灰度值大小，所以这里就只选取图中几个像素点灰度大小的中值）；
第五步：经过前面几步，至此，M的大小已全部确定，但是这里作者没有直接取M值作为最后处理结果，作者对原始灰度图像计算33矩形框内与中心点灰度值绝对值相差最大值；
第六步：通过一个函数计算F值大小，可参考下图理解：

结合图中那个求F公式可以知道，其主要思想是：当原始图像33矩形框内有像素灰度差异比较大时，认为这个点灰度可信度比较低，也即是F为1，对应于Y值就是不考虑X（原始灰度）大小，当矩形框内像素灰度值差异很小时，认为这点可信度很高（因为如果框内存在椒盐噪声，则必为0或者225，大概率两者都有，所以这个时候值只会很大，不会很小），这个时候F就为0，也就是相信X值大小，同理可理解中间那段（如上图中间那段呈线性变化）。我在matlab程序里选取的T1和T2值分别为30和60，F值是用来衡量M与真实值之间误差；
第七步：最后通过公式求出Y，即最后结果。

二、算法实现
这是matlab上运行部分结果图：
加10%椒盐噪声：

加20%椒盐噪声

加50%椒盐噪声：
加70%椒盐噪声：
加90%椒盐噪声：

可以看到，NAFSM滤波器效果非常好，碾压普通中值滤波器。

三、MATLAB代码：
1、原始代码

clear
[Im,map]=imread('lena.gif');
J=ind2gray(Im, map);
I= imnoise(J,'salt & pepper',0.5);%加入椒盐噪声
[m n]=size(I);
N=zeros(m,n);%定义一个二值图空矩阵
for i=1:1:m %确定图像二值图，将灰度值为0或者255的像素定义为椒盐噪声，二值图大小为0for j=1:1:nif I(i,j)==255||I(i,j)==0N(i,j)=0;elseN(i,j)=1;endend
end
s=input('请输入大小为s*s的矩形窗，s=');
G=zeros(m,n);%定义一用来保存矩形窗内二值总大小的矩阵
for i=1:1:m%计算矩形窗内各像素大小，并且求预恢复值M（i，j） for j=1:1:nW=zeros(s,s);%定义一个空的矩形窗for li=1:1:sfor lj=1:1:sif (i+li-(s+1)/2)<=0||(j+lj-(s+1)/2)<=0||(i+li-(s+1)/2)>m||(j+lj-(s+1)/2)>n%判断是否超出原始图像范围continue;endW(li,lj)=N(i+li-(s+1)/2,j+lj-(s+1)/2);%获取矩形窗内数据endendG(i,j)=sum(sum(W));z=s;for k=s:2:5if G(i,j)==0z=k+2;W1=zeros(z,z);%定义一个空的矩形窗for li=1:1:zfor lj=1:1:zif (i+li-(z+1)/2)<=0||(j+lj-(z+1)/2)<=0||(i+li-(z+1)/2)>m||(j+lj-(z+1)/2)>n%判断是否超出原始图像范围continue;endW1(li,lj)=N(i+li-(z+1)/2,j+lj-(z+1)/2);%获取矩形窗内数据endendG(i,j)=sum(sum(W1));if G(i,j)==0continue;elsebreak;endelsebreak;endendif z==7M(i,j)=median(I(i-1,j-1),I(i,j-1),I(i+1,j-1),I(i-1,j));
%            M(i,j)=I(i,j);continue;endW2=zeros(1,G(i,j));nn=1;for i1=1:1:zfor j1=1:1:zif (i+i1-(z+1)/2)<=0||(j+j1-(z+1)/2)<=0||(i+i1-(z+1)/2)>m||(j+j1-(z+1)/2)>n%判断是否超出原始图像范围continue;      elseif N(i+i1-(z+1)/2,j+j1-(z+1)/2)==1W2(nn)=I(i+i1-(z+1)/2,j+j1-(z+1)/2);nn=nn+1;endendendM(i,j)=median(W2);end
end
W3=zeros(3,3);
D=zeros(m,n);
for i=1:1:mfor j=1:1:nL=zeros(3,3);%定义一个m*n的矩形窗for li=1:1:3for lj=1:1:3if (i+li-(3+1)/2)<=0||(j+lj-(3+1)/2)<=0||(i+li-(3+1)/2)>m||(j+lj-(3+1)/2)>n%判断是否超出原始图像范围continue;endL(li,lj)=I(i+li-(3+1)/2,j+lj-(3+1)/2);%获取矩形窗内数据endendd=zeros(1,9);ii=1;for i1=1:1:3for j1=1:1:3d(ii)=abs(L(i1,j1)-L(2,2));ii=ii+1;endendD(i,j)=max(d);end
end
F=zeros(m,n);
T1=20;
T2=60;
for i=1:1:mfor j=1:1:nif D(i,j)<T1F(i,j)=0;elseif D(i,j)>=T1&&D(i,j)<T2F(i,j)=(D(i,j)-T1)/(T2-T1);elseF(i,j)=1;endend
end
Y=zeros(m,n);
for i=1:1:mfor j=1:1:nY(i,j)=(1-F(i,j))*I(i,j)+F(i,j)*M(i,j);end
end
figure(1)
subplot(2,2,1)
imshow(J);
title('原始图像');
subplot(2,2,2)
imshow(I);
title('加入椒盐噪声后的图像');
subplot(2,2,3)
T=Median_filter(I,s,s);
T=uint8(T);
imshow(T);
title('普通中值滤波函数滤波后图像')
subplot(2,2,4)
Y=uint8(Y);
imshow(Y);
title('NAFSM滤波器滤波后的图像');
function [T]=Median_filter(I,m,n)
[M,N]=size(I);
T=zeros(M,N);
for i=1:1:Mfor j=1:1:NL=zeros(m,n);%定义一个m*n的矩形窗for li=1:1:mfor lj=1:1:nif (i+li-(m+1)/2)<=0||(j+lj-(n+1)/2)<=0||(i+li-(m+1)/2)>M||(j+lj-(n+1)/2)>N%判断是否超出原始图像范围continue;endL(li,lj)=I(i+li-(m+1)/2,j+lj-(n+1)/2);%获取矩形窗内数据endendT(i,j)=median(median(L));%求矩形窗内中值end
end
end

经过qq_50872053提醒，之前代码有点不完美，具体：
（1）、步骤4中，G=0时，应该取的是每次都更新后的值，但是这里一直取得是原始未处理图像数据
（2）、之前没有考虑i,j循环顺序，即考虑先从图片横轴开始处理还是纵轴开始处理
新的修改后代码如下（实际效果有进一步提升）：

2、最新版代码

clear
[Im,map]=imread('lena.gif');
J=ind2gray(Im, map);
J_noise= imnoise(J,'salt & pepper',0.90);%加入椒盐噪声
I = J_noise;
[m n]=size(I);
N=zeros(m,n);%定义一个二值图空矩阵
for i=1:1:m %确定图像二值图，将灰度值为0或者255的像素定义为椒盐噪声，二值图大小为0for j=1:1:nif I(i,j)==255||I(i,j)==0N(i,j)=0;elseN(i,j)=1;endend
end
s=input('请输入大小为s*s的矩形窗，s=');
G=zeros(m,n);%定义一用来保存矩形窗内二值总大小的矩阵
for j=1:1:n%计算矩形窗内各像素大小，并且求预恢复值M（i，j） for i=1:1:mW=zeros(s,s);%定义一个空的矩形窗for li=1:1:sfor lj=1:1:sif (i+li-(s+1)/2)<=0||(j+lj-(s+1)/2)<=0||(i+li-(s+1)/2)>m||(j+lj-(s+1)/2)>n%判断是否超出原始图像范围continue;endW(li,lj)=N(i+li-(s+1)/2,j+lj-(s+1)/2);%获取矩形窗内数据endendG(i,j)=sum(sum(W));z=s;for k=s:2:5if G(i,j)==0z=k+2;W1=zeros(z,z);%定义一个空的矩形窗for li=1:1:zfor lj=1:1:zif (i+li-(z+1)/2)<=0||(j+lj-(z+1)/2)<=0||(i+li-(z+1)/2)>m||(j+lj-(z+1)/2)>n%判断是否超出原始图像范围continue;endW1(li,lj)=N(i+li-(z+1)/2,j+lj-(z+1)/2);%获取矩形窗内数据endendG(i,j)=sum(sum(W1));if G(i,j)==0continue;elsebreak;endelsebreak;endendif z==7 && G(i,j)==0if (i-1)<=0||(j-1)<=0||(i+1)>m||(j+1)>n %判断是否超出原始图像范围continue;  elseM(i,j)=median([I(i-1,j-1),I(i,j-1),I(i+1,j-1),I(i-1,j)]);endendW2=zeros(1,G(i,j));nn=1;for i1=1:1:zfor j1=1:1:zif (i+i1-(z+1)/2)<=0||(j+j1-(z+1)/2)<=0||(i+i1-(z+1)/2)>m||(j+j1-(z+1)/2)>n%判断是否超出原始图像范围continue;      elseif N(i+i1-(z+1)/2,j+j1-(z+1)/2)==1W2(nn)=I(i+i1-(z+1)/2,j+j1-(z+1)/2);nn=nn+1;endendendM(i,j)=median(W2);W3=zeros(3,3);
D=zeros(m,n);
L=zeros(3,3);%定义一个m*n的矩形窗
for li=1:1:3for lj=1:1:3if (i+li-(3+1)/2)<=0||(j+lj-(3+1)/2)<=0||(i+li-(3+1)/2)>m||(j+lj-(3+1)/2)>n%判断是否超出原始图像范围continue;endL(li,lj)=I(i+li-(3+1)/2,j+lj-(3+1)/2);%获取矩形窗内数据end
end
d=zeros(1,9);
ii=1;
for i1=1:1:3for j1=1:1:3d(ii)=abs(L(i1,j1)-L(2,2));ii=ii+1;end
end
D(i,j)=max(d);F=zeros(m,n);
T1=20;
T2=60;
if D(i,j)<T1F(i,j)=0;
elseif D(i,j)>=T1&&D(i,j)<T2F(i,j)=(D(i,j)-T1)/(T2-T1);
elseF(i,j)=1;
end% Y=zeros(m,n);
% Y(i,j)=(1-F(i,j))*I(i,j)+F(i,j)*M(i,j);
I(i,j) = (1-F(i,j))*I(i,j)+F(i,j)*M(i,j);end
end        figure(1)
subplot(2,2,1)
imshow(J);
title('原始图像');
subplot(2,2,2)
imshow(J_noise);
title('加入椒盐噪声后的图像');
subplot(2,2,3)
T=Median_filter(J_noise,s,s);
T=uint8(T);
imshow(T);
title('普通中值滤波函数滤波后图像')
subplot(2,2,4)
Y=uint8(I);
imshow(I);
title('NAFSM滤波器滤波后的图像');
function [T]=Median_filter(I,m,n)
[M,N]=size(I);
T=zeros(M,N);
for i=1:1:Mfor j=1:1:NL=zeros(m,n);%定义一个m*n的矩形窗for li=1:1:mfor lj=1:1:nif (i+li-(m+1)/2)<=0||(j+lj-(n+1)/2)<=0||(i+li-(m+1)/2)>M||(j+lj-(n+1)/2)>N%判断是否超出原始图像范围continue;endL(li,lj)=I(i+li-(m+1)/2,j+lj-(n+1)/2);%获取矩形窗内数据endendT(i,j)=median(median(L));%求矩形窗内中值end
end
end

加90%椒盐噪声（修改后）：

可以发现，明显比之前效果要更好

对于本次实验报告、图片和源matlab代码，可以去我主页下载附件！！！！！！！！！！！
转载请注明出处，谢谢

NAFSM中值滤波器讲解与实现相关推荐

【STM32F407的DSP教程】第48章 STM32F407的中值滤波器实现，适合噪声和脉冲过滤（支持逐个数据的实时滤波）
完整版教程下载地址:http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547 第48章 STM32F407的中值滤波器实现,适 ...
【STM32F429的DSP教程】第48章 STM32F429的中值滤波器实现，适合噪声和脉冲过滤（支持逐个数据的实时滤波）
完整版教程下载地址:http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547 第48章 STM32F429的中值滤波器实现,适 ...
从命令行到IDE，版本管理工具Git详解（远程仓库创建+命令行讲解+IDEA集成使用）
首先,Git已经并不只是GitHub,而是所有基于Git的平台,只要在你的电脑上面下载了Git,你就可以通过Git去管理"基于Git的平台"上的代码,常用的平台有GitHub.Gi ...
详细通俗重点CRF层讲解
本文翻译自GitHub博客上的原创文章,结尾有原文链接.文章没有晦涩的数学公式,而是通过实例一步一步讲解CRF的实现过程,是入门CRF非常非常合适的资料. 相关项目代码: BERT-BiLSMT-CR ...
高级数据结构讲解与案例分析
然而,仅仅掌握好它们不足以应付大厂的算法面试的.为了达到对时间和空间复杂度的理想要求,本节课探究高级数据结构,它们的实现要比那些常用的数据结构要复杂得多.其中重点介绍: 优先队列图前缀树线段树 ...
php 伪静态 page-18.html,PHP 伪静态实现技术原理讲解
PHP 伪静态实现技术原理讲解发布于 2015-01-18 23:52:58 | 129 次阅读 | 评论: 0 | 来源: 网友投递 PHP开源脚本语言PHP(外文名: Hypertext Pre ...
ssm开发框架原理_SSM 单体框架 - 前端开发：视频讲解
视频讲解知乎视频www.zhihu.com 知乎视频www.zhihu.com 知乎视频www.zhihu.com 知乎视频www.zhihu.com 知乎视频www.zhihu.com ...
python计算wav的语谱图_Python实现电脑录音（含音频基础知识讲解）
前言今天开始进入近期系列文章的第一篇,如何用 Python 来实现录音功能. 在开始"造轮子"之前,个人一直强调一个观点,如果有些东西已经有了,不妨直接去 github 上搜,用 ...
PCL：k-d tree 1 讲解
1.简介 kd-tree简称k维树,是一种空间划分的数据结构.常被用于高维空间中的搜索,比如范围搜索和最近邻搜索.kd-tree是二进制空间划分树的一种特殊情况.(在激光雷达SLAM中,一般使用的是三 ...
基础矩阵，本质矩阵，单应性矩阵讲解
ORB-SLAM点云地图中相机的位姿初始化,无论算法工作在平面场景,还是非平面场景下,都能够完成初始化的工作.其中主要是使用了适用于平面场景的单应性矩阵H和适用于非平面场景的基础矩阵F,程序中通过一个 ...

NAFSM中值滤波器讲解与实现

NAFSM中值滤波器讲解与实现相关推荐

最新文章

热门文章