高斯滤波器参数的确定

来源：互联网发布：手机如何淘宝客赚佣金编辑：程序博客网时间：2024/04/29 12:21

以下文字引用自：http://blog.csdn.net/cay22/archive/2010/05/17/5600836.aspx

http://vipbase.net/ipbook/chap03.htm

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

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

一般的模板为3×3或5×5大小，其权值分布如下图：

若使用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)为该点经过高斯滤波后的值。

http://vipbase.net/ipbook/chap03.htm

高斯函数有两个特性：
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

http://blog.csdn.net/jianxiong8814/archive/2007/04/12/1562728.aspx

/////////////////////////////////////////////////////////////

// 代码 1

/**************************************************
* 功能: 使用模板对彩色图邻域进行运算
* 参数: imageBuf为目标图像 w、h为图像大小
*       templt为模板 tw为邻域大小
* x + y为要取得像素的坐标
*       cn为颜色分量编号 0为蓝色 1为绿色 2为红色
**************************************************/
int TempltExcuteCl(BYTE** imageBuf0 + int w + int h + int* templt + int tw + int x + int y + int cn)
{
int i + j;                        //循环变量
int m=0;                      //用来存放加权和
int px + py;
//依次对邻域中每个像素进行运算
for(i=0; i<tw; i++)
{
for(j=0; j<tw; j++)
{
   // 计算对应模板上位置的像素在原图像中的位置
   py=y-tw/2+i;
   px=x-tw/2+j;
   // 加权求和 + 加权求和
   m+=imageBuf0[py][px*4+cn] * templt[i*tw+j];
}
}
return m;                     //返回结果
}

/******************************************************************
* 功能: 彩色图像的高斯平滑处理
* 参数: image0为原图形，image1平滑结果，
* w、h为图像的宽和高
******************************************************************/
void SmoothGaussCl(BYTE* image0 + BYTE* image1 + unsigned int w + unsigned int h)
{
//将图像转化为矩阵形式
BYTE** imageBuf0 = CreatImage(image0 + w + h);
BYTE** imageBuf1 = CreatImage(image1 + w + h);
//设定模板
int templt[9]={1 + 2 + 1 + 2 + 4 + 2 + 1 + 2 + 1};
int x + y + c;
int a;
int scale;

//设定衰减因子
scale = 16; // 因为模板的衰减因子为16 (1 + 2 + 1 + 2 + 4 + 2 + 1 + 2 + 1)

//依次对原图像的每个像素进行处理
for(y=1; y<h-1; y++)
for(x=1; x<w-1; x++)
   for(c=0; c<3; c++)
   {
    //利用高斯模板对邻域进行处理
    a=TempltExcuteCl(imageBuf0 + w + h + templt + 3 + x + y + c);
    a/= scale;
    //过限处理
    a = a>255?255:a;
    a = a<0?0:a;
    imageBuf1[y][x*4+c]=a;
   }

//清理内存
free(imageBuf0);
free(imageBuf1);
}

///////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////

// 代码 2

/*************************************************************************
*
* \函数名称：
*   GaussianSmooth()
*
* \输入参数:
*   unsigned char * pUnchImg    - 指向图象数据的指针
*   int nWidth           - 图象数据宽度
*   int nHeight          - 图象数据高度
*   double dSigma         - 高斯函数的标准差
*   unsigned char * pUnchSmthdImg - 指向经过平滑之后的图象数据
*
* \返回值:
*   无
*
* \说明:
*   为了抑止噪声，采用高斯滤波对图象进行滤波，滤波先对x方向进行，然后对
*   y方向进行。
*
*************************************************************************
*/
void GaussianSmooth(unsigned char *pUnchImg, int nWidth, int nHeight,
          double sigma, unsigned char * pUnchSmthdImg)
{
// 循环控制变量
    int y;
int x;

int i;

// 高斯滤波器的数组长度

int nWindowSize;

// 窗口长度的1/2
int nHalfLen;

// 一维高斯数据滤波器
double *pdKernel ;

// 高斯系数与图象数据的点乘
double dDotMul     ;

// 高斯滤波系数的总和
double dWeightSum     ;

// 中间变量
double * pdTmp ;

// 分配内存
pdTmp = new double[nWidth*nHeight];

// 产生一维高斯数据滤波器
// MakeGauss(sigma, &dKernel, &nWindowSize);
MakeGauss(sigma, &pdKernel, &nWindowSize) ;     // 这里是生成高斯模板, 但这个函数好像有点问题的

// 也可以像代码1 一样固定模板，不用动态生成。

// MakeGauss返回窗口的长度，利用此变量计算窗口的半长
nHalfLen = nWindowSize / 2;

// 这里为了加快速度, 先对行再对列

// x方向进行滤波
for(y=0; y<nHeight; y++)
{
for(x=0; x<nWidth; x++)
{
   dDotMul = 0;
   dWeightSum = 0;
   for(i=(-nHalfLen); i<=nHalfLen; i++)
   {
    // 判断是否在图象内部
    if( (i+x) >= 0 && (i+x) < nWidth )
    {
     dDotMul += (double)pUnchImg[y*nWidth + (i+x)] * pdKernel[nHalfLen+i];
     dWeightSum += pdKernel[nHalfLen+i];
    }
   }
   pdTmp[y*nWidth + x] = dDotMul/dWeightSum ;
}
}

// y方向进行滤波
for(x=0; x<nWidth; x++)
{
for(y=0; y<nHeight; y++)
{
   dDotMul = 0;
   dWeightSum = 0;
   for(i=(-nHalfLen); i<=nHalfLen; i++)
   {
    // 判断是否在图象内部
    if( (i+y) >= 0 && (i+y) < nHeight )
    {
     dDotMul += (double)pdTmp[(y+i)*nWidth + x] * pdKernel[nHalfLen+i];
     dWeightSum += pdKernel[nHalfLen+i];
    }
   }
   pUnchSmthdImg[y*nWidth + x] = (unsigned char)(int)dDotMul/dWeightSum ;
}
}

// 释放内存
delete []pdKernel;
pdKernel = NULL ;

delete []pdTmp;
pdTmp = NULL;
}

/*************************************************************************
*
* \函数名称：
*   MakeGauss()
*
* \输入参数:
*   double sigma                 - 高斯函数的标准差
*   double **pdKernel          - 指向高斯数据数组的指针
*   int *pnWindowSize          - 数据的长度
*
* \返回值:
*   无
*
* \说明:
*   这个函数可以生成一个一维的高斯函数的数字数据，理论上高斯数据的长度应
*   该是无限长的，但是为了计算的简单和速度，实际的高斯数据只能是有限长的
*   pnWindowSize就是数据长度
*
*************************************************************************
*/
void MakeGauss(double sigma, double **pdKernel, int *pnWindowSize)
{
// 循环控制变量
int i   ;

// 数组的中心点
int nCenter;

// 数组的某一点到中心点的距离
double dDis ;

double PI = 3.14159;
// 中间变量
double dValue;
double dSum ;
dSum = 0 ;

// 数组长度，根据概率论的知识，选取[-3*sigma, 3*sigma]以内的数据。
// 这些数据会覆盖绝大部分的滤波系数
// ceil 的功能是返回不小于其参数的最小整数
// 例如 ceil(6.4) 的值为 7;
*pnWindowSize = 1 + 2 * ceil(3 * sigma);

// 中心
nCenter = (*pnWindowSize) / 2;

// 分配内存
*pdKernel = new double[*pnWindowSize] ;

// g(x)=exp( -x^2/(2 sigma^2)

for(i=0; i< (*pnWindowSize); i++)
{
dDis = (double)(i - nCenter);
dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma)) / (sqrt(2 * PI) * sigma );
(*pdKernel)[i] = dValue ;
dSum += dValue;
}

// 归一化
for(i=0; i<(*pnWindowSize) ; i++)
{
(*pdKernel)[i] /= dSum;
}
}
//////////////////////////////////////////////////////////////////////

http://www.cvchina.info/tag/%E9%AB%98%E6%96%AF%E6%BB%A4%E6%B3%A2/

高斯滤波（高斯平滑）被以下几个问题困扰：

给定sigma，即标准偏差，怎么确定离散化后滤波器的窗口大小？
给定窗口大小，怎么计算高斯核的sigma，即标准方差？
怎么实现可分离滤波器？

我尝试结合三份源码，做个小小的总结。

三份源码分别是：

opencv中的cvfilter.cpp
autopano-sift-c中的GaussianConvolution.c
GIMP中的blur-gauss.c和unsharp-mask.c

在图像处理中，高斯滤波一般有两种实现方式，一是用离散化窗口滑窗卷积，另一种通过傅里叶变换。最常见的就是第一种滑窗实现，只有当离散化的窗口非常大，用滑窗计算量非常大（即使用可分离滤波器的实现）的情况下，可能会考虑基于傅里叶变化的实现方法。

这里我们只讨论第一种方法。二维高斯函数的形式是这样的：

有着如下的形状，形状很激凸，怎么会叫高斯平滑呢，分明是高斯激凸嘛：

基本上，离散化的主旨就是保留高斯函数中心能量最集中的中间部分，忽略四周能量很小的平坦区域。下面结合三份源码，看看现实世界里的高斯平滑到底长的什么样子。

首先是第一个问题：给定sigma，怎么计算窗口大小？

直接上opencv的源码，在cvFilter函数中：

param1 = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;

opencv认为半径为3*sigma的窗口就是高斯函数能量最集中的区域。（为什么在图像深度不是8U的时候，使用4*sigma半径的窗口就不得而知了，有高人指点一下？）

autopan0-sift-c是图像拼接软件hugin里面的sift实现，在实现DoG的时候需要做不同尺度的高斯平滑，实现如下，在GaussianConvolution_new1函数中：

dim = 1 + 2 * ((int) (3.0 * sigma));

可见autopano也是实现的3*sigma半径的窗口。

在GIMP里，实现比较奇特，在blur_gauss.c的make_rle_curve函数里面，

const gdouble sigma2 = 2 * sigma * sigma;
const gdouble l = sqrt (-sigma2 * log (1.0 / 255.0));
int n = ceil (l) * 2;
if ((n % 2) == 0)
n += 1;

我是没看懂那个log (1.0 / 255.0)是干嘛的。。。惭愧。效果来看，这个实现的窗口半径是约等于2.2*sigma。

然后是第二个问题：给定窗口大小，怎么计算sigma？

opencv的实现，在cvFilter.cpp的init_gaussian_kernel函数中：

sigmaX = sigma > 0 ? sigma : (n/2 – 1)*0.3 + 0.8;

再次不可解。。。乘以0.3还可以接受，加上0.8是为嘛啊？

autopano没有实现这个特性。

GIMP的实现：

/* we want to generate a matrix that goes out a certain radius
* from the center, so we have to go out ceil(rad-0.5) pixels,
* inlcuding the center pixel. Of course, that’s only in one direction,
* so we have to go the same amount in the other direction, but not count
* the center pixel again. So we double the previous result and subtract
* one.
* The radius parameter that is passed to this function is used as
* the standard deviation, and the radius of effect is the
* standard deviation * 2. It’s a little confusing.
*/
radius = fabs (radius) + 1.0;

std_dev = radius;
radius = std_dev * 2;
/* go out ‘radius’ in each direction */
matrix_length = 2 * ceil (radius – 0.5) + 1;

注释讲的很清楚了，基本上就是认为sigma应该等于窗口半径的一半。

看完这三分源码，结论就是，关于sigma和半径，你爱怎么算就怎么算吧，差不多就行。。。

第三个问题是可分离滤波器：