目錄
前言
高斯函數
高斯濾波過程
高斯濾波實現
總結
參考資料
附錄
對一幅圖像而言,低頻部分對應整體灰度級的顯示,高頻部分對應著圖像的細節部分.因此去掉低頻部分(或者增強高頻部分)可以銳化圖像,去掉高頻部分(或者增強低頻部分)可以實現模糊/平滑圖像的作用.
去除低頻或者高頻部分可以通過濾波器來完成,但是在時域空間並不能直觀看出所謂的高頻低頻成分,所以就需要做個傅立葉變換,轉換成頻域,然後直接在頻域上把需要去除的頻率對應的「基」的權值置零即可.
高斯濾波器是一種線性低通濾波器,能夠有效的抑制噪聲,平滑圖像.其本質是帶權值的加權均值濾波,權值大小與濾波器中元素到濾波中心距離有關,卷積核中心權重數值最大,並向四周減小,減小的幅度並不是隨意的,而是要求整個卷積核近似高斯函數的圖像.高斯濾波器對於抑制服從正態分布的噪聲非常有效.
線性濾波器(linear filter) :輸出圖像上每個像素點的值都是由輸入圖像各像素點值加權求和的結果.其原始數據與濾波結果是一種算術運算,即用加減乘除等運算實現,如方框濾波(BoxFilter)、均值濾波(MeanBlur)、高斯濾波(GaussianBlur)等.由於線性濾波器是算術運算,有固定的模板,因此濾波器的轉移函數是可以確定並且是唯一的.
非線性濾波器(non-linear filter):非線性濾波的算子中包含了取絕對值、置零等非線性運算.其原始數據與濾波結果是一種邏輯關係,即用邏輯運算實現,如最大值濾波器、最小值濾波器、中值濾波（medianBlur）和雙邊濾波（bilateralFilter）等,是通過比較一定鄰域內的灰度值大小來實現的,沒有固定的模板,因而也就沒有特定的轉移函數（因為沒有模板作傅立葉變換）,另外,膨脹和腐蝕也是通過最大值、最小值濾波器實現的.
思考1:為什麼要進行 加權均值 ?一幅圖像中輪廓和邊緣這種變化比較劇烈的地方高頻信號(噪聲和噪聲周邊的像素相比,其灰度變化同樣比較劇烈),其他的部分(可反應整副圖像的強度)是低頻信號,所以經過低通濾波器後圖像的輪廓和邊緣信號會被濾掉一部分,直觀感受是圖片變模糊了,這是平滑.從另一方面講,加權是將一個像素用鄰域幾個像素的均值代替,這樣會使信號的變化變得平緩(拉大小的值,拉小大的值),信號變化劇烈的地方是高頻信號,所以,這樣就達到了低通的效果,而低通可以濾除高頻噪聲.因此加權的目的在於實現低通.
思考2:如何實現 加權均值 中的求平均目的? 濾波核前面的係數等於矩陣中所有元素之和的倒數
計算高斯濾波的濾波核之前先了解一下高斯函數.
高斯函數在圖像處理中出現的頻率相當高.
一維高斯函數其中,μ是x的均值,σ是x的方差.本質上,x,μ都是空間中的坐標,x是卷積核內任一點的坐標,μ是卷積核中心的坐標.
由於每次計算時,都是以當前計算點為原點(中心點就是原點),因此均值μ等於0.所以公式進一步簡化為:
其函數曲線如下圖所示:
可以看到其是一鐘形曲線,σ決定了分布的形狀或者說表徵了鐘的寬度,σ越小形狀越瘦高,σ越大越矮胖
二維高斯函數在計算機視覺中,圖像是二維的,因此引入下面的二維高斯函數.
二維高斯函數為X,Y兩個方向的一維高斯函數的乘積:
同樣,在圖像濾波中,一般情況下
其中,(x,y)為點的坐標,在圖像處理中,其為整數.
其函數曲線如下圖所示:
從函數曲線可以看到,在整個定義域內,高斯函數值都是非負的,
要想得到一個高斯核,可以對高斯函數進行離散化,得到的高斯函數值作為高斯核的元素.例如:要產生一個3x3的高斯核,以高斯核的中心位置為坐標原點進行離散取樣, 那麼高斯核中心的8鄰域坐標為:
將各個位置的坐標代入到高斯函數中,得到的值就是初步的高斯核:
假設σ=1,則上面的高斯核求解結果如下:
上面初步計算3x3高斯核的權重和為0.7794836.考慮這樣一個問題,假若某一鄰域內所有像素的灰度值為255,那麼通過該高斯核進行卷積之後,模板中心像素的灰度值為0.7794836×255=198.77 < 255,偏離了實際的灰度值,使得圖像亮度相比原圖偏暗,產生了誤差.因此有必要對該高斯核進行歸一化處理,保證權重和為1以保證圖像的均勻灰度區域不受影響.
歸一化的原因更加正確的表述應該為:對灰度級為常數的圖像區域,高斯模板的響應該為1.這樣才能確保灰度級為常數的圖像區域經過高斯模板處理之後,依舊為其本身,而不是被改變為其他灰度級.
將求得的高斯核中各權重值除以權重和進行歸一化後如下:
作為對比,在matlab中,z = fspecial('gaussian', [3 3], 1);得到的結果和上面結果一致.
上面歸一化後求解的高斯核為小數形式的.
除了小數形式,高斯核也可以寫為整數形式.只用對求得的初始高斯核進行取整處理,就能得到整數形式的高斯核.從高斯函數圖像可以看到:1)從曲線中心位置到兩端,高斯函數值是單調遞減的;2)曲線位於x軸上方,其值總是正的.因此高斯核中距離中心最遠的位置(如高斯核左上角)的元素總是最小的,如果進行取整,那麼左上角元素的最小值也要保證≥1.所以對求解整數形式的高斯核過程如下:
1). 高斯核左上角的值歸一化為1.即高斯核中所有元素均除以左上角的元素.在本例中,左上角的元素為0.058549833,歸一化後高斯核為:
2).取整(對於取整方式,未查找到具體的原則.有文章選擇了選擇向下取整的方式,也許有向上取整,或者四捨五入等方式,但是從不同取整方式的結果來看,四捨五入的取整方式更接近於小數形式的結果).
3).取整後,高斯核所有元素的和明顯大於1,為了保證圖像的均勻灰度區域不受影響,因此需要進行權重歸一化,給高斯核一係數,以保證高斯和所有元素的和與係數的乘積為1.顯然,該係數為高斯核和的倒數,該例中取整、歸一化後高斯核為:
總結高斯核求解過程圖示如下:
濾波計算過程和CNN中的單通道卷積過程一致,只是此處濾波時的步長為1,濾波核為單通道,而CNN中卷積步長可以指定為1,2等值.
有了高斯核,就可以對圖像進行高斯濾波了.具體濾波過程和其他濾波方法一致.以一個3x3數據為例,利用上面求解的高斯核進行濾波.
其中心點的高斯濾波過程如下:
圖中中心像素值為50,經過濾波後其值仍為50,這只是由於數字選擇的原因導致其湊巧在數值上相等而已.如果將中心像素值修改為52,其濾波後值也為50.
在邊界附近的點,其鄰域內沒有足夠的點進行濾波處理,因此在濾波之前需要進行上下、左右邊界擴充操作,各邊界擴充的寬度為濾波半徑大小.具體的邊界擴充類型可以參考OpenCV中BorderTypes指定的類型進行.
考慮此,對上述3x3輸入進行濾波的過程如下:
易計算得濾波前輸入像素的和為450,濾波後輸出矩陣各元素相加後的和仍為450,即圖像的亮度在經過該高斯核濾波前後沒有發生變化.
高斯濾波步驟按照上面的流程,可以歸納高斯濾波的步驟如下:
1.生成高斯核.根據選擇的濾波半徑以及標準差,確定高斯核鄰域內各點到中心點的距離,代入二維高斯函數,求解高斯核各位置值,然後進行歸一化等處理;
2.邊界擴充.根據高斯核大小對輸入圖像進行擴充邊界操作;
3.卷積.針對圖像矩陣中的某個元素P,把高斯核的中心和待處理元素P空間對齊,將高斯核中各權值和P的鄰域內元素相乘後相加,作為P的輸出值;
4.對圖像中的每個元素(不包含邊界擴充的元素)進行步驟3的操作,得到的結果就是高斯濾波的結果.
從上面的步驟可以看到,高斯濾波僅考慮了圖像中各元素空間位置之間的關係,沒有將像素值之間的關係考慮進去. 而之前講到的雙邊濾波就不僅考慮像素空間之間的關係,同時考慮像素值之間的關係,它是空間域和像素值域之間的綜合結果.並且空間域和像素值域的求解都是套用高斯函數公式得到的.只是f(x,y)中x,y,μ的含義不同:針對空間域,它和這裡的高斯核是一樣的,x,y,μ表示的是空間位置,而μ=0;針對值域,x,y表示每一點的像素值,μ表示中心的像素值.所以它不光考慮了像素在空間中位置遠近程度的影響,還考慮了像素亮度相近程度的影響.從這個角度來看,雙邊濾波就很好理解了.
高斯濾波實現根據上面的濾波步驟,可以比較方便的實現高斯濾波代碼.
高斯濾波標準差與窗口大小的換算如下圖所示為高斯函數的分布特點.一般3σ外的數值已接近於0,可忽略.數值分布在(μ—3σ,μ+3σ)中的概率為0.997,或者說以均值為中心,半徑為3σ的範圍內已經包含了高斯函數99.7%以上的信息。如果高斯核的尺寸大,那麼相應的標準差σ也更大;相對應,如果σ大,那麼意味著高斯核的尺寸也要相應的進行增加.因此當我們知道高斯核的大小後,可以推斷出標準差的值.
我們大概可以感受到高斯核的半徑大約為3σ的大小.據此,在OpenCV中,有下面公式:
推導標準差的取值.公式中0.3比較好理解,但0.8的來歷並不是很清楚.
當知道標準差而不知道高斯核尺寸時,同樣可以由上式進行反向求解,得到高斯核尺寸,公式如下:
公式為先取整,然後和1進行按位或運算,以保證高斯核尺寸為奇數.
OpenCV中該部分代碼如下:
代碼中半徑為σ的3或4倍.
4倍的來源不明,無資料說明
實現常規實現常規實現過程按照上面的公式計算即可.其中高斯核或標準差的計算按上節公式實現.
代碼如下:
#include <iostream>
#include <chrono>
#include <opencv2/opencv.hpp>
#include <numeric>
/*!
* 生成小數形式的二維高斯核
* @param ksize 高斯核大小,奇數.當ksize<=0,並且sigma>0時,自動根據sigma計算ksize
* @param sigma 高斯濾波標準差.當sigma<=0,並且ksize>0為奇數時,自動根據ksize計算sigma
* @return 小數形式的高斯核.CV_32FC1.尺寸:ksize * ksize
* @author liheng
*/
cv::Mat GetGaussianKernel( int ksize, float sigma)
{
//根據sigma自動計算ksize
if( ksize <= 0 && sigma > 0 )
ksize = cvRound((sigma*3)*2 + 1)|1;// | 1操作以保證ksize為奇數
//校驗ksize
CV_Assert( ksize > 0 && ksize % 2 == 1);
//未設定sigma值時,根據窗口大小進行計算
sigma = (sigma>0 ? sigma : 0.3f*((ksize-1)*0.5f - 1) + 0.8f);
auto sigma2_inverse = 1.0/(2*sigma*sigma);// 1/(2*sigma^2)
cv::Mat kernel(ksize, ksize,CV_32FC1);
const int center = ksize/2;//高斯核中心坐標
for(int row=0;row<ksize;++row)
{
int y=row-center; //以center為中心,第row行的 y坐標
auto y2 = y*y;//y^2
for(int col=0;col<ksize;++col)
{
int x=col-center; //以center為中心,第col列的x坐標
//float f_xy = 1.0/(2*CV_PI*sigma*sigma) * std::exp( -(x*x+y*y)/(2*sigma*sigma) );
float f_xy = std::exp(-(x*x+y2)*sigma2_inverse );//減少不必要計算;係數在後面歸一化時會被消去;
kernel.at<float>(row,col) = f_xy;
}
}
//小數形式高斯核
auto sum = cv::sum(kernel)[0];
kernel /= sum;
////整數形式高斯核
////左上角係數歸一化為1
//kernel /= kernel.at<float>(0,0);
////向下取整
//kernel.convertTo(kernel,CV_32SC1);//向上取整?
//kernel.convertTo(kernel,CV_32FC1);//係數未計算
return kernel;
}
/*!
* 利用卷積核對輸入圖像進行二維卷積(濾波)
* @param src 輸入.CV_8UC1
* @param dst 輸出.CV_8UC1
* @param kernel 卷積核(相關核).CV_32FC1. kernel.size()需要為奇數
* @note 該函數在OpenCV中的對應函數:cv::filter2D(...)
* @author liheng
*/
void Filter2D(const cv::Mat& src,cv::Mat& dst,const cv::Mat& kernel)
{
CV_Assert( kernel.cols % 2 == 1 && kernel.rows % 2 == 1);
//邊界填充
int border_lr = kernel.cols/2;//左右填充寬度
int border_tb = kernel.rows/2;//上下填充寬度
cv::Mat borderedSrc;
cv::copyMakeBorder(src,borderedSrc,border_tb,border_tb,border_lr,border_lr,cv::BORDER_REFLECT);
dst.create(src.size(),src.type());
//對圖像中各元素(不包含填充元素)分別進行濾波
int rows = borderedSrc.rows;
int cols = borderedSrc.cols;
for(int row=border_tb; row<rows-border_tb; ++row)
{
for(int col=border_lr; col<cols-border_lr; ++col)
{
//遍歷卷積核鄰域求加權和作為該點輸出值
float sum = 0;
for(int y=-border_tb; y<=border_tb; ++y)//鄰域高度:2*border_tb+1
{
for(int x=-border_lr;x<=border_lr; ++x)//鄰域寬度:2*border_lr+1
sum += kernel.at<float>(y+border_tb,x+border_lr)*borderedSrc.at<uchar>(row+y,col+x);//權重與元素相乘
}
//輸出
dst.at<uchar>(row-border_tb,col-border_lr) = cv::saturate_cast<uchar>(sum);
}
}
}
/*!
* 自定義實現高斯濾波.原始版本,未採用分離卷積進行加速
* @param src 輸入.CV_8UC1
* @param dst 輸出.CV_8UC1
* @param ksize 高斯核大小,奇數.當ksize<=0,並且sigma>0時,自動根據sigma計算ksize
* @param sigma 高斯濾波標準差.當sigma<=0,並且ksize>0為奇數時,自動根據ksize計算sigma
* @author liheng
*/
void MyGaussianFilter(const cv::Mat& src,cv::Mat& dst,int ksize,float sigma)
{
//生成二維高斯核
cv::Mat kernel = GetGaussianKernel(ksize,sigma);
//濾波
Filter2D(src,dst,kernel);
}
int main(int argc, char *argv[])
{
cv::Mat kernel = GetGaussianKernel(3,1);
std::cout<<"二維高斯核:"<<std::endl;
std::cout<<kernel<<std::endl;
cv::Mat data = (cv::Mat_<uchar>(3,3)<<
10,20,30,
40,50,60,
70,80,90);
//
cv::Mat dst;
MyGaussianFilter(data,dst,3,1);
std::cout<<"自定義實現結果::"<<std::endl;
std::cout<<dst<<std::endl;
cv::Mat dst2;
cv::filter2D(data,dst2,CV_8UC1,kernel,cv::Point(-1,-1),0,cv::BORDER_REFLECT);
std::cout<<"filter2D結果::"<<std::endl;
std::cout<<dst2<<std::endl;
cv::GaussianBlur(data,data,cv::Size(3,3),1,1,cv::BORDER_REFLECT);
std::cout<<"cv::GaussianBlur結果::"<<std::endl;
std::cout<<data<<std::endl;
//輸出結果如下:
/**
二維高斯核:
[0.075113609, 0.1238414, 0.075113609;
0.1238414, 0.20417996, 0.1238414;
0.075113609, 0.1238414, 0.075113609]
自定義實現結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
filter2D結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
cv::GaussianBlur結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
**/
return 0;
}可以參考上面圖片中流程以輔助理解.
上面的濾波過程,其循環次數為row*col*(ksize*ksize),其中row,col為圖像的尺寸.其時間複雜度為
分離實現高斯濾波從二維高斯函數的公式可以看到,其為兩個一維高斯函數的乘積.因此二維高斯核具有可分離性,可以將高斯核分為兩個一維高斯核的乘積,如下圖所示:
同樣,利用高斯函數進行卷積(高斯濾波)的過程也具有可分離性. 其示意圖如下:
證明過程如下.
高斯濾波函數如上面所示的兩者之間為 離散二維卷積.下面公式即為離散二維卷積的定義.
記
可視為圖像I與f(y)之間的卷積運算.
記
在上面推導過程中,f(y),f(x)並不代表具體的行列順序,f(x)可以表示行方向，也可以代表列方向,f(y)同理.因此公式(3)可整理為:
從而:高斯函數進行卷積(高斯濾波)的過程同樣具有可分離性得證.
上面結論可以進行推廣.如果卷積核kernel是可分離的,並且
需要注意的是,上面的卷積核kernel為一維水平方向和一維垂直方向共2個卷積核的全卷積,這2個卷積核是滿足交換律的,但是對於多個卷積核的全卷積,其並不滿足交換律.對於均值濾波,其卷積核也是可分離的,如下:
用於邊緣檢測的Sobel卷積核也是可分離的,如下:
CNN中非對稱卷積(Asymmetric Convolution)/空間可分離卷積(Spatial Separable Convolution) 的思想和這裡有點相似.
按照公式的推導結論,可對原始高斯濾波進行優化,改寫為可分離卷積的形式.
為方便,上一小節代碼同樣拷貝到該代碼中.#include <iostream>
#include <chrono>
#include <opencv2/opencv.hpp>
#include <numeric>
/*!
* 生成小數形式的二維高斯核
* @param ksize 高斯核大小,奇數.當ksize<=0,並且sigma>0時,自動根據sigma計算ksize
* @param sigma 高斯濾波標準差.當sigma<=0,並且ksize>0為奇數時,自動根據ksize計算sigma
* @return 小數形式的高斯核.CV_32FC1.尺寸:ksize * ksize
* @author liheng
*/
cv::Mat GetGaussianKernel( int ksize, float sigma)
{
//根據sigma自動計算ksize
if( ksize <= 0 && sigma > 0 )
ksize = cvRound((sigma*3)*2 + 1)|1;// | 1操作以保證ksize為奇數
//校驗ksize
CV_Assert( ksize > 0 && ksize % 2 == 1);
//未設定sigma值時,根據窗口大小進行計算
sigma = (sigma>0 ? sigma : 0.3f*((ksize-1)*0.5f - 1) + 0.8f);
auto sigma2_inverse = 1.0/(2*sigma*sigma);// 1/(2*sigma^2)
cv::Mat kernel(ksize, ksize,CV_32FC1);
const int center = ksize/2;//高斯核中心坐標
for(int row=0;row<ksize;++row)
{
int y=row-center; //以center為中心,第row行的 y坐標
auto y2 = y*y;//y^2
for(int col=0;col<ksize;++col)
{
int x=col-center; //以center為中心,第col列的x坐標
//float f_xy = 1.0/(2*CV_PI*sigma*sigma) * std::exp( -(x*x+y*y)/(2*sigma*sigma) );
float f_xy = std::exp(-(x*x+y2)*sigma2_inverse );//減少不必要計算;係數在後面歸一化時會被消去;
kernel.at<float>(row,col) = f_xy;
}
}
//小數形式高斯核
auto sum = cv::sum(kernel)[0];
kernel /= sum;
////整數形式高斯核
////左上角係數歸一化為1
//kernel /= kernel.at<float>(0,0);
////向下取整
//kernel.convertTo(kernel,CV_32SC1);//向上取整?
//kernel.convertTo(kernel,CV_32FC1);//係數未計算
return kernel;
}
/*!
* 利用卷積核對輸入圖像進行卷積(濾波)
* @param src 輸入.CV_8UC1
* @param dst 輸出.CV_8UC1
* @param kernel 卷積核(相關核).CV_32FC1. kernel.size()需要為奇數
* @note 該函數在OpenCV中的對應函數:cv::filter2D(...)
* @author liheng
*/
void Filter2D(const cv::Mat& src,cv::Mat& dst,const cv::Mat& kernel)
{
CV_Assert( kernel.cols % 2 == 1 && kernel.rows % 2 == 1);
//邊界填充
int border_lr = kernel.cols/2;//左右填充寬度
int border_tb = kernel.rows/2;//上下填充寬度
cv::Mat borderedSrc;
cv::copyMakeBorder(src,borderedSrc,border_tb,border_tb,border_lr,border_lr,cv::BORDER_REFLECT);
dst.create(src.size(),src.type());
//對圖像中各元素(不包含填充元素)分別進行濾波
int rows = borderedSrc.rows;
int cols = borderedSrc.cols;
for(int row=border_tb; row<rows-border_tb; ++row)
{
for(int col=border_lr; col<cols-border_lr; ++col)
{
//遍歷卷積核鄰域求加權和作為該點輸出值
float sum = 0;
for(int y=-border_tb; y<=border_tb; ++y)//鄰域高度:2*border_tb+1
{
for(int x=-border_lr;x<=border_lr; ++x)//鄰域寬度:2*border_lr+1
sum += kernel.at<float>(y+border_tb,x+border_lr)*borderedSrc.at<uchar>(row+y,col+x);//權重與元素相乘
}
//輸出
dst.at<uchar>(row-border_tb,col-border_lr) = cv::saturate_cast<uchar>(sum);
}
}
}
/*!
* 自定義實現高斯濾波.原始版本,未採用分離卷積進行加速
* @param src 輸入.CV_8UC1
* @param dst 輸出.CV_8UC1
* @param ksize 高斯核大小,奇數.當ksize<=0,並且sigma>0時,自動根據sigma計算ksize
* @param sigma 高斯濾波標準差.當sigma<=0,並且ksize>0為奇數時,自動根據ksize計算sigma
* @author liheng
*/
void MyGaussianFilter(const cv::Mat& src,cv::Mat& dst,int ksize,float sigma)
{
//生成二維高斯核
cv::Mat kernel = GetGaussianKernel(ksize,sigma);
//濾波
Filter2D(src,dst,kernel);
}
//
/*!
* 生成小數形式的 "一維" 高斯核
* @param ksize 高斯核大小,奇數.當ksize<=0,並且sigma>0時,自動根據sigma計算ksize
* @param sigma 高斯濾波標準差.當sigma<=0,並且ksize>0為奇數時,自動根據ksize計算sigma
* @return 小數形式的高斯核.CV_32FC1.尺寸: 1 * ksize
* @author liheng
*/
cv::Mat GetSingleDimGaussianKernel( int ksize, double sigma)
{
//根據sigma自動計算ksize
if( ksize <= 0 && sigma > 0 )
ksize = cvRound((sigma*3)*2 + 1)|1;// | 1操作以保證ksize為奇數
//校驗ksize
CV_Assert( ksize > 0 && ksize % 2 == 1);
//未設定sigma值時,根據窗口大小進行計算
sigma = (sigma>0 ? sigma : 0.3f*((ksize-1)*0.5f - 1) + 0.8f);
auto sigma2_inverse = 1.0/(2*sigma*sigma);// 1/(2*sigma^2)
cv::Mat kernel(1, ksize,CV_32FC1);
const int center = ksize/2;//高斯核中心坐標
for(int col=0;col<ksize;++col)
{
int x=col-center; //以center為中心,第col列的 x坐標
//float f_x = 1.0/(std::sqrt(2*CV_PI)*sigma) * std::exp( -(x*x)/(2*sigma*sigma) );
float f_x = std::exp( -(x*x)*sigma2_inverse );//減少不必要計算;係數在後面歸一化時會被消去;
kernel.at<float>(0,col) = f_x;
}
auto sum = cv::sum(kernel)[0];
kernel /= sum;
return kernel;
}
/*!
* 自定義實現高斯濾波.改進版.利用高斯函數的可分離性以加速卷積過程
* @param src 輸入.CV_8UC1
* @param dst 輸出.CV_8UC1
* @param ksize 高斯核大小,奇數.當ksize<=0,並且sigma>0時,自動根據sigma計算ksize
* @param sigma 高斯濾波標準差.當sigma<=0,並且ksize>0為奇數時,自動根據ksize計算sigma
* @author liheng
*/
void separateGaussianFilter(const cv::Mat& src,cv::Mat& dst,int ksize,float sigma)
{
//獲取水平和豎直高斯核
cv::Mat kernelX = GetSingleDimGaussianKernel(ksize,sigma);
cv::Mat kernelY = kernelX.t();
cv::Mat Z;//中間變量.存儲水平方向高斯濾波後結果.
//水平方向高斯濾波
Filter2D(src,Z,kernelX);
//對水平方向濾波後結果進行豎直方向高斯濾波
Filter2D(Z,dst,kernelY);
}
int main(int argc, char *argv[])
{
cv::Mat kernel = GetGaussianKernel(3,1);
std::cout<<"二維高斯核:"<<std::endl;
std::cout<<kernel<<std::endl;
cv::Mat data = (cv::Mat_<uchar>(3,3)<<
10,20,30,
40,50,60,
70,80,90);
//
cv::Mat dst;
MyGaussianFilter(data,dst,3,1);
std::cout<<"自定義實現結果::"<<std::endl;
std::cout<<dst<<std::endl;
cv::Mat dst1;
separateGaussianFilter(data,dst1,3,1);
std::cout<<"自定義分離實現結果::"<<std::endl;
std::cout<<dst1<<std::endl;
cv::Mat dst2;
cv::filter2D(data,dst2,CV_8UC1,kernel,cv::Point(-1,-1),0,cv::BORDER_REFLECT);
std::cout<<"filter2D結果::"<<std::endl;
std::cout<<dst2<<std::endl;
cv::GaussianBlur(data,data,cv::Size(3,3),1,1,cv::BORDER_REFLECT);
std::cout<<"cv::GaussianBlur結果::"<<std::endl;
std::cout<<data<<std::endl;
//輸出結果如下:
/**
二維高斯核:
[0.075113609, 0.1238414, 0.075113609;
0.1238414, 0.20417996, 0.1238414;
0.075113609, 0.1238414, 0.075113609]
自定義實現結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
自定義分離實現結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
filter2D結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
cv::GaussianBlur結果::
[ 21, 28, 35;
43, 50, 57;
65, 72, 79]
**/
return 0;
}從代碼可以看到,其實現原理還是比較簡單的.首先得到一維的高斯核,在卷積（濾波）的過程中,保持行不變,列變化,在水平方向上做卷積運算;接著在上述得到的結果上,保持列不邊,行變化,在豎直方向上做卷積運算. 這樣分解開來,算法的時間複雜度為
總結高斯函數性質高斯函數具有五個重要的性質,這些性質使得它在早期圖像處理中特別有用．這些性質表明,高斯平滑濾波器無論在空間域還是在頻率域都是十分有效的低通濾波器,且在實際圖像處理中得到了工程人員的有效使用．高斯函數具有五個十分重要的性質,它們是:
1).二維高斯函數具有旋轉對稱性,即濾波器在各個方向上的平滑程度是相同的．一般來說,一幅圖像的邊緣方向是事先不知道的,因此,在濾波前是無法確定一個方向上比另一方向上需要更多的平滑．旋轉對稱性意味著高斯平滑濾波器在後續邊緣檢測中不會偏向任一方向．
2).高斯函數是單值函數．這表明,高斯濾波器用像素鄰域的加權均值來代替該點的像素值,而每一鄰域像素點權值是隨該點與中心點的距離單調增減的．這一性質是很重要的,因為邊緣是一種圖像局部特徵,如果平滑運算對離算子中心很遠的像素點仍然有很大作用,則平滑運算會使圖像失真．
3).高斯函數的傅立葉變換頻譜是單瓣的．如下圖所示,這一性質是高斯函數傅立葉變換等於高斯函數本身這一事實的直接推論．圖像常被不希望的高頻信號所汙染(噪聲和細紋理)．而所希望的圖像特徵（如邊緣）,既含有低頻分量,又含有高頻分量．高斯函數傅立葉變換的單瓣意味著平滑圖像不會被不需要的高頻信號所汙染,同時保留了大部分所需信號．高斯函數傅立葉變換等於高斯函數本身的證明可自行查找資料
一維高斯函數頻譜4).高斯濾波器寬度(決定著平滑程度)是由參數σ表徵的,而且σ和平滑程度的關係是非常簡單的．σ越大,高斯濾波器的頻帶就越寬,平滑程度就越好．通過調節平滑程度參數σ,可在圖像特徵過分模糊(過平滑)與平滑圖像中由於噪聲和細紋理所引起的過多的不希望突變量(欠平滑)之間取得折衷．
高斯濾波應用
5).由於高斯函數的可分離性,較大尺寸的高斯濾波器可以得以有效地實現．二維高斯函數卷積可以分兩步來進行,首先將圖像與一維高斯函數進行卷積,然後將卷積結果與方向垂直的相同一維高斯函數卷積．因此,二維高斯濾波的計算量隨濾波模板寬度成線性增長而不是成平方增長．高斯濾波後圖像被平滑的程度取決於標準差.它的輸出是領域像素的加權平均,同時離中心越近的像素權重越高.因此,相對於均值濾波(mean filter)它的平滑效果更柔和,而且邊緣保留的也更好.
高斯濾波被用作為平滑濾波器的本質原因是因為它是一個低通濾波器(其頻譜圖如上節所示).而且,大部份基於卷積平滑濾波器都是低通濾波器.
參考資料1.[圖像濾波之高斯濾波介紹]https://www.cnblogs.com/qiqibaby/p/5289977.html
附錄高斯函數及頻譜繪圖代碼
2.[圖像處理基礎(4):高斯濾波器詳解]https://www.cnblogs.com/wangguchangqing/p/6407717.html
3.[簡單易懂的高斯濾波]https://www.jianshu.com/p/73e6ccbd8f3f
4.[線性濾波、非線性濾波區別]https://blog.csdn.net/zqx951102/article/details/82967360
5.[如何確定高斯濾波的標準差和窗口大小]https://www.cnblogs.com/shine-lee/p/9671253.html
6.[二維圖像處理中的可分離卷積核]https://blog.csdn.net/jgj123321/article/details/95348795
7.[可分離卷積的運算量比較與分析（一維、二維卷積）]https://blog.csdn.net/Ocean_waver/article/details/106887264
8.[關於高斯濾波的一些理解]https://blog.csdn.net/lz0499/article/details/54015150繪製高斯函數曲線代碼:
一維高斯函數:from matplotlib import pyplot as plt
import numpy as np
def gaussian(x, mu, sig):
a = 1 / (np.sqrt(2 * np.pi) * sig)
return a*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
x_values = np.linspace(-11, 11, 120)
for mu, sig in [(0, 1), (0, 2), (0, 3)]:
plt.plot(x_values, gaussian(x_values, mu, sig),label='σ={}'.format(sig))
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.grid(linestyle='--')
plt.show()二維高斯函數:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
x,y = np.meshgrid(np.arange(-5,5,0.1),np.arange(-5,5,0.1))
sigma = 1
z = 1/(2 * np.pi * (sigma**2)) * np.exp(-(x**2+y**2)/(2 * sigma**2))
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap='rainbow',alpha = 0.9)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
# plt.title("二維高斯分布")
plt.show()繪製高斯傅立葉變換頻譜:
# Reference:https://www.cnblogs.com/jingsupo/p/9989559.html
import numpy as np#導入一個數據處理模塊
import pylab as pl#導入一個繪圖模塊,matplotlib下的模塊
def gaussian(x, mu, sig):
a = 1 / (np.sqrt(2 * np.pi) * sig)
return a*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
sampling_rate = 8000#採樣頻率為8KHz
fft_size = 8000 #FFT處理的取樣長度
t = np.arange(0, 4.0, 1.0/sampling_rate)#np.arange(起點,終點,間隔)產生4s長的取樣時間
x = gaussian(t,0,1)
xs = x[:fft_size]# 從波形數據中取樣fft_size個點進行運算
xf = np.fft.rfft(xs)# 利用np.fft.rfft()進行FFT計算,rfft()是為了更方便對實數信號進行變換,由公式可知/fft_size為了正確顯示波形能量
xf = xf / len(xf)
# rfft函數的返回值是N/2+1個複數,分別表示從0(Hz)到sampling_rate/2(Hz)的分.
#於是可以通過下面的np.linspace計算出返回值中每個下標對應的真正的頻率:
freqs = np.linspace(0, sampling_rate/2, fft_size//2+1)#頻率
xfp = np.abs(xf)#幅值
#繪圖顯示結果
pl.figure(figsize=(8,4))
pl.subplot(211)
pl.plot(t, x)
pl.xlim(0,)
pl.ylim(0,)
pl.xlabel('t')
pl.title('Gaussias\'s WaveForm And Freq')
pl.subplot(212)
freqs = freqs[:100]
xfp = xfp[:100]
_range = np.max(xfp) - np.min(xfp)
xfp=(xfp - np.min(xfp)) / _range
pl.plot(freqs, xfp)
pl.xlim(0,)
pl.ylim(0,)
pl.xlabel('Freq/Hz')
pl.ylabel('Amplitude')
pl.subplots_adjust(hspace=0.4)
pl.show()