图像处理--高斯滤波

原文http://blog.csdn.net/l_inyi/article/details/8915116

部分内容整理自：

http://blog.csdn.net/jianxiong8814/article/details/1562728

http://www.cnblogs.com/pegasus/archive/2011/05/20/2052031.html

http://blog.sina.com.cn/s/blog_640577ed0100yz8v.html

http://blog.csdn.net/pearl333/article/details/8696307

http://www.cnblogs.com/tornadomeet/tag/opencv/

关于滤波

滤波通常是用卷积或者相关来描述，而线性滤波一般是通过卷积来描述的。他们非常类似，但是还是会有不同。下面我们来根据相关和卷积计算过程来体会一下他们的具体区别：

卷积的计算步骤：

1）卷积核绕自己的核心元素顺时针旋转180度

2）移动卷积核的中心元素，使它位于输入图像待处理像素的正上方

3）在旋转后的卷积核中，将输入图像的像素值作为权重相乘

4）第三步各结果的和做为该输入像素对应的输出像素

相关的计算步骤：

1）移动相关核的中心元素，使它位于输入图像待处理像素的正上方

2）将输入图像的像素值作为权重，乘以相关核

3）将上面各步得到的结果相加做为输出

可以看出他们的主要区别在于计算卷积的时候，卷积核要先做旋转。而计算相关过程中不需要旋转相关核。

例如：magic(3) =[8 1 6;3 5 7;4 9 2]，旋转180度后就成了[2 9 4;7 5 3;6 1 8]

 

高斯滤波

高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。

高斯滤波的具体操作是：用一个模板(或称卷积、掩模)扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值。

若使用3×3模板，则计算公式如下g(x,y)={f(x-1,y-1)+f(x-1,y+1)+f(x+1,y-1)+f(x+1,y+1)+[f(x-1,y)+f(x,y-1)+f(x+1,y)+f(x,y+1)]*2+f(x,y)*4}/16;

其中，f(x,y)为图像中(x,y)点的灰度值，g(x,y)为该点经过高斯滤波后的值。

 

 

1、高斯分布:

一维高斯分布:

 

 

二维高斯分布:

 

2、高斯核

理论上，高斯分布在所有定义域上都有非负值，这就需要一个无限大的卷积核。实际上，仅需要取均值周围3倍标准差内的值，以外部份直接去掉即可。如下图为一个标准差为1.0的整数值高斯核。

3、应用：

高斯滤波后图像被平滑的程度取决于标准差。它的输出是领域像素的加权平均，同时离中心越近的像素权重越高。因此，相对于均值滤波（mean filter）它的平滑效果更柔和，而且边缘保留的也更好。

高斯滤波被用作为平滑滤波器的本质原因是因为它是一个低通滤波器（让某一频率以下的信号分量通过，而对该频率以上的信号分量大大抑制）

 

4、特征：   

       1 ）：一个高斯函数跟另外一个高斯函数的卷积仍然是一个高斯函数，A*B=C   C的标准差的平方是A和B的标准差的平方和，也就是说卷积后的高斯函数更宽，模糊的效果更明显（直观上看，连续做高斯模糊运算，图像会越来越模糊） 。

  二维高斯函数卷积可以分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积结果与方向垂直的相同一维高斯函数卷积．因此，二维高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长．

      2）：高斯函数的傅立叶变换仍然是一个高斯函数，如果原来的高斯函数越宽（标准差越大），变换后的高斯函数就越窄（标准差越小），也就是说一个越宽的高斯函数，低通（高阻）滤波的效果越明显，处理后的图像的细节就越不清楚（更模糊）。   
    
  要对数字图像做高斯模糊，就是用一个符合高斯函数分布的卷积核对数字图像做卷积运算。   
  要确定的有标准差的大小，卷积核的大小，最后的比例系数的大小。   

σ越大，高斯滤波器的频带就越宽，平滑程度就越好
    
  一个标准差为1.4的高斯5x5的卷积核：   
    
  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         
    
  最后乘以比例系数   1/115   

 

在opencv中使用高斯滤波

 

复合形式函数原型为：

[cpp]  view plain copy

1. void cvSmooth(const CvArr* src, CvArr* dst,intsmoothtype=CV_GAUSSIAN,int param1=3, int param2=0,double param3=0, doubleparam4=0 );  

该函数前三个参数很容易理解，至于后四个参数以下进行分析。

      1） 如果指定param1和param2，则代表核函数的行列，即为滤波窗口的宽度和高度；

      2） Param3：高斯卷积的Sigma值

      3） 如果用户希望采用非对称的高斯核，则引入param4，最后两个参数分别代表水平核以及垂直核维数；

      4） 如果param3没有给出，则有前两个参数param1和param2计算出Sigma。这里的根据是高斯分布的特点（如图所示，数值分布在（μ—3σ,μ+3σ)中的概率为0.9974），如果核矩阵更大，那么相应的Sigma也更大，相反，如果Sigma更大，那么核矩阵覆盖范围也更大。具体到OpenCv下，用如下公式进行计算(根据其源代码显示)。

  

 5）同样的根据这个公式可知，如果param1和param2为0（或者没有给出），那么滤波窗口的尺寸，则有后两个参数代表的Sigma来确定。

 单独的高斯滤波函数声明为

[cpp]  view plain copy

1. void GaussianBlur(InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY=0, int borderType=BORDER_DEFAULT ) ;  

      功能：对输入的图像src进行高斯滤波后用dst输出。

      参数：src和dst当然分别是输入图像和输出图像。

                  Ksize为高斯滤波器模板大小

                  sigmaX和sigmaY分别为高斯滤波在横线和竖向的滤波系数。

                  borderType为边缘扩展点插值类型。

       接下来的工作就是进入GaussianBlur函数内部，跟踪其函数代码，经过分析，在该函数内部调用了很多其他的函数，其调用的函数层次结构如下图所示：

      这里我们分析源代码不需要深入到最底层，我们只需分析到函数createSeparableLinearFilter和getGaussianKernel这一层。

 

      从函数调用层次结构图可以看出，要分析函数GaussianBlur，必须先分析其调用过的内部函数。

      

因此首先分析函数getGaussianKernel。

      功能：返回一个ksize*1的数组，数组元素满足高斯公式：

     其中只有系数alpha和参数sigma未知，sigma的求法为：

     如果输入sigma为非正，则计算公式为：sigma =0.3*((ksize-1)*0.5 - 1) + 0.8 .

     如果输入sigma为正，则就用该输入参数sigma。

     最后alpha为归一化系数，即计算出的ksize个数之和必须为1，所以后面只需求ksize个数，计算其和并求倒即可。

    其源码与注释如下：

[cpp]  view plain copy

1. <span style="font-size:18px;">cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )  
2. {  
3.     const int SMALL_GAUSSIAN_SIZE = 7;  
4.     static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =  
5.     {  
6.         {1.f},  
7.         {0.25f, 0.5f, 0.25f},  
8.         {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},  
9.         {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}  
10.     };  
11.   
12.   
13.     const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?  
14.         small_gaussian_tab[n>>1] : 0;  
15.     CV_Assert( ktype == CV_32F || ktype == CV_64F );//确保核元素为32位浮点数或者64位浮点数  
16.     Mat kernel(n, 1, ktype);//建立一个n*1的数组kernel，一个Mat矩阵包括一个矩阵头和一个指向矩阵元素的指针  
17.     float* cf = (float*)kernel.data;//定义指针cf指向kernel单精度浮点型数据  
18.     double* cd = (double*)kernel.data;//定义指针cd指向kernerl双精度浮点型数据  
19.     double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;//当sigma小于0时，采用公式得到sigma(只与n有关)  
20.     double scale2X = -0.5/(sigmaX*sigmaX);//高斯表达式后面要用到  
21.     double sum = 0;  
22.     int i;  
23.     for( i = 0; i < n; i++ )  
24.     {  
25.         double x = i - (n-1)*0.5;  
26.         //如果自己算其核的话，就常用公式exp(scale2X*x*x)计算，否则就用固定系数的核  
27.         double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);  
28.         if( ktype == CV_32F )  
29.         {  
30.             cf[i] = (float)t;//单精度要求时存入cf数组中  
31.             sum += cf[i];//进行归一化时要用到  
32.         }  
33.         else  
34.         {  
35.             cd[i] = t;//双精度时存入cd数组中  
36.             sum += cd[i];  
37.         }  
38.     }  
39.     sum = 1./sum;//归一化时核中各元素之和为1  
40.     for( i = 0; i < n; i++ )  
41.     {  
42.         if( ktype == CV_32F )  
43.             cf[i] = (float)(cf[i]*sum);//归一化后的单精度核元素  
44.         else  
45.             cd[i] *= sum;//归一化后的双精度核元素  
46.     }  
47.     return kernel;//返回n*1的数组，其元素或是单精度或是双精度，且符合高斯分布  
48. }</span>  

下面该分析函数createSeparableLinearFilter了。

    功能为：创建一个图像滤波其引擎类，其主要处理的是原图像和目标图像数据格式的统以及滤波器核的合成。

其源码及注释如下：

[cpp]  view plain copy

1. cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize,  
2.     double sigma1, double sigma2,  
3.     int borderType )  
4. {  
5.     int depth = CV_MAT_DEPTH(type);//取数组元素的深度  
6.     if( sigma2 <= 0 )  
7.         sigma2 = sigma1;//当第3个参数为非正时，取其与第二个参数相同的值  
8.     // automatic detection of kernel size from sigma  
9.   
10.     if( ksize.width <= 0 && sigma1 > 0 )//当滤波器核的宽非正时，其宽要重新经过计算  
11.   
12.         ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;  
13.     if( ksize.height <= 0 && sigma2 > 0 )  
14.   
15.         ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1;  
16.     CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 &&  
17.         ksize.height > 0 && ksize.height % 2 == 1 );//确保核宽和核高为正奇数  
18.     sigma1 = std::max( sigma1, 0. );//sigma最小为0  
19.     sigma2 = std::max( sigma2, 0. );  
20.     Mat kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) );//得到x方向一维高斯核  
21.     Mat ky;  
22.     if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )  
23.         ky = kx;//如果核宽和核高相等，且两个sigma相差很小的情况下，y方向的高斯核去与x方向一样，减少计算量  
24.     else  
25.         ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) );//否则计算y方向的高斯核系数  
26.     return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );//返回2维图像滤波引擎  
27. }  

最后来看真正的高斯滤波函数GaussianBlur：功能：对输入图像_src进行滤波得到输出图像_dst,滤波核大小为ksize，滤波参数由sigma1和sigma2计算出，边缘扩展模式为borderType.其源代码和注释如下：

[cpp]  view plain copy

1. void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,  
2. double sigma1, double sigma2,  
3. int borderType )  
4.   
5. Mat src = _src.getMat();//创建一个矩阵src，利用_src的矩阵头信息  
6. _dst.create( src.size(), src.type() );//构造与输入矩阵同大小的目标矩阵  
7. Mat dst = _dst.getMat();//创建一个目标矩阵  
8.   
9. if( ksize.width == 1 && ksize.height == 1 )  
10. {  
11.     src.copyTo(dst);//如果滤波器核的大小为1的话，则说明根本就不用滤波，输出矩阵与输入矩阵完全相同  
12.     return;  
13. }  
14. if( borderType != BORDER_CONSTANT )//当边缘扩展不是常数扩展时  
15. {  
16.     if( src.rows == 1 )  
17.         ksize.height = 1;//如果输入矩阵是一个行向量，则滤波核的高强制为1  
18.     if( src.cols == 1 )  
19.         ksize.width = 1;//如果输入矩阵是一个列向量，则滤波核的宽强制为1  
20. }  
21.   
22. Ptr<FilterEngine> f = createGaussianFilter( src.type(), ksize, sigma1, sigma2, borderType );  
23. f->apply( src, dst );//调用引擎函数，完成将输入矩阵src高斯滤波为输出矩阵dst  

具体的可实现的高斯滤波

源码和注释如下：

[cpp]  view plain copy

1. #include "highgui.h"    
2. #include "cv.h"    
3.   
4. #pragma comment(linker, "/subsystem:windows /entry:mainCRTStartup")  
5.   
6. void gaussianFilter(uchar* data, int width, int height)     
7. {       
8.     int i, j, index, sum;     
9.     int templates[9] = { 1, 2, 1,     
10.         2, 4, 2,     
11.         1, 2, 1 };//模板的值     
12.     sum = height * width * sizeof(uchar);//图像所占内存的大小     
13.     uchar *tmpdata = (uchar*)malloc(sum);     
14.     memcpy((char*)tmpdata,(char*)data, sum);     
15.     for(i = 1;i < height - 1;i++)     
16.     {     
17.         for(j = 1;j < width - 1;j++)     
18.         {               
19.             index = sum = 0;     
20.             for(int m = i - 1;m < i + 2;m++)     
21.             {     
22.                 for(int n = j - 1; n < j + 2;n++)     
23.                 {     
24.                     sum += tmpdata[m * width + n] * templates[index++];  //处理     
25.                 }     
26.             }     
27.             data[i * width + j] = sum / 16;     
28.   
29.         }     
30.     }     
31.     free(tmpdata);     
32. }     
33.   
34. void imgOperate( IplImage* image )    
35. {    
36.     cvNamedWindow( "image-in", CV_WINDOW_AUTOSIZE );    
37.     cvNamedWindow( "image-access", CV_WINDOW_AUTOSIZE);    
38.     cvNamedWindow( "image-out", CV_WINDOW_AUTOSIZE);    
39.     cvShowImage( "image-in", image );    
40.   
41.     //将色彩图像强制转化为灰度图像    
42.     IplImage* pGrayImg=NULL;   
43.     IplImage* pImg=NULL;  
44.   
45.     pGrayImg=cvCreateImage(cvGetSize(image),8,1);    
46.     pImg=cvCreateImage(cvGetSize(image),8,1);   
47.     cvCvtColor(image,pGrayImg,CV_RGB2GRAY);    
48.     //  cvSmooth( image, out, CV_GAUSSIAN, 5,5 );  
49.     //添加高斯噪声  
50.     cvZero(pImg);  
51.     CvRNG rng = cvRNG(-1);  //初始化随机数发生器  
52.     cvRandArr(&rng, pImg, CV_RAND_NORMAL, cvScalarAll(0), cvScalarAll(15));  
53.     cvAdd(pGrayImg, pImg,pImg);  
54.     cvShowImage("image-access",pImg);  
55.     gaussianFilter((unsigned char*)pImg->imageData,pImg->width,pImg->height);    
56.     cvShowImage( "image-out", pImg );    
57.   
58.     cvReleaseImage( &pGrayImg );    
59.     cvWaitKey( 0 );     
60.     cvDestroyWindow("image-in" );    
61.     cvDestroyWindow("image-out" );        
62. }    
63.   
64. int main()    
65. {    
66.     IplImage* img = cvLoadImage("d:\\demo.jpg");    
67.     imgOperate( img );    
68.     cvReleaseImage( &img );    
69.     return 0;    
70. }