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算法理解：
首先发一张我的总结：

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