圖像噪聲，通常指圖像中除了成像物體之外的其它資訊，比如斑點和顆粒，這些額外的錯誤資訊幹擾了成像物體的顯示，影響成像品質，是以往往需要通過圖像濾波(也稱為圖像去噪)來消除這些噪點。常見的圖像濾波算法有均值濾波、高斯濾波、中值濾波、雙邊濾波、非局部均值濾波，以及近幾年火熱的基于深度學習的圖像濾波等。本章節将詳細講解均值濾波算法的原理，以及C++實作和優化。首先膜拜一下那些寫Opencv代碼的大佬們，他們寫的代碼不僅穩定性良好，運作效率也超級高，很多時候我們費盡心思寫了一個相同的算法，發現性能與Opencv的接口函數相比還是差了許多，是以會有一丢丢的心理落差，但是抱着學習的态度，追趕大佬的腳步，精益求精，相信我們自己也是可以的！均值濾波，也就是計算每一個像素點周圍像素點(包括該點)的平均值，作為該像素點濾波之後的值，通常取以該像素點為中心的矩形視窗内的所有像素點來計算平均值，矩形視窗的大小一般為3*3，5*5，9*9，...，(2n+1)* (2n+1)。視窗越大，濾波效果越好，但是圖像也變得更加模糊，是以需要根據實際情況設定矩形視窗的大小。比如3*3視窗的均值濾波如下圖所示，點(x,y)的濾波值由其周圍9個點(包括其自身)計算平均值得到。
上圖中點(x,y)的濾波值用公式表示為：
對于 (2n+1)* (2n+1)視窗，點(x,y)的平均濾波值可根據如下公式計算：
為了解決圖像邊緣像素點取不到完整矩形視窗的問題，通常先把圖像的上、下邊界擴充n行，左右邊界擴充n列，實際計算時，隻計算圖像原有像素點的視窗平均值。比如當矩形視窗為3*3，則n的值為1，這種情況下擴充邊界的示意圖如下圖所示：
根據以上原理，基于Opencv和C++的均值濾波實作代碼如下：
void blur_mean(Mat src, Mat &dst, int winsize){ if(winsize&1) //如果視窗的邊不是奇數，則加1使其為奇數，因為隻有視窗的邊為奇數的時候目前像素點才是視窗的中心點 { winsize += 1; } const int winsize_2 = winsize/2; //winsize_2 就是上述公式中的n const float winsize_num = winsize*winsize; //(2n+1)*(2n+1) Mat src_board; //調用Opencv的copyMakeBorder函數擴充邊界 copyMakeBorder(src, src_board, winsize_2, winsize_2, winsize_2, winsize_2, BORDER_REFLECT); const int row = src_board.rows; //行 const int col = src_board.cols; Mat dst_tmp(src.size(), CV_8UC1); //列 for(int i = winsize_2; i < row-winsize_2; i++) //行循環，隻計算圖像的原有行 { for(int j = winsize_2; j //列循環，隻計算圖像的原有列 { float sum = 0.0; //計算每一個像素點周圍矩形區域内所有像素點的累加和 for(int y = 0; y < winsize; y++) { for(int x = 0; x < winsize; x++) { sum += src_board.ptr(i-winsize_2+y)[j-winsize_2+x]; } } //求得累加和之後再求視窗像素的平均值。作為目前像素點的濾波值 dst_tmp.ptr(i-winsize_2)[j-winsize_2] = (uchar)(sum/winsize_num + 0.5); } } dst_tmp.copyTo(dst);}
運作以上代碼，對一幀1024*1024的圖像進行均值濾波，得到的結果如下圖所示。可以看到，随着視窗尺寸的增加，濾波圖像變得更加模糊，整體耗時也大幅度增加，由此可以判斷在濾波的計算過程中，主要耗時點為計算視窗内像素點的累加和。
原圖
3*3視窗，耗時17 ms
11*11視窗，耗時147.8 ms
21*21視窗，耗時535.4 ms
前面兩篇文章中主要講解了積分圖的計算，根據積分圖的特點，可以使用積分圖來簡化矩形視窗内像素累加和的計算，如下圖所示，藍色矩形區域内點的像素累加和，可以使用其四個頂點A、B、C、D的積分值來計算。
假設 A、B、C、D的積分值分别為I(A)、 I(B )、 I(C )、 I(D )，那麼藍色區域的像素累加和可以按照下式計算：
将上述積分圖計算矩形區域内像素和的原理應用于均值濾波中，可以大大簡化運算。對于每一個像素點，其計算濾波值的計算量由原本的(2n+1)*(2n+1)次加法，簡化為2次加法和1次減法。使用積分圖優化加速的均值濾波代碼如下：
void blur_mean_integral(Mat src, Mat &dst, int winsize){ if(winsize&1) { winsize += 1; } const int winsize_2 = winsize/2; const float winsize_num = 1.0/(winsize*winsize); Mat src_board; copyMakeBorder(src, src_board, winsize_2, winsize_2, winsize_2, winsize_2, BORDER_REFLECT); Mat integral; Integal_row(src_board, integral); //計算積分圖 const int row = src_board.rows; const int col = src_board.cols; dst.create(src.size(), CV_8UC1); for(int i = winsize_2; i < row-winsize_2; i++) { for(int j = winsize_2; j < col-winsize_2; j++) { //使用矩形區域四個頂點的積分值來計算區域内的x像素累加和 float sum = integral.ptr<float>(i+winsize_2)[j+winsize_2] - integral.ptr<float>(i-winsize_2)[j+winsize_2] - integral.ptr<float>(i+winsize_2)[j-winsize_2] + integral.ptr<float>(i-winsize_2)[j-winsize_2]; //得到累加和之後再計算平均值 dst.ptr(i-winsize_2)[j-winsize_2] = (uchar)(sum*winsize_num + 0.5); } }}
運作以上代碼，取視窗大小為21*21，同樣對一幀1024*1024的圖像進行均值濾波，耗時由535.4 ms減小到6.78 ms，可以看到，整體耗時大大減小。在使用積分圖的基礎上，再進行SSE指令的優化，可以進一步減小計算耗時：
void blur_mean_integral(Mat src, Mat &dst, int winsize){ if(winsize&1) { winsize += 1; } const int winsize_2 = winsize/2; const float winsize_num = 1.0/(winsize*winsize); Mat src_board; copyMakeBorder(src, src_board, winsize_2, winsize_2, winsize_2, winsize_2, BORDER_REFLECT); Mat integral; Integal_row(src_board, integral); const int row = src_board.rows; const int col = src_board.cols; dst.create(src.size(), CV_8UC1); for(int i = winsize_2; i < row-winsize_2; i++) { float *p1 = integral.ptr<float>(i-winsize_2); float *p2 = integral.ptr<float>(i+winsize_2); uchar *p3 = dst.ptr(i-winsize_2); int j = winsize_2; for(; j -4; j+= { /*float sum = p2[j+winsize_2] - p1[j+winsize_2] - p2[j-winsize_2] + p1[j-winsize_2]; p3[j-winsize_2] = (uchar)(sum*winsize_num + 0.5); sum = p2[j+winsize_2+1] - p1[j+winsize_2+1] - p2[j-winsize_2+1] + p1[j-winsize_2+1]; p3[j-winsize_2+1] = (uchar)(sum*winsize_num + 0.5); sum = p2[j+winsize_2+2] - p1[j+winsize_2+2] - p2[j-winsize_2+2] + p1[j-winsize_2+2]; p3[j-winsize_2+2] = (uchar)(sum*winsize_num + 0.5); sum = p2[j+winsize_2+3] - p1[j+winsize_2+3] - p2[j-winsize_2+3] + p1[j-winsize_2+3]; p3[j-winsize_2+3] = (uchar)(sum*winsize_num + 0.5);*/ //以下的SSE指令代碼對應上方的C++代碼 __m128 a3, a2, a1, a0; //将4個點的同一方向頂點的像素值加載到__m128變量中 //a3 : p2[j+winsize_2+3] p2[j+winsize_2+2] p2[j+winsize_2+1] p2[j+winsize_2] //a2 : p1[j+winsize_2+3] p1[j+winsize_2+2] p1[j+winsize_2+1] p1[j+winsize_2] //a1 : p2[j-winsize_2+3] p2[j-winsize_2+2] p2[j-winsize_2+1] p2[j-winsize_2] //a0 : p1[j-winsize_2+3] p1[j-winsize_2+2] p1[j-winsize_2+1] p1[j-winsize_2] a3 = _mm_loadu_ps(&p2[j+winsize_2]); a2 = _mm_loadu_ps(&p1[j+winsize_2]); a1 = _mm_loadu_ps(&p2[j-winsize_2]); a0 = _mm_loadu_ps(&p1[j-winsize_2]); //(a3-a2)+(a0-a1) __m128 sum = _mm_add_ps(_mm_sub_ps(a3, a2), _mm_sub_ps(a0, a1)); __m128 winnum = _mm_set1_ps(winsize_num); //winsize_num winsize_num winsize_num winsize_num sum = _mm_mul_ps(sum, winnum); //sum中有4個浮點數，winnum中也有4個浮點數，兩者對應位置的浮點數相乘 __m128i sum_i = _mm_cvtps_epi32(sum); //四舍五入取整 p3[j-winsize_2+3] = (uchar)sum_i.m128i_i32[3]; p3[j-winsize_2+2] = (uchar)sum_i.m128i_i32[2]; p3[j-winsize_2+1] = (uchar)sum_i.m128i_i32[1]; p3[j-winsize_2] = (uchar)sum_i.m128i_i32[0]; } for(; j < col-winsize_2; j++) { float sum = p2[j+winsize_2] - p1[j+winsize_2] - p2[j-winsize_2] + p1[j-winsize_2]; p3[j-winsize_2] = (uchar)(sum*winsize_num + 0.5); } }}
使用SSE指令優化之後，同樣對1024*1024的圖像，取視窗21*21進行濾波，耗時由 6.78 ms減少到3.98 ms。由此可知還是有一定優化效果的。