本文主要介紹了高斯濾波器的原理及其實現過程算法
高斯濾波器是一種線性濾波器,可以有效的抑制噪聲,平滑圖像。其做用原理和均值濾波器相似,都是取濾波器窗口內的像素的均值做爲輸出。其窗口模板的係數和均值濾波器不一樣,均值濾波器的模板係數都是相同的爲1;而高斯濾波器的模板係數,則隨着距離模板中心的增大而係數減少。因此,高斯濾波器相比於均值濾波器對圖像個模糊程度較小。數組
既然名稱爲高斯濾波器,那麼其和高斯分佈(正態分佈)是有必定的關係的。一個二維的高斯函數以下:
\[ h(x,y) = e ^ {- \frac{x^2 + y^2}{2\sigma ^ 2}} \]
其中\((x,y)\)爲點座標,在圖像處理中可認爲是整數;\(\sigma\)是標準差。要想獲得一個高斯濾波器的模板,能夠對高斯函數進行離散化,獲得的高斯函數值做爲模板的係數。例如:要產生一個\(3 \times 3\)的高斯濾波器模板,以模板的中心位置爲座標原點進行取樣。模板在各個位置的座標,以下所示(x軸水平向右,y軸豎直向下)
函數
這樣,將各個位置的座標帶入到高斯函數中,獲得的值就是模板的係數。
對於窗口模板的大小爲 \((2k + 1) \times (2k + 1)\),模板中各個元素值的計算公式以下:
\[ H_{i,j} = \frac{1}{2\pi \sigma ^ 2}e ^{-\frac{(i - k - 1)^2 + (j - k - 1)^2}{2 \sigma ^ 2}} \]
這樣計算出來的模板有兩種形式:小數和整數。es5
知道模板生成的原理,實現起來也就不困難了spa
void generateGaussianTemplate(double window[][11], int ksize, double sigma) { static const double pi = 3.1415926; int center = ksize / 2; // 模板的中心位置,也就是座標的原點 double x2, y2; for (int i = 0; i < ksize; i++) { x2 = pow(i - center, 2); for (int j = 0; j < ksize; j++) { y2 = pow(j - center, 2); double g = exp(-(x2 + y2) / (2 * sigma * sigma)); g /= 2 * pi * sigma; window[i][j] = g; } } double k = 1 / window[0][0]; // 將左上角的係數歸一化爲1 for (int i = 0; i < ksize; i++) { for (int j = 0; j < ksize; j++) { window[i][j] *= k; } } }
須要一個二維數組,存放生成的係數(這裏假設模板的最大尺寸不會超過11);第二個參數是模板的大小(不要超過11);第三個參數就比較重要了,是高斯分佈的標準差。
生成的過程,首先根據模板的大小,找到模板的中心位置ksize/2
。 而後就是遍歷,根據高斯分佈的函數,計算模板中每一個係數的值。
須要注意的是,最後歸一化的過程,使用模板左上角的係數的倒數做爲歸一化的係數(左上角的係數值被歸一化爲1),模板中的每一個係數都乘以該值(左上角係數的倒數),而後將獲得的值取整,就獲得了整數型的高斯濾波器模板。
下面截圖生成的是,大小爲\(3 \times 3,\sigma = 0.8\)的模板
code
對上述解結果取整後獲得以下模板:
\[ \frac{1}{16}\left[ \begin{array}{c} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{array} \right] \]
這個模板就比較熟悉了,其就是根據\(\sigma = 0.8\)的高斯函數生成的模板。blog
至於小數形式的生成也比較簡單,去掉歸一化的過程,而且在求解過程後,模板的每一個係數要除以全部係數的和。具體代碼以下:it
void generateGaussianTemplate(double window[][11], int ksize, double sigma) { static const double pi = 3.1415926; int center = ksize / 2; // 模板的中心位置,也就是座標的原點 double x2, y2; double sum = 0; for (int i = 0; i < ksize; i++) { x2 = pow(i - center, 2); for (int j = 0; j < ksize; j++) { y2 = pow(j - center, 2); double g = exp(-(x2 + y2) / (2 * sigma * sigma)); g /= 2 * pi * sigma; sum += g; window[i][j] = g; } } //double k = 1 / window[0][0]; // 將左上角的係數歸一化爲1 for (int i = 0; i < ksize; i++) { for (int j = 0; j < ksize; j++) { window[i][j] /= sum; } } }
\(3 \times 3,\sigma = 0.8\)的小數型模板。圖像處理
經過上述的實現過程,不難發現,高斯濾波器模板的生成最重要的參數就是高斯分佈的標準差\(\sigma\)。標準差表明着數據的離散程度,若是\(\sigma\)較小,那麼生成的模板的中心繫數較大,而周圍的係數較小,這樣對圖像的平滑效果就不是很明顯;反之,\(\sigma\)較大,則生成的模板的各個係數相差就不是很大,比較相似均值模板,對圖像的平滑效果比較明顯。opencv
來看下一維高斯分佈的機率分佈密度圖:
橫軸表示可能得取值x,豎軸表示機率分佈密度F(x),那麼不難理解這樣一個曲線與x軸圍成的圖形面積爲1。\(\sigma\)(標準差)決定了這個圖形的寬度,能夠得出這樣的結論:\(\sigma\)越大,則圖形越寬,尖峯越小,圖形較爲平緩;\(\sigma\)越小,則圖形越窄,越集中,中間部分也就越尖,圖形變化比較劇烈。這其實很好理解,若是sigma也就是標準差越大,則表示該密度分佈必定比較分散,因爲面積爲1,因而尖峯部分減少,寬度越寬(分佈越分散);同理,當\(\sigma\)越小時,說明密度分佈較爲集中,因而尖峯越尖,寬度越窄!
因而能夠獲得以下結論:
\(\sigma\)越大,分佈越分散,各部分比重差異不大,因而生成的模板各元素值差異不大,相似於平均模板;
\(\sigma\)越小,分佈越集中,中間部分所佔比重遠遠高於其餘部分,反映到高斯模板上就是中心元素值遠遠大於其餘元素值,因而天然而然就至關於中間值得點運算。
在生成高斯模板好,其簡單的實現和其餘的空間濾波器沒有區別,具體代碼以下:
void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma) { CV_Assert(src.channels() || src.channels() == 3); // 只處理單通道或者三通道圖像 const static double pi = 3.1415926; // 根據窗口大小和sigma生成高斯濾波器模板 // 申請一個二維數組,存放生成的高斯模板矩陣 double **templateMatrix = new double*[ksize]; for (int i = 0; i < ksize; i++) templateMatrix[i] = new double[ksize]; int origin = ksize / 2; // 以模板的中心爲原點 double x2, y2; double sum = 0; for (int i = 0; i < ksize; i++) { x2 = pow(i - origin, 2); for (int j = 0; j < ksize; j++) { y2 = pow(j - origin, 2); // 高斯函數前的常數能夠不用計算,會在歸一化的過程當中給消去 double g = exp(-(x2 + y2) / (2 * sigma * sigma)); sum += g; templateMatrix[i][j] = g; } } for (int i = 0; i < ksize; i++) { for (int j = 0; j < ksize; j++) { templateMatrix[i][j] /= sum; cout << templateMatrix[i][j] << " "; } cout << endl; } // 將模板應用到圖像中 int border = ksize / 2; copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT); int channels = dst.channels(); int rows = dst.rows - border; int cols = dst.cols - border; for (int i = border; i < rows; i++) { for (int j = border; j < cols; j++) { double sum[3] = { 0 }; for (int a = -border; a <= border; a++) { for (int b = -border; b <= border; b++) { if (channels == 1) { sum[0] += templateMatrix[border + a][border + b] * dst.at<uchar>(i + a, j + b); } else if (channels == 3) { Vec3b rgb = dst.at<Vec3b>(i + a, j + b); auto k = templateMatrix[border + a][border + b]; sum[0] += k * rgb[0]; sum[1] += k * rgb[1]; sum[2] += k * rgb[2]; } } } for (int k = 0; k < channels; k++) { if (sum[k] < 0) sum[k] = 0; else if (sum[k] > 255) sum[k] = 255; } if (channels == 1) dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]); else if (channels == 3) { Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) }; dst.at<Vec3b>(i, j) = rgb; } } } // 釋放模板數組 for (int i = 0; i < ksize; i++) delete[] templateMatrix[i]; delete[] templateMatrix; }
只處理單通道或者三通道圖像,模板生成後,其濾波(卷積過程)就比較簡單了。不過,這樣的高斯濾波過程,其循環運算次數爲\(m \times n \times ksize^2\),其中m,n爲圖像的尺寸;ksize爲高斯濾波器的尺寸。這樣其時間複雜度爲\(O(ksize^2)\),隨濾波器的模板的尺寸呈平方增加,當高斯濾波器的尺寸較大時,其運算效率是極低的。爲了,提升濾波的運算速度,能夠將二維的高斯濾波過程分解開來。
因爲高斯函數的可分離性,尺寸較大的高斯濾波器能夠分紅兩步進行:首先將圖像在水平(豎直)方向與一維高斯函數進行卷積;而後將卷積後的結果在豎直(水平)方向使用相同的一維高斯函數獲得的模板進行卷積運算。具體實現代碼以下:
// 分離的計算 void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma) { CV_Assert(src.channels()==1 || src.channels() == 3); // 只處理單通道或者三通道圖像 // 生成一維的高斯濾波模板 double *matrix = new double[ksize]; double sum = 0; int origin = ksize / 2; for (int i = 0; i < ksize; i++) { // 高斯函數前的常數能夠不用計算,會在歸一化的過程當中給消去 double g = exp(-(i - origin) * (i - origin) / (2 * sigma * sigma)); sum += g; matrix[i] = g; } // 歸一化 for (int i = 0; i < ksize; i++) matrix[i] /= sum; // 將模板應用到圖像中 int border = ksize / 2; copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT); int channels = dst.channels(); int rows = dst.rows - border; int cols = dst.cols - border; // 水平方向 for (int i = border; i < rows; i++) { for (int j = border; j < cols; j++) { double sum[3] = { 0 }; for (int k = -border; k <= border; k++) { if (channels == 1) { sum[0] += matrix[border + k] * dst.at<uchar>(i, j + k); // 行不變,列變化;先作水平方向的卷積 } else if (channels == 3) { Vec3b rgb = dst.at<Vec3b>(i, j + k); sum[0] += matrix[border + k] * rgb[0]; sum[1] += matrix[border + k] * rgb[1]; sum[2] += matrix[border + k] * rgb[2]; } } for (int k = 0; k < channels; k++) { if (sum[k] < 0) sum[k] = 0; else if (sum[k] > 255) sum[k] = 255; } if (channels == 1) dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]); else if (channels == 3) { Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) }; dst.at<Vec3b>(i, j) = rgb; } } } // 豎直方向 for (int i = border; i < rows; i++) { for (int j = border; j < cols; j++) { double sum[3] = { 0 }; for (int k = -border; k <= border; k++) { if (channels == 1) { sum[0] += matrix[border + k] * dst.at<uchar>(i + k, j); // 列不變,行變化;豎直方向的卷積 } else if (channels == 3) { Vec3b rgb = dst.at<Vec3b>(i + k, j); sum[0] += matrix[border + k] * rgb[0]; sum[1] += matrix[border + k] * rgb[1]; sum[2] += matrix[border + k] * rgb[2]; } } for (int k = 0; k < channels; k++) { if (sum[k] < 0) sum[k] = 0; else if (sum[k] > 255) sum[k] = 255; } if (channels == 1) dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]); else if (channels == 3) { Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) }; dst.at<Vec3b>(i, j) = rgb; } } } delete[] matrix; }
代碼沒有重構較長,不過其實現原理是比較簡單的。首先獲得一維高斯函數的模板,在卷積(濾波)的過程當中,保持行不變,列變化,在水平方向上作卷積運算;接着在上述獲得的結果上,保持列不邊,行變化,在豎直方向上作卷積運算。 這樣分解開來,算法的時間複雜度爲\(O(ksize)\),運算量和濾波器的模板尺寸呈線性增加。
在OpenCV也有對高斯濾波器的封裝GaussianBlur
,其聲明以下:
CV_EXPORTS_W void GaussianBlur( InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY = 0, int borderType = BORDER_DEFAULT );
二維高斯函數的標準差在x和y方向上應該分別有一個標準差,在上面的代碼中一直設其在x和y方向的標準是相等的,在OpenCV中的高斯濾波器中,能夠在x和y方向上設置不一樣的標準差。
下圖是本身實現的高斯濾波器和OpenCV中的GaussianBlur
的結果對比
上圖是\(5\times5,\sigma = 0.8\)的高斯濾波器,能夠看出兩個實現獲得的結果沒有很大的區別。
高斯濾波器是一種線性平滑濾波器,其濾波器的模板是對二維高斯函數離散獲得。因爲高斯模板的中心值最大,四周逐漸減少,其濾波後的結果相對於均值濾波器來講更好。
高斯濾波器最重要的參數就是高斯分佈的標準差\(\sigma\),標準差和高斯濾波器的平滑能力有很大的能力,\(\sigma\)越大,高斯濾波器的頻帶就較寬,對圖像的平滑程度就越好。經過調節\(\sigma\)參數,能夠平衡對圖像的噪聲的抑制和對圖像的模糊。