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

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

1. 基本原理

非局部均值降噪算法（Non-Local Means）是空间降噪算法的一种，和中值滤波、高斯滤波这些局部滤波算法不同的是，非局部均值降噪算法是一种全局的算法，思路是利用整幅图像中相似像素的灰度值来代替当前像素的灰度值$${\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$之间的权重；$B(p,r)$为噪声图像中，以像素$p$为中心，宽为$2r+1$的区域，$C(p)$为权重归一化系数，计算公式为：$$C(p)=\sum _{q\in B(p,r)}w(p,q)$$

公式很好理解，中间比较重要的就是权重如何设计，权重需要描述两个像素之间的相似度，而这个相似度通常是通过这两个像素邻域像素间的欧拉距离来描述：$$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$的区域，在这个基础上，添加指数核函数来计算权值：$$w(p,q)=e^{-\frac{\max \left(d^2-2{\sigma }^2,0,0\right)}{h^2}}$$其中，$\sigma$和$h$是我们人为设定的参数，以上就完成了非局部均值降噪算法的理论介绍。

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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

第二份代码

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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

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

添加高斯噪声后：

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

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

3. 结论

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

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