中值、均值、高斯、双边滤波（第五天）

--------滤波简介--------

滤波的是图像处理之中必备的手段，也是必经之路。模式识别，深度学习都会用到滤波的相关内容，没有好的图片你怎么识别？

个人理解的“滤波”二字：

     滤波没有什么神奇的地方，就是通过旁边像素对比来做判断（其他的下面再说）举个例子：

 

 

                                    假如上面的图形1是瑕疵，通过滤波给过滤掉

上面是一个两张简单的图片，两个问题： 1.把左边的'形状1'清除(过滤掉)。

                                                  2.把右边的'形状1'清除(过滤掉)。

对于第一个问题：那很简单了，前面的博文已经知道，图像是一个个像素点构成的，那对比蓝色和白色的像素点就可以知道，直接把蓝色排除了。

对于第二个问题：那就复杂了，因为图形通过旁边的像素来做判断根本无法识别，这就的用到一些算法知识了。

感觉说的不着边际，简单的意思就是滤波没什么神秘的，接下来举几个滤波实际例子就能说明：

--------均值滤波--------

从另一篇博文可以知道，对于“核”这个概念的知识，还有怎么操作定义核，这里就不啰嗦了。不知道的直接去看另一篇博文：掩膜操作实例分析

均值就是平均的意思，对邻近的像素做个平均：x[i][j]=1/9(y[i-1][j-1]+y[i-1][j]+y[i-1][j+1]+y[i][j-1]+y[i][j]+y[i][j+1]+y[i+1][j-1]+y[i+1][j]+y[i+1][j+1])

手写实例：

 其实这些例子都不难，看懂理论基础和opencv源代码，自己再去写很简单，有时间这部分会补上，现在为了想快点把opencv学一遍。。。。所以偷懒了

API：

 

--------中值滤波--------

 

 

 

 中值滤波就是旁边的值进行排序取中间值x[i][j]=Middle(y[i-1][j-1]，y[i-1][j]，y[i-1][j+1]，y[i][j-1]，y[i][j]，y[i][j+1]，y[i+1][j-1]，y[i+1][j]，y[i+1][j+1]);

 

手写实例：

 

 

 

API：

 

--------高斯滤波--------

 

这个得重点说明一下，刚开始学的时候真的感觉好难，数学公式也看不懂，源代码看着就烦。。。

公式大家自己看书本或者百度百科吧，这里只说一下我自己的理解

 

 

二维图像 

三维图像

 不说专业名词，看着这幅图就是对称的，并且两遍递减，想象一下这个函数去表示“核”是什么样子？

这是随便画的一个图像，从内到外逐渐减小，如果运用到图像“卷积”之后，那就是距离处理点越近那就权值越大了。

卷积：在数字图像上就是，简单理解就是图像和内核的乘积！

至于这个权值是怎么计算的，那就是直接把坐标带入高斯公式就行了，想深入了解的就去看一些数学的书籍。

 

 

原公式

简化公式

 具体公式不介绍了，看百度、或者其他资料。。。

opencv源码分析：

 

/****************************************************************************************\
                                     Gaussian Blur
\****************************************************************************************/

cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )
{
    //----自定义了几个简单的内核数据-------//
    const int SMALL_GAUSSIAN_SIZE = 7;
    static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
    {
        {1.0f},
        {0.25f, 0.5f, 0.25f},
        {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
        {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
    };
   //-------当 n是奇数 && n<7 && sigma<=0 --->>>直接取small_gaussian_tab里的数据-----//
   //-------当不满足上面的情况时候，就自己按照公式计算，其实上面的值也是按照公式计算的----//
    const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
        small_gaussian_tab[n>>1] : 0;

    CV_Assert( ktype == CV_32F || ktype == CV_64F );
    Mat kernel(n, 1, ktype);//计算的核是n行1列，后面有程序会合并成二维数组
    //----保证程序可以执行float、double两种类型都可以----//
    float* cf = (float*)kernel.data;
    double* cd = (double*)kernel.data;
//-------当sigma取值，这个是高斯函数规定的------//
    double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
    double scale2X = -0.5/(sigmaX*sigmaX);//高斯公式里面的一个系数
    double sum = 0;

    int i;
    for( i = 0; i < n; i++ )
    {
        double x = i - (n-1)*0.5;
        //------就是上面说的看自己定义还是直接读取------//
        double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);
        if( ktype == CV_32F )
        {
            cf[i] = (float)t;
            sum += cf[i];
        }
        else
        {
            cd[i] = t;
            sum += cd[i];
        }
    }

    sum = 1./sum;//系数归一化
    //---归一化系数之后再把系数放在内核里面---//
    for( i = 0; i < n; i++ )
    {
        if( ktype == CV_32F )
            cf[i] = (float)(cf[i]*sum);
        else
            cd[i] *= sum;
    }

    return kernel;
}

//-----函数是处理x核、y核、size核，其中之间的关系和数据变换--------//
cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize,
                                        double sigma1, double sigma2,
                                        int borderType )
{
    int depth = CV_MAT_DEPTH(type);
    if( sigma2 <= 0 )
        sigma2 = sigma1;

    // automatic detection of kernel size from sigma
    if( ksize.width <= 0 && sigma1 > 0 )
        ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
    if( ksize.height <= 0 && sigma2 > 0 )
        ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1;

    CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 &&
        ksize.height > 0 && ksize.height % 2 == 1 );

    sigma1 = std::max( sigma1, 0. );
    sigma2 = std::max( sigma2, 0. );

    Mat kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) );
    Mat ky;
    if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )
        ky = kx;
    else
        ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) );

    return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );
}


void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
                   double sigma1, double sigma2,
                   int borderType )
{
    Mat src = _src.getMat();
    _dst.create( src.size(), src.type() );
    Mat dst = _dst.getMat();
//-----处理边界问题----//
    if( borderType != BORDER_CONSTANT )
    {
        if( src.rows == 1 )
            ksize.height = 1;
        if( src.cols == 1 )
            ksize.width = 1;
    }

    if( ksize.width == 1 && ksize.height == 1 )
    {
        src.copyTo(dst);
        return;
    }

#ifdef HAVE_TEGRA_OPTIMIZATION
    if(sigma1 == 0 && sigma2 == 0 && tegra::gaussian(src, dst, ksize, borderType))
        return;
#endif

    Ptr<FilterEngine> f = createGaussianFilter( src.type(), ksize, sigma1, sigma2, borderType );
    f->apply( src, dst );
}

总结的注意点：

A.高斯核的计算是两种方法，第一种是经过计算放在一个数组里面，第二种是直接用高斯公式计算。

 

 static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
 10     {
 11         {1.0f},
 12         {0.25f, 0.5f, 0.25f},
 13         {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
 14         {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
 15     };

 

 

 

B.高斯的内核是合成的，而不是一次型的计算，具体合成过程我没有细看，我们自己手写API的时候基本    都不会去合成，直接计算便可。

 Mat kernel(n, 1, ktype);//计算的核是n行1列，后面有程序会合并成二维数组 

 

 Mat kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) );
 84     Mat ky;
 85     if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )
 86         ky = kx;
 87     else
 88         ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) );
 89 
 90     return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );

 

 

 

C.滤波器的引擎里面有很多函数，通过一些参数的匹配来分别解决哪个滤波。就是一个封装多个滤波的函    数，看源代码就会发现。源代码在filter.cpp文件，想看的可以自己去看

 Ptr<FilterEngine> f = createGaussianFilter( src.type(), ksize, sigma1, sigma2, borderType ); 2 f->apply( src, dst );   

 

--------双边滤波--------

公式百度吧。。。。。

记得“高斯滤波”是考虑了空间问题，使图像的滤波提升了一个档次（相对于均值滤波），但是高斯滤波对比中值滤波不一定哪个更好，看下面例子：

 

 

 这图片有两个突出的点，一个是0，一个是255，黑白点，其它的像素点都是110左右，这一看就知道是噪声，下面用滤波函数操作这个图像。

均值：0->80+

中值：0->110

高斯：0->90

这个例子举的不明显，但是意图就是高斯函数只能过滤空间的部分，但是距离的部分不能过滤，就是像素差别特别大的时候无法过滤，看网上有个a-截断滤波

高斯+截断=双边    这是简单表示的。

 

 双边滤波就是把那些怪异的点去除之后在进行高斯滤波：

 opencv源代码分析：

/****************************************************************************************\
                                   Bilateral Filtering
\****************************************************************************************/

namespace cv
{

static void
bilateralFilter_8u( const Mat& src, Mat& dst, int d,
                    double sigma_color, double sigma_space,
                    int borderType )
{
    int cn = src.channels();//原图通道数
    int i, j, k, maxk, radius;
    Size size = src.size();//原图大小

    CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
        src.type() == dst.type() && src.size() == dst.size() &&
        src.data != dst.data );//只能处理单通道和三通道的uchar
//-------限定sigma的取值--------//
    if( sigma_color <= 0 )
        sigma_color = 1;
    if( sigma_space <= 0 )
        sigma_space = 1;
//------高斯因子,计算方便-------//
    double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
    double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
//-----d 表示滤波时像素邻域直径，d为负时由 sigaColor计算得到；d>5时不能实时处理--//
    if( d <= 0 )
        radius = cvRound(sigma_space*1.5);//float->int
    else
        radius = d/2;
//---保证直径是基整数---//
    radius = MAX(radius, 1);
    d = radius*2 + 1;

    Mat temp;
    copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
//---直接用公式计算权值，只有用到直接调用----//
    vector<float> _color_weight(cn*256);//距离权值范围为channels*256
    vector<float> _space_weight(d*d);//空间权值存储
    vector<int> _space_ofs(d*d);//位置便宜存储->用于计算点的位置
    float* color_weight = &_color_weight[0];
    float* space_weight = &_space_weight[0];
    int* space_ofs = &_space_ofs[0];

    // initialize color-related bilateral filter coefficients
    for( i = 0; i < 256*cn; i++ )//把空间可能的权值权值全部存储起来，用的时候直接找
        color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);//省的到时候计算

    // ------计算核的值和偏移值----//
    for( i = -radius, maxk = 0; i <= radius; i++ )
        for( j = -radius; j <= radius; j++ )
        {
            double r = std::sqrt((double)i*i + (double)j*j);//高斯公式系数之一
            if( r > radius )
                continue;
            space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);//空间权值
            space_ofs[maxk++] = (int)(i*temp.step + j*cn);//偏移值计算
        }
//------------------进行滤波操作----------//
    for( i = 0; i < size.height; i++ )
    {
        //-----sptr核中心的像素地址，同时是输入像素变量----//
        //-----dptr表示第i行的第一个像素地址，同时是输出像素变量-----//
        const uchar* sptr = temp.data + (i+radius)*temp.step + radius*cn;
        uchar* dptr = dst.data + i*dst.step;

        if( cn == 1 )
        {
            for( j = 0; j < size.width; j++ )
            {
                float sum = 0, wsum = 0;
                int val0 = sptr[j];
                for( k = 0; k < maxk; k++ )//maxk=d就是直径的大小
                {
                    int val = sptr[j + space_ofs[k]];//核对应的像素值
                    float w = space_weight[k]*color_weight[std::abs(val - val0)];//总权值
                    sum += val*w;//权值*像素值总和
                    wsum += w;//权值和，为了下面归一化
                }
                // overflow is not possible here => there is no need to use CV_CAST_8U
                dptr[j] = (uchar)cvRound(sum/wsum);//权值归一化
            }
        }
        else
        {
            assert( cn == 3 );
            for( j = 0; j < size.width*3; j += 3 )//三个通道(一个像素)同时操作
            {
                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];//三个通道代表一个像素
                for( k = 0; k < maxk; k++ )//直径的大小
                {
                    const uchar* sptr_k = sptr + j + space_ofs[k];//核对应的像素地址
                    int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
                    //-----三个通道的差值和作为计算、所以三个通道的权重是相同的------//
                    float w = space_weight[k]*color_weight[std::abs(b - b0) +
                        std::abs(g - g0) + std::abs(r - r0)];
                    sum_b += b*w; sum_g += g*w; sum_r += r*w;
                    wsum += w;
                }
                //----下面就是BGR三个通道的调整了-------//
                wsum = 1.f/wsum;
                b0 = cvRound(sum_b*wsum);
                g0 = cvRound(sum_g*wsum);
                r0 = cvRound(sum_r*wsum);
                dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
            }
        }
    }
}

 

总结的注意点：

A.双边滤波如果是三通道的BGA的话，那是三个通道同时进行的，三个通道的差值求和，所以三个通道是一个权重。

为对于其他其他的滤波方式：

 

 B.这个滤波函数使用了指针访问像素，但是这个指针是:*(M.data+i*M.step[0]+j*M.channels())，记得上次我们使用的指针是     M.ptr<Vec3b>(i)[j]访问的。下面就详细说明一下这个函数：

        (a).其中M.data是M这个Mat类的首地址，也就是一幅图像的首地址。

        (b).M.step[],这个函数比较绕人，三维空间咱们不说，就说我们使用的二维数据。

              step[0] :表示图片一行数据字节。

              step[1] :表示图片一个像素的字节。

　　　　　扩展几个知识点：

              elemSize()  :一个像素的字节　== step[1]

  　　　　 elemSize1():一个通道的字节

              比如：我们图片是16X16,CV_8UC3-->>>elemSize1()=1，elemSize()=step[1]=1*3，step[0]=1*3*16

                      我们图片是25*25,CV_16UC3-->>>elemSize1()=2，elemSize()=step[1]=2*3，step[0]=2*3*25

              8位的uchar是一个字节，16位的uchar是2个字节。。。

              16X16,CV_8UC3的图片 *(M.data+i*M.step[0]+j*M.channels()):*(0+1*3*16*I+3*J)，这里省略step[1]=1

 　　　　  我们之前写的Mat.ptr和Mat.at内部已经帮我们做了这些事，现在直接写的话只有通过字节的地址操作，

　　　　　　　　就像char a[100];char*p=a;p++就可以操作a[100],int a[100];int* p=a;p+=4;

 

 参考：

        http://blog.csdn.net/qianqing13579/article/details/45318279#comments

 　　　  http://blog.csdn.net/xiaowei_cqu/article/details/7785365

 　　　  http://blog.csdn.net/zhaocj/article/details/39520187

还有其它参考，没来及记录，入您发现侵犯您的权利，请告之，立马清除！