图像线性平滑空间滤波（加权均值滤波器，几何均值滤波，谐波均值滤波，逆谐波均值滤波），非线性平滑空间滤波（中值滤波，最大值滤波，最小值滤波）MATLAB自写函数实现

平滑空间域滤波

1 理论

1.1 空间域滤波

滤波是信号处理中的一个概念，是将信号中特定波段频率滤除的操作，在数字信号处理中通常通过傅里叶变换及其逆变换实现。由于空间域滤波实际上和通过傅里叶变换实现的频域下的滤波是等效的，故而也称为滤波。不同的是空间域滤波主要直接基于邻域（空间域）对图像中像素执行计算。

1.1.1 空间域滤波和邻域处理

对图像中的每一点$(x,y)$，重复下面的操作。

（1）对预先定义的以$（x,y）$为中心的邻域内的像素进行运算。

（2）将（1）中运算的结果作为$（x,y）$点新的响应。

上述过程就称为邻域处理或空间域滤波。一幅数字图像可以看成一个二维函数$f(x,y)$，而$x-y$平面表明了空间位置信息，称为空间域，基于$x-y$空间邻域的滤波操作叫作空间域滤波。如果对于邻域中的像素计算为线性运算，则又称为线性空间域滤波，否则称为非线性空间域滤波。

下图直观地展示了用一个3× 3的模板（又称为滤波器、模板、掩模、核或者窗口）进行空间滤波的过程，模板为$w$，用黑笔圈出的$w(0,0)$是其中心。

滤波过程就是在图像$f(x,y)$中逐点地移动模板，使模板中心和点$(x,y)$重合。在每一点$（x,y）$处，滤波器在该点的响应是根据模板的具体内容并通过预先定义的关系来计算。一般来说模板中的非0元素指出了邻域处理的范围，只有那些当模板中心与点$（x,y）$重合时，图像$f$中和模板中非0像素重合的像素参与了决定点$（x,y）$像素值的操作。

在线性空间滤波中模板的系数则给出了一种加权模式，即$（x,y）$处的响应由模板系数与模板下面区域的相应$f$的像素值的乘积之和给出。例如，在上图中，对于模板的响应$R$为：
$$R=w(-1,-1)f(x-1,y-1)+w(-1,0)f(x-1,y)+\dots +w(0,0)f(x,y)+\dots +w(1,0)f(x+1,y)+w(1,1)f(x+1,y+1)$$

更一般的情况，对于一个大小为$m\times n$的模板，其中$m=2a+1,n=2b+1$, $a,b$均为正整数，即模板长与宽均为基数，且可能的最小尺寸为3× 3（偶数尺寸的模板由于其不具有对称性因而很少被使用，而1× 1大小的模板的操作不考虑邻域信息，退化为图像点运算），可以将滤波操作形式化地表示为：
$$g(x,y)=\sum _{s=-a}^a\sum _{t=-b}^{b}w(s,t)f(x+s,y+t)$$

对于大小为$M\times N$的图像$f(0,\dots ,M-1,0,\dots ,N-1)$，对$x=0,1,2,\dots ,M-1$和$y=0,1,2,\dots ,N-1$依次应用公式，从而完成了对于图像$f$所有像素的处理，得到新的图像$g$。

1.1.2 边界处理

执行滤波操作要注意的一点是当模板位于图像边缘时，会产生模板的某些元素很可能会位于图像之外的情况，这时，对于在边缘附近执行滤波操作需要单独处理，以避免引用到本不属于图像的无意义的值。

以下3种策略都可以用来解决边界问题。

（1）收缩处理范围——处理时忽略位于图像$f$边界附近会引起问题的那些点。
设原图像为：
$$\left[\begin{array}{lllll}1& 1& 1& 1& 1\\ 2& 2& 2& 2& 2\\ 3& 3& 3& 3& 3\\ 4& 4& 4& 4& 4\end{array}\right]$$
对于原图像中所使用的模板，处理时忽略图像四周一圈1个像素宽的边界：
$$\left[\begin{array}{ccccc}-& -& -& -& -\\ -& 2& 2& 2& -\\ -& 3& 3& 3& -\\ -& -& -& -& -\end{array}\right]$$

" -"表示无法进行模板操作的像素点。

从而确保了滤波过程中模板始终不会超出图像的边界。

2）使用常数填充图像——根据模板形状为图像$f$虚拟出边界，虚拟边界像素值为指定的常数，如0，得到虚拟图像$f^{\prime }$，保证模板在移动过程中始终不会超出$f^{\prime }$的边界。

（3）使用复制像素的方法填充图像——和（2）基本相同，只是用来填充虚拟边界像素值的不是固定的常数，而是复制图像$f$本身边界的模式。

1.1.3 互相关和卷积

上述的滤波操作实际是信号处理中的互相关操作，互相关和卷积的略有不同，卷积表示如下。
$$g(x,y)=\sum _{s=-a}^a\sum _{t=-b}^{b}w(-s,-t)f(x+s,y+t)$$

两者尽管差别细微，但有本质不同，卷积在卷积时模板是相对其中心点做镜像后再对$f$位于模板下的子图像做加权和的，或者说在做加权和之前模板先要以其中心点为原点旋转180°。如果忽略了这一细微差别将导致完全错误的结果，只有当模板本身是关于中心点对称时，互相关和卷积的结果才会相同。

1.2 线性空间滤波

1.2.1 加权均值滤波

$$g(x,y)=\frac{{\sum }_{s=-a}^a{\sum }_{s=-b}^{b}w(s,t)f(x+s,y+t)}{{\sum }_{s=-a}^a{\sum }_{s=-b}^{b}w(s,t)}$$

模板大小为$m\times n$($m,n$为奇数)，则$a=\frac{m-1}{2},a=\frac{n-1}{2}$。$f(x,y)$为原图像在$(x,y)$的像素点；$g(x,y)$为滤波后得到的在$(x,y)$的像素点。

• 均值滤波模板示例

3x3	5x5
$\frac{1}{9}\left[\begin{array}{lll}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$	$\frac{1}{25}\left[\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1\end{array}\right]$

• 高斯滤波模板示例

3x3	5x5
$\frac{1}{16}\left[\begin{array}{lll}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]$	$\frac{1}{159}\times \left[\begin{array}{ccccc}2& 4& 5& 4& 2\\ 4& 9& 12& 9& 4\\ 5& 12& 15& 12& 5\\ 4& 9& 12& 9& 4\\ 2& 4& 5& 4& 2\end{array}\right]$

1.2.2 均值滤波的变形

1.2.2.1 几何均值滤波

$$g(x,y)=[\prod _{(s,t)\in S_{xy}}f(s,t){]}^{\frac{1}{m\cdot n}}$$

模板大小为$m\times n$，$f(x,y)$为原图像在$(x,y)$的像素点；$g(x,y)$为滤波后得到的在$(x,y)$的像素点。$S_{x,y}$为以$(x,y)$为中心的含有奇数个像素点的邻域窗口；

和算术均值滤波器相比，几何均值滤波器能够更好的去除高斯噪声，并且能够更多的保留图像的边缘信息。但其对0值是非常敏感的，在滤波器的窗口内只要有一个像素的灰度值为0，就会造成滤波器的输出结果为0。

1.2.2.2 谐波均值滤波

$$g(x,y)=\frac{mn}{{\sum }_{s,t\in S_{x,y}}\frac{1}{f(x,y)}}$$

模板大小为$m\times n$，$f(x,y)$为原图像在$(x,y)$的像素点；$g(x,y)$为滤波后得到的在$(x,y)$的像素点。$S_{x,y}$为以$(x,y)$为中心的含有奇数个像素点的邻域窗口；

谐波均值滤波器对盐粒噪声（白噪声）效果较好，不适用于胡椒噪声；也比较适合处理高斯噪声。

1.2.2.3 逆谐波均值滤波

$$g(x,y)=\frac{{\sum }_{s,t\in S_{x,y}}f(s,t{)}^{Q+1}}{{\sum }_{s,t\in S_{x,y}}f(s,t{)}^Q}$$

其中$Q$称为滤波器的阶数。当$Q$是正数时，滤波器用于消除“椒”噪声;当$Q$是负数时，滤波器用于消除“盐”噪声。但它不能同时消除这两种噪声。
当$Q=0$，逆谐波均值滤波器转变为算术均值滤波器;当$Q=-1$，逆谐波均值滤波器转变为谐波均值滤波器。

1.3 非线性空间滤波

1.3.1 中值滤波

g ( x , y ) = M e d ( s , t ) ∈ S x , y { f ( s , t ) } g(x,y)=Med_{(s,t)\in S_{x,y}}\{f(s,t)\} g(x,y)=Med(s,t)∈Sx,y​​{f(s,t)}

式中$f(x,y)$为原图像在$(x,y)$的像素点；$S_{xy}$为以像素$f(x,y)$为中心的含有奇数个像素点的邻域窗口；$g(x,y)$为滤波后得到的在$(x,y)$的像素点,$Med$为取中值。

适于处理椒盐噪声，通过多次使用小模板，可以获得很好的去噪效果，但多次应用中值滤波器,会使图像模糊。

1.3.2 最大值滤波

g ( x , y ) = max ⁡ ( s , t ) ∈ S x , y { f ( s , t ) } g(x,y)=\max_{(s,t)\in S_{x,y}}\{f(s,t)\} g(x,y)=(s,t)∈Sx,y​max​{f(s,t)}

式中$f(x,y)$为原图像在$(x,y)$的像素点；$S_{xy}$为以像素$f(x,y)$为中心的含有奇数个像素点的邻域窗口；$g(x,y)$为滤波后得到的在$(x,y)$的像素点,$max$为取最大值。

最大值滤波器对于 “椒” 噪声具有较好的消除效果，对于 “盐” 噪声则效果不好。

1.3.3最小值滤波

g ( x , y ) = min ⁡ ( s , t ) ∈ S x , y { f ( s , t ) } g(x,y)=\min_{(s,t)\in S_{x,y}}\{f(s,t)\} g(x,y)=(s,t)∈Sx,y​min​{f(s,t)}

式中$f(x,y)$为原图像在$(x,y)$的像素点；$S_{xy}$为以像素$f(x,y)$为中心的含有奇数个像素点的邻域窗口；$g(x,y)$为滤波后得到的在$(x,y)$的像素点,$min$为取最小值。

最小值滤波器对于 “盐” 噪声具有较好的消除效果，对于 “椒” 噪声则效果不好。

2 上代码！

2.1添加噪声

本文只使用椒盐噪声和高斯噪声，更多噪声添加可参见：
图像噪声模型

function imgOut = addNoise(imgIn,type,x,y)
% 函数功能：为图像添加噪声
% 关注微信公众号【二进制人工智能】，回复LB获取demo全部代码
% 输入：
%      imgIn:输入图形矩阵
%      type:噪声类型
%           高斯噪声：gaussian,参数为(x,y)，默认值为(0,10)
%           椒盐噪声： salt & pepper,强度为x，默认值为0.02
% 输出：
%      imgOut：处理后的图像
if ndims(imgIn)>=3
    imgIn=rgb2gray(imgIn);
end
[M,N]=size(imgIn);  % 输入图像的大小
% 设置默认噪声类型
if nargin==1
    type='gaussian';
end
switch lower(type)
    case 'gaussian'  
        if nargin<4
            y=10;        % 标准差默认值
        end
        if nargin<3
            x=0;          % 均值默认值
        end
        R=normrnd(x,y,M,N);        % 产生高斯分布随机数
        imgOut=double(imgIn)+R;            % 添加噪声
        imgOut=uint8(round(imgOut));
    case 'salt & pepper'
        if nargin<3
            x=0.02;
        end
        tem1=rand(M,N)<x;   % 椒噪声
        tem2=rand(M,N)<x;   % 盐噪声
        imgOut=imgIn;
        imgOut(tem1)=0;    
        imgOut(tem2)=255;
    otherwise
        error('Unkown distribution type')    
end
end

2.2 图像边界处理

function imgOut = padding(imgIn,fsize,type)
% 函数功能：边界处理函数
% 关注微信公众号【二进制人工智能】，回复LB获取demo全部代码
% 输入：
%      imgIn:输入图片
%      fsize：模板大小
%      type:padding类型
% 输出：
%      imgOut：处理后的图像
for i=1:(fsize-1)/2  % 模板为3，填充一层 模板为5填充两层，以此类推
    [m,n]=size(imgIn);
    imgTem=zeros(m+2,n+2);
    switch(lower(type))
        case 'zeros'        
            imgTem(2:m+1,2:n+1)=imgIn;
        case 'borden_value'  
            imgTem(2:m+1,2:n+1)=imgIn;
            % 边沿
            imgTem(1,2:n+1)=imgIn(1,:); imgTem(m+2,2:n+1)=imgIn(m,:); 
            imgTem(2:m+1,1)=imgIn(:,1); imgTem(2:m+1,n+2)=imgIn(:,n);
            % 边角
            imgTem(1,1)=imgIn(1,1); imgTem(1,n+2)=imgIn(1,n);
            imgTem(m+2,1)=imgIn(m,1);imgTem(m+2,n+2)=imgIn(m,n);
        otherwise
            disp('please input ''zeros''/''borden_value''') 
            break;
    end
    imgIn=imgTem;
end
imgOut=imgIn; 
end

2.3 生成高斯模板

function w = wGaussian(size,sigma)
% 函数功能：生成高斯模板
% 关注微信公众号【二进制人工智能】，回复LB获取demo全部代码
% 输入：
%     size:模板大小
%     sigma:高斯分布标准差
% 输出：
%      w:高斯模板
    k=(size-1)/2;
    [x,y]=meshgrid([-k:k],[-k:k]);
    w=exp(-(x.*x+y.*y)./(2*sigma^2))/(2*pi*sigma^2);
    w=w/sum(w(:));
end

2.4 线性滤波函数

function imgOut = linearFilter(imgIn,w,paddingType,filterType)
% 函数功能：线性滤波
% 关注微信公众号【二进制人工智能】，回复LB获取demo全部代码
% 输入：
%      imgIn：输入图像
%      w：模板或模板大小
%      paddingType：边缘填充类型
%      filterType：滤波类型
% 输出：
%      imgOut：处理后的图像
[m,n]=size(imgIn); 
if length(w(1,:))>1  % 如果w是模板
    wSize=size(w,1);
else                 % w是模板大小
    wSize=w;
end                    
imgPad=padding(imgIn,wSize,paddingType); 
imgOut=imgIn;
switch(filterType)
    case 'weighted_average'    % 加权均值滤波       
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1);  
                %中心点的值用子图像的算术均值代替
                imgOut(i,j) = sum(sum(childImg.*w))/sum(sum(w));                
            end
        end
    case 'Geometric_average'   % 几何均值滤波
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1);  
                %中心点的值用子图像的几何均值代替
                imgOut(i,j)=(prod(childImg(:)))^(1/wSize^2); 
            end
        end
    case 'Harmonic_average'  % 谐波均值滤波
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1);  
                %中心点的值用子图像的谐波均值代替
                imgOut(i,j)=wSize^2/(sum(sum(1./childImg))); 
            end
        end
    case 'Inverse_harmonic' % 逆谐波均值滤波
        Q=input('input the order：');
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1); 
                %中心点的值用子图像的逆谐波均值代替
                imgOut(i,j)=sum(sum(childImg.^(Q+1)))/sum(sum(childImg.^(Q)));
            end 
        end
    otherwise
        error('The filter type is invalid ')
end
end

2.5 非线性滤波函数

function imgOut = nonlinearFilter(imgIn,wSize,paddingType,filterType)
% 函数功能：非线性滤波
% 关注微信公众号【二进制人工智能】，回复LB获取demo全部代码
% 输入：
%      imgIn：输入图像
%      wSize: 模板大小
%      paddingType：边缘填充类型
%      filterType：滤波类型
% 输出：
%      imgOut：处理后的图像
[m,n]=size(imgIn);                    
imgPad=padding(imgIn,wSize,paddingType); 
imgOut=imgIn;
switch(filterType)
    case 'median'    % 中值滤波
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1);
                 %中心点的值用子图像的中值代替
                imgOut(i,j)=median(childImg(:));
            end
        end
    case 'max'       % 最大值滤波
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1);  
                %中心点的值用子图像的最大值代替
                imgOut(i,j)=max(childImg(:));
            end
        end
    case 'min'      % 最小值滤波
        for i = 1:m
            for j = 1:n
                childImg = imgPad(i:i+wSize-1,j:j+wSize-1);  
                %中心点的值用子图像的最小值代替
                imgOut(i,j)=min(childImg(:));
            end
        end
    otherwise
        error('This filter type is invalid ')
end
end

3 demo结果

3.1 高斯噪声

3.2 椒盐噪声

参考：

[1]https://blog.csdn.net/shushi6969/article/details/102992931

[2]https://blog.csdn.net/weixin_44378835/article/details/109005041

[3]https://zhuanlan.zhihu.com/p/187517189

关注微信公众号【二进制人工智能】，回复【LB】，获取demo所有代码。