高斯濾波是一種線性平滑濾波，適用於消除高斯噪聲，廣泛應用於圖像處理的減噪過程。通俗的講，高斯濾波就是對整幅圖像進行加權平均的過程，每一個像素點的值，都由其本身和鄰域內的其他像素值經過加權平均後得到。
高斯濾波的具體操作是：用一個模板(或稱卷積、掩模)掃描圖像中的每一個像素，用模板確定的鄰域內像素的加權平均灰度值去替代模板中心像素點的值。
一般的模板為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)為該點經過高斯濾波後的值。
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：一個高斯函數跟另外一個高斯函數的卷積仍然是一個高斯函數，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
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
/////////////////////////////////////////////////////////////
// 代碼 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 {
for(j=0; j {
// 計算對應模板上位置的像素在原圖像中的位置
py=y-tw/2+i;
px=x-tw/2+j;
// 加權求和 + 加權求和
m+=imageBuf0[py][px*4+cn] * templt[i*tw+j];
}
}
return m; //返回結果
}
* 功能: 使用模板對彩色圖鄰域進行運算
* 參數: 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
for(j=0; 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;
* 功能: 彩色圖像的高斯平滑處理
* 參數: 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)
scale = 16; // 因為模板的衰減因子為16 (1 + 2 + 1 + 2 + 4 + 2 + 1 + 2 + 1)
//依次對原圖像的每個像素進行處理
for(y=1; y for(x=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;
}
for(y=1; y
{
//利用高斯模板對鄰域進行處理
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);
}
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) ; // 這裡是生成高斯模板, 但這個函數好像有點問題的
*
* /函數名稱：
* 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;
// MakeGauss返回窗口的長度，利用此變量計算窗口的半長
nHalfLen = nWindowSize / 2;
// 這裡為了加快速度, 先對行再對列
// x方向進行濾波
for(y=0; y {
for(x=0; 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 {
for(y=0; 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 ;
}
}
// x方向進行濾波
for(y=0; y
for(x=0; 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
for(y=0; 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;
}
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;
*
* /函數名稱：
* MakeGauss()
*
* /輸入參數:
* double sigma - 高斯函數的標准差
* double **pdKernel - 指向高斯數據數組的指針
* int *pnWindowSize - 數據的長度
*
* /返回值:
* 無
*
* /說明:
* 這個函數可以生成一個一維的高斯函數的數字數據，理論上高斯數據的長度應
* 該是無限長的，但是為了計算的簡單和速度，實際的高斯數據只能是有限長的
* pnWindowSize就是數據長度
*
*************************************************************************
*/
void MakeGauss(double sigma, double **pdKernel, int *pnWindowSize)
{
// 循環控制變量
int i ;
// 數組的中心點
int nCenter;
// 數組的某一點到中心點的距離
double dDis ;
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)
// 中間變量
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;
}
}
//////////////////////////////////////////////////////////////////////
{
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;
}
}
//////////////////////////////////////////////////////////////////////
高斯濾波（高斯平滑） 被以下幾個問題困擾：
- 給定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;
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;
* 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;
radius = std_dev * 2;
/* go out 『radius』 in each direction */
matrix_length = 2 * ceil (radius – 0.5) + 1;
注釋講的很清楚了，基本上就是認為sigma應該等於窗口半徑的一半。
看完這三分源碼，結論就是，關於sigma和半徑，你愛怎麼算就怎麼算吧，差不多就行。。。
第三個問題是可分離濾波器：
由於高斯函數可以寫成可分離的形式，因此可以采用可分離濾波器實現來加速。所謂的可分離濾波器，就是可以把多維的卷積化成多個一維卷積。具體到二維的高斯濾波，就是指先對行做一維卷積，再對列做一維卷積。這樣就可以將計算復雜度從O(M*M*N*N)降到O(2*M*M*N)，M，N分別是圖像和濾波器的窗口大小。問題是實現時候怎麼計算一維的卷積核呢？
其實很簡單，按照前面計算出來的窗口大小，計算所有離散點上一維高斯函數的權值，最後別忘了將權值之和歸一化到1.有碼有真相，來自opencv：
for( i = 0; i <= n/2; i++ )
{
double t = fixed_kernel ? (double)fixed_kernel[i] : exp(scale2X*i*i);
if( type == CV_32FC1 )
{
cf[(n/2+i)*step] = (float)t;
sum += cf[(n/2+i)*step]*2;
}
else
{
cd[(n/2+i)*step] = t;
sum += cd[(n/2+i)*step]*2;
}
}
{
double t = fixed_kernel ? (double)fixed_kernel[i] : exp(scale2X*i*i);
if( type == CV_32FC1 )
{
cf[(n/2+i)*step] = (float)t;
sum += cf[(n/2+i)*step]*2;
}
else
{
cd[(n/2+i)*step] = t;
sum += cd[(n/2+i)*step]*2;
}
}
sum = 1./sum;
for( i = 0; i <= n/2; i++ )
{
if( type == CV_32FC1 )
cf[(n/2+i)*step] = cf[(n/2-i)*step] = (float)(cf[(n/2+i)*step]*sum);
else
cd[(n/2+i)*step] = cd[(n/2-i)*step] = cd[(n/2+i)*step]*sum;
}
for( i = 0; i <= n/2; i++ )
{
if( type == CV_32FC1 )
cf[(n/2+i)*step] = cf[(n/2-i)*step] = (float)(cf[(n/2+i)*step]*sum);
else
cd[(n/2+i)*step] = cd[(n/2-i)*step] = cd[(n/2+i)*step]*sum;
}
沒有留言:
張貼留言