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

1. 基本原理

非局部均值降噪算法（Non-Local Means）是空间降噪算法的一种，和中值滤波、高斯滤波这些局部滤波算法不同的是，非局部均值降噪算法是一种全局的算法，思路是利用整幅图像中相似像素的灰度值来代替当前像素的灰度值u^i(p)=1C(p)∑q∈B(p,r)ui(q)w(p,q)\hat{u}_{i}(p)=\frac{1}{C(p)} \sum_{q \in B(p, r)} u_{i}(q) w(p, q)u^i(p)=C(p)1q∈B(p,r)∑ui(q)w(p,q)其中，ui(q)u_{i}(q)ui(q)是噪声图像像素qqq的灰度值；u^i(p)\hat{u}_{i}(p)u^i(p)是降噪后图像像素ppp的灰度值；w(p,q)w(p, q)w(p,q)是像素ppp和qqq之间的权重；B(p,r)B(p, r)B(p,r)为噪声图像中，以像素ppp为中心，宽为2r+12r+12r+1的区域，C(p)C(p)C(p)为权重归一化系数，计算公式为：C(p)=∑q∈B(p,r)w(p,q)C(p)=\sum_{q \in B(p, r)} w(p, q)C(p)=q∈B(p,r)∑w(p,q)

公式很好理解，中间比较重要的就是权重如何设计，权重需要描述两个像素之间的相似度，而这个相似度通常是通过这两个像素邻域像素间的欧拉距离来描述：d2(B(p,f),B(q,f))=13(2f+1)2∑i=13∑j∈B(0,Ω)(ui(p+j)−ui(q+j))2d^{2}(B(p, f), B(q, f))=\frac{1}{3(2 f+1)^{2}} \sum_{i=1}^{3} \sum_{j \in B(0, \Omega)}\left(u_{i}(p+j)-u_{i}(q+j)\right)^{2}d2(B(p,f),B(q,f))=3(2f+1)21i=1∑3j∈B(0,Ω)∑(ui(p+j)−ui(q+j))2其中，333次求和是对于彩色图而言的，B(p,f)B(p, f)B(p,f)为噪声图像中，以像素ppp为中心，宽为2f+12f+12f+1的区域，在这个基础上，添加指数核函数来计算权值：w(p,q)=e−max⁡(d2−2σ2,0,0)h2w(p, q)=e^{-\frac{\max \left(d^{2}-2 \sigma^{2}, 0,0\right)}{h^{2}}}w(p,q)=e−h2max(d2−2σ2,0,0)其中，σ\sigmaσ和hhh是我们人为设定的参数，以上就完成了非局部均值降噪算法的理论介绍。

2. C++代码实现

这里我基于OpenCV完成了两份代码，其中第一份是我根据上面公式自己实现，比较容易理解，但是运行速度较慢。因为太慢了，所以我尝试写了第二份代码。第二份是参考他人的代码基于Mat指针实现的，因为是指针操作，所以运行速度会相对较快。

第一份代码

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);//构建边界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;//遍历搜索区域每一个像素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);//构建边界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);//搜索区域中心指针uchar* pUpLeft = pad.data + pad.step * i + j;//搜索区域左上角指针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;
}

下面是输出结果
首先是原图：

添加高斯噪声后：

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

可以看出这个效果还是非常感人的

3. 结论

可以看到非局部均值降噪算法的效果视觉上是非常可以接受的，但是它的缺点是时间复杂度太高，可以清楚看到程序里面有六个for循环，当我把搜索尺寸设为21，模板尺寸设为7是：
第一份代码运行时间：135.034秒
第二份代码运行时间：29.6425秒
OpenCV库函数运行时间：0.512942秒
通过对比可以发现，Mat通过at函数操作速度要慢于指针操作，而我写的指针操作的代码要远慢于OpenCV库函数，OpenCV中的库函数应该是采用多线程操作。
2014年有一篇paper讲的是通过积分图的方式优化NLM，我对其进行的C++代码实现，单线程下速度能提高到7秒左右，同时我也简单读了下OpenCV的原理，同样是使用了空间换时间的思想，不过相对积分图的方式会更加巧妙，值得进一步学习。

此外，这里我写一个各种算法的总结目录图像降噪算法——图像降噪算法总结，对图像降噪算法感兴趣的同学欢迎参考