圖像降噪算法——非局部均值降噪算法

• 圖像降噪算法——非局部均值降噪算法
• • 1. 基本原理
  • 2. C++代碼實現(xiàn)
  • 3. 結(jié)論

圖像降噪算法——非局部均值降噪算法

1. 基本原理

非局部均值降噪算法（Non-Local Means）是空間降噪算法的一種，和中值濾波、高斯濾波這些局部濾波算法不同的是，非局部均值降噪算法是一種全局的算法，思路是利用整幅圖像中相似像素的灰度值來代替當(dāng)前像素的灰度值$${\hat{u}}_i(p)=\frac{1}{C(p)}\sum _{q\in B(p,r)}{u}_i(q)w(p,q)$$其中，$u_i(q)$是噪聲圖像像素$q$的灰度值；${\hat{u}}_i(p)$是降噪后圖像像素$p$的灰度值；$w(p,q)$是像素$p$和$q$之間的權(quán)重；$B(p,r)$為噪聲圖像中，以像素$p$為中心，寬為$2r+1$的區(qū)域，$C(p)$為權(quán)重歸一化系數(shù)，計算公式為：$$C(p)=\sum _{q\in B(p,r)}w(p,q)$$

公式很好理解，中間比較重要的就是權(quán)重如何設(shè)計，權(quán)重需要描述兩個像素之間的相似度，而這個相似度通常是通過這兩個像素鄰域像素間的歐拉距離來描述：$$d^2(B(p,f),B(q,f))=\frac{1}{3(2f+1{)}^2}\sum _{i=1}^3\sum _{j\in B(0,\Omega )}{\left(u_i(p+j)?u_i(q+j)\right)}^2$$其中，$3$次求和是對于彩色圖而言的，$B(p,f)$為噪聲圖像中，以像素$p$為中心，寬為$2f+1$的區(qū)域，在這個基礎(chǔ)上，添加指數(shù)核函數(shù)來計算權(quán)值：$$w(p,q)=e^{?\frac{\max \left(d^2?2{\sigma }^2,0,0\right)}{h^2}}$$其中，$\sigma$和$h$是我們?nèi)藶樵O(shè)定的參數(shù)，以上就完成了非局部均值降噪算法的理論介紹。

2. C++代碼實現(xiàn)

這里我基于OpenCV完成了兩份代碼，其中第一份是我根據(jù)上面公式自己實現(xiàn)，比較容易理解，但是運行速度較慢。因為太慢了，所以我嘗試寫了第二份代碼。第二份是參考他人的代碼基于Mat指針實現(xiàn)的，因為是指針操作，所以運行速度會相對較快。

第一份代碼

Mat Denoise::NonLocalMeansFilter(const Mat &src, int searchWindowSize, int templateWindowSize, double sigma, double h) {Mat dst, pad;dst = Mat::zeros(src.rows, src.cols, CV_8UC1);//構(gòu)建邊界int padSize = (searchWindowSize+templateWindowSize)/2;copyMakeBorder(src, pad, padSize, padSize, padSize, padSize, cv::BORDER_CONSTANT);int tN = templateWindowSize*templateWindowSize;int sN = searchWindowSize*searchWindowSize;vector<double> gaussian(256*256, 0);for(int i = 0; i<256*256; i++){double g = exp(-max(i-2.0*sigma*sigma, 0.0))/(h*h);gaussian[i] = g;if(g<0.001)break;}//遍歷圖像上每一個像素for(int i = 0; i<src.rows; i++){for(int j = 0; j<src.cols; j++){cout<<i<<" "<<j<<endl;//遍歷搜索區(qū)域每一個像素int pX = i+searchWindowSize/2;int pY = j+searchWindowSize/2;vector<vector<double>> weight(searchWindowSize, vector<double>(searchWindowSize, 0));double weightSum = 0;for(int m = searchWindowSize-1; m>=0; m--){for(int n = searchWindowSize-1; n>=0; n--){int qX = i+m;int qY = j+n;int w = 0;for(int x = templateWindowSize-1; x>=0; x--){for(int y = templateWindowSize-1; y>=0; y--){w += pow(pad.at<uchar>(pX+x, pY+y) - pad.at<uchar>(qX+x, qY+y), 2);}}weight[m][n] = gaussian[(int)(w/tN)];weightSum += weight[m][n];}}dst.at<uchar>(i,j) = 0;double sum = 0;for(int m = 0; m<searchWindowSize; m++){for(int n = 0; n<searchWindowSize; n++){sum += pad.at<uchar>(i+templateWindowSize/2+m, j+templateWindowSize/2+n)*weight[m][n];}}dst.at<uchar>(i,j) = (uchar)(sum/weightSum);}}return dst; }

第二份代碼

Mat Denoise::NonLocalMeansFilter2(const Mat &src, int searchWindowSize, int templateWindowSize, double sigma, double h) {Mat dst, pad;dst = Mat::zeros(src.rows, src.cols, CV_8UC1);//構(gòu)建邊界int padSize = (searchWindowSize+templateWindowSize)/2;copyMakeBorder(src, pad, padSize, padSize, padSize, padSize, cv::BORDER_CONSTANT);int tN = templateWindowSize*templateWindowSize;int sN = searchWindowSize*searchWindowSize;int tR = templateWindowSize/2;int sR = searchWindowSize/2;vector<double> gaussian(256*256, 0);for(int i = 0; i<256*256; i++){double g = exp(-max(i-2.0*sigma*sigma, 0.0))/(h*h);gaussian[i] = g;if(g<0.001)break;}double* pGaussian = &gaussian[0];const int searchWindowStep = (int)pad.step - searchWindowSize;const int templateWindowStep = (int)pad.step - templateWindowSize;for(int i = 0; i < src.rows; i++){uchar* pDst = dst.ptr(i);for(int j = 0; j < src.cols; j++){cout<<i<<" "<<j<<endl;int *pVariance = new int[sN];double *pWeight = new double[sN];int cnt = sN-1;double weightSum = 0;uchar* pCenter = pad.data + pad.step * (sR + i) + (sR + j);//搜索區(qū)域中心指針uchar* pUpLeft = pad.data + pad.step * i + j;//搜索區(qū)域左上角指針for(int m = searchWindowSize; m>0; m--){uchar* pDownLeft = pUpLeft + pad.step * m;for(int n = searchWindowSize; n>0; n--){uchar* pC = pCenter;uchar* pD = pDownLeft + n;int w = 0;for(int k = templateWindowSize; k>0; k--){for(int l = templateWindowSize; l>0; l--){w += (*pC - *pD)*(*pC - *pD);pC++;pD++;}pC += templateWindowStep;pD += templateWindowStep;}w = (int)(w/tN);pVariance[cnt--] = w;weightSum += pGaussian[w];}}for(int m = 0; m<sN; m++){pWeight[m] = pGaussian[pVariance[m]]/weightSum;}double tmp = 0.0;uchar* pOrigin = pad.data + pad.step * (tR + i) + (tR + j);for(int m = searchWindowSize, cnt = 0; m>0; m--){for(int n = searchWindowSize; n>0; n--){tmp += *(pOrigin++) * pWeight[cnt++];}pOrigin += searchWindowStep;}*(pDst++) = (uchar)tmp;delete pWeight;delete pVariance;}}return dst; }

下面是輸出結(jié)果
首先是原圖：

添加高斯噪聲后：

通過非局部均值降噪算法降噪效果：

可以看出這個效果還是非常感人的

3. 結(jié)論

可以看到非局部均值降噪算法的效果視覺上是非常可以接受的，但是它的缺點是時間復(fù)雜度太高，可以清楚看到程序里面有六個for循環(huán)，當(dāng)我把搜索尺寸設(shè)為21，模板尺寸設(shè)為7是：
第一份代碼運行時間：135.034秒
第二份代碼運行時間：29.6425秒
OpenCV庫函數(shù)運行時間：0.512942秒
通過對比可以發(fā)現(xiàn)，Mat通過at函數(shù)操作速度要慢于指針操作，而我寫的指針操作的代碼要遠(yuǎn)慢于OpenCV庫函數(shù)，OpenCV中的庫函數(shù)應(yīng)該是采用多線程操作。

2014年有一篇paper講的是通過積分圖的方式優(yōu)化NLM，我對其進行的C++代碼實現(xiàn)，單線程下速度能提高到7秒左右，同時我也簡單讀了下OpenCV的原理，同樣是使用了空間換時間的思想，不過相對積分圖的方式會更加巧妙，值得進一步學(xué)習(xí)。

此外，這里我寫一個各種算法的總結(jié)目錄圖像降噪算法——圖像降噪算法總結(jié)，對圖像降噪算法感興趣的同學(xué)歡迎參考

總結(jié)

以上是生活随笔為你收集整理的图像降噪算法——非局部均值降噪算法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。