基于Opencv的自适应中值滤波函数selfAdaptiveMedianBlur（）

7.3.3 自适应滤波器

自适应中值滤波器

对于7.3.2节所讨论的中值滤波器，只要脉冲噪声的空间密度不大，性能还是可以的（根据经验需P_a和P_b小于0.2）。本节将证明，自适应中值滤波器可以处理更大概率的脉冲噪声。自适应中值滤波器的另一个优点是平滑非脉冲噪声时，试图保留细节，这是传统中值滤波器所做不到的。正如前面几节中所讨论的所有滤波器一样，自适应中值滤波器也工作于矩形窗口区域S_xy内。然而，与这些滤波器不同的是自适应中值滤波器在进行滤波处理时，会根据本节列举的某些条件而改变（增大或缩小）S_xy的尺寸。记住，滤波器的输出是一个单值，该值用于替代点(x,y)处的像素值，位于(x,y)处的点是S_xy的中心也是它的锚点。

考虑如下符号：

Zmin=S_xy中的最小灰度值；

Zmax=S_xy中的最大灰度值；

Zmed=S_xy中的中值；

Zxy=坐标(x,y)处的灰度值；

Smax=S_xy所允许的最大尺寸(在程序中，用kernal_size表示)；

自适应中值滤波算法以两个进程工作，分别为进程A和B，如下所示：

进程A：

如果Zmin<Zmed<Zmax,则转到进程B
否则增大窗口尺寸；
增大后的窗口尺寸（程序中用ks表示），
如果ks<Smax，则重复A。（在程序中，kernal_size即Smax）
否则输出Zmed

进程B：

如果Zmin<Zxy<Zmax，则输出Zxy
否则，输出Zmed.

该算法的设计意图是要实现3个目的：①去除椒盐噪声；②平滑其它非脉冲噪声；③减少物体边界细化或粗化的失真。Zmin和Zmax和在算法统计上，认为是类脉冲噪声分量，既使它们不是图像中的最大值或最小值。

进程A的目的是确定中值滤波器的输出Zmed，是否是一个脉冲（黑或白）。如果条件Zmin<Zmed<Zmax有效，则根据前节提到的原因，Zmed不可能是脉冲。这种情况下，转到进程B，检验窗口S_xy的中心点Zxy（即锚点）是否是一个脉冲。如果满足条件，就不是脉冲，原因与前同。这时，算法输出一个未修改的像素值Zxy。通过不修改这些“中间灰度级”的点，减少图像中的失真。如果Zmin<Zxy<Zmax为假，则Zxy=Zmin或Zxy=Zmax。在任何一种情况下，像素值都是一个极端值，且算法输出中值Zmed，从进程A可知Zmed不是脉冲噪声。最后一步是执行标准的中值滤波。问题是，标准中值滤波器使用图像中相应邻域的中值代替图像中的每一点，这会引起不必要的细节损失。

继续上面的说明，假设进程A确实找到了一个脉冲（若不是，则转到进程B），算法会增大窗口尺寸，并重复进程A。该循环会一直继续，直到算法找到一个非脉冲的中值，并转到进程B。如果循环中窗口达到了最大尺寸，则算法会返回值Zmed。注意，这并不能保证该值不是脉冲。噪声的概率Pa或Pb越小，或者在S_xy允许的范围内越大，退出条件也就越难满足。这是合理的，随着脉冲密度的增大，我们会需要更大的窗口来消除尖峰噪声。

算法没输出一个值，窗口S_xy就被移动到图像中的下一个位置。然后，算法重新初始化并应用到新位置的像素。仅使用新像素就可以反复更新中值，因而减少了计算开销。

#include <iostream>
#include <opencv2/opencv.hpp>
#include <vector>
#include <algorithm>
using namespace cv;
using namespace std;
//下面的宏，定义了在矩阵src的第m行、n列，ks*ks覆盖的矩形区域内的像素，并将像素压到矢量v中
//该覆盖区域的左上角坐标为（m，n），宽为ks，高为ks，要求src必须是单通道，数据类型为CV_8UC1
#define CV_ROI_ELEM(src,vector,m,n,ks)  
{                                      
    uchar* kn;                         
    int st0=src.step[0];
    int st1=src.step[1];
    for(int k=0;k<(ks);k++)            
    {                                  
        for(int s=0;s<(ks);s++)        
        {                              
            kn =src.data+(k+m)*st0+(s+n)*st1;   
            vector.push_back(*kn);              
        }                                       
    }                                           
}

#define CV_MAT_ELEM2(src,dtype,y,x) 
    (dtype*)(src.data+src.step[0]*(y)+src.step[1]*(x))
/*********************自适应中值滤波********************************/
void selfAdaptiveMedianBlur(Mat&src,Mat&dst,int kernal_size)
{
    CV_Assert(src.type()==CV_8UC1||src.type()==CV_8U);
    if(dst.empty())
    {
        dst.create(src.rows,src.cols,CV_8UC1);
    }
    uchar* pdst=dst.data;
    uchar Zmin,Zmax,Zmed,Zxy;
    int step0=src.step[0];
    int step1=src.step[1];
    for(int i=kernal_size/2;i<src.rows-kernal_size/2;i++)
    {
        for(int j=kernal_size/2;j<src.cols-kernal_size/2;j++)
        {
            int ks=3;//kernal_size;
            int count=0;
            Zxy=*CV_MAT_ELEM2(src,uchar,i,j);//Sxy覆盖区域的中心点像素值，即锚点像素值
            vector<uchar> v;//将模板覆盖区域的像素，压入矢量v中
            do{
                if(cout==0)
                {//获取模板ks*ks覆盖区域的像素，压入矢量v中
                    CV_ROI_ELEM(src,v,i-ks/2,j-ks/2,ks);
                }
               else
                {
               /****************下面的for循环，将外扩的四个边的像素添加到v中**************/
                    uchar* p=src.data+(i-ks/2)*step0+(j-ks/2)*step1;
                    for(int u=0;u<ks;u++)
                    {
                        v.push_back(*(p+u*step1));//向外扩展的四个边的上边
                        v.push_back(*(p+(ks-1)*step0+u*step1));//向外扩展的四个边的下边
                        if(u!=0&&u!=ks-1)
                        {
                          v.push_back( *(p+u*step0));//向外扩展的四个边的左边
                          v.push_back(*(p+u*step0+(ks-1)*step1));//向外扩展的四个边的右边
                        }
                    }
                }

                //对v的元素排序
                //排序后，Sxy覆盖区域内，最大值为Zmax=v[v.size-1],最小值为Zmin=v[0]
                std::sort(v.begin(),v.end());
                Zmin=v[0],Zmax=v[v.size()-1],Zmed=v[ks*ks/2];
                pdst =CV_MAT_ELEM2(dst,uchar,i,j);
                if(Zmin<Zmed&&Zmed<Zmax)
                {
                    if(Zmin<Zxy&&Zxy<Zmax)
                        {*pdst=Zxy;break;}
                    else
                        {*pdst=Zmed;break;}
                }
                else
                {
                    ks +=2;
                }
                count++;
            }while(ks<=kernal_size);

            *pdst=Zmed;
        }
    }
}


int main()
{
    Mat src=imread("D:\Qt\MyImage\3.bmp",0);
    imshow("src image",src);

    Mat dst;
    selfAdaptiveMedianBlur(src,dst,7);

    imshow("adaptive median filter",dst);

    waitKey();

    return 0;
}

原图像：

自适应中值滤波后的图像：