1. 非局部均值滤波

非局部均值滤波（Non-Local Means，NL-Means）是一种非线性的图像去噪算法。它基于图像中的像素具有相似结构这一假设，利用图像的全局信息来对图像进行去噪。

1.1 全局算法 VS 局部算法

非局部均值滤波在计算每个像素点的估计值时，会考虑图像中所有与该像素点具有相似邻域结构的像素点。因此，非局部均值滤波是一种全局算法。

那么相对于全局算法的局部算法是什么呢？局部算法是指仅利用图像局部信息进行处理的算法。例如其邻域窗口内的信息，来计算该像素点的估计值。常用的局部算法包括：

• 均值滤波
• 中值滤波
• 高斯滤波
• 双边滤波

局部算法的优点是计算量小，速度快。但其缺点是去噪效果有限，容易造成图像细节丢失。

1.2 均值滤波和非局部均值滤波的区别

均值滤波：对于图像中的每个像素点，其滤波值是其邻域窗口内所有像素点的平均值。

非局部均值滤波：该算法需要计算图像中所有像素与当前像素之间的相似性。首先需要定义一个大的搜索窗口和一个小的邻域窗口。搜索窗口用于搜索与当前像素点具有相似邻域结构的像素点，邻域窗口用于计算像素点之间的相似度。邻域窗口在搜索窗口中滑动，对于搜索窗口内的每个像素点，计算其与当前像素点的邻域窗口的相似度，并将其作为权重。相似度越大则权重越大。

非局部均值滤波的基本原理与均值滤波类似都要取平均值，但是非局部均值滤波在计算中加入了每一个点的权重。

非局部均值滤波是一种比均值滤波更先进的图像去噪方法，但其计算量也更大。

1.3 非局部均值滤波的原理

非局部均值滤波的公式如下：

其中，w(x,y) 是一个权重，表示在原始图像 v 中，像素 x 和像素 y 的相似度。是像素 x 的搜索窗口。是滤波后的图像。

非局部均值滤波的示意.png

常用的相似度度量方法包括：欧式距离、高斯相似度、L1 范数、L2 范数、MSE 等等。

衡量两个图像块的相似度最常用的方法是计算他们之间的欧氏距离：

其中：

• n(x) 是一个归一化的因子，是所有权重的和。对每个权重除以该因子后，使得权重满足和为1的条件。
• a 是高斯核的标准差。在求欧式距离的时候，不同位置的像素的权重是不一样的，距离块的中心越近，权重越大，距离中心越远，权重越小，权重服从高斯分布。实际计算中常常采用均匀分布的权重。
• h 是滤波系数。控制指数函数的衰减从而改变欧氏距离的权重，h >0 。

非局部均值滤波的复杂度跟图像的大小、搜索窗口的大小、相似度计算方法、权重计算方法密切相关。

2. 非局部均值滤波的实现

下面的例子，是在图像中添加斑点噪声，然后用非局部均值滤波消除噪声。

#include <opencv2/opencv.hpp>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <random>

using namespace std;
using namespace cv;

void addSpeckleNoise(Mat& image, double scale, Mat &dst) {
    dst = image.clone();
    RNG rng;

    dst.forEach<Pixel>([&](Pixel &p, const int * position) -> void {
        int row = position[0];
        int col = position[1];

        double random_value = rng.uniform(0.0, 1.0);
        double noise_intensity = random_value * scale;
        dst.at<Vec3b>(row, col) = dst.at<Vec3b>(row, col) + Vec3b(noise_intensity * 255, noise_intensity * 255, noise_intensity * 255);
    });
}

//NL-means 算法的实现
void nonlocalMeansFilter(Mat& src, Mat& dst, int templeteWindowSize, int searchWindowSize, double h, double sigma = 0.0){
    //邻域的大小不能超过搜索窗口的大小
    if (templeteWindowSize > searchWindowSize){
        cout << "searchWindowSize should be larger than templeteWindowSize" << endl;
        return;
    }

    if (dst.empty())
        dst = Mat::zeros(src.size(), src.type());

    const int tr = templeteWindowSize >> 1;//tr为邻域的中心位置
    const int sr = searchWindowSize >> 1;  //sr为搜索域的中心位置
    const int bb = sr + tr;//需增加的边界宽度
    const int D = searchWindowSize*searchWindowSize;//搜索域中的元素个数
    const int H = D / 2 + 1;//搜索域中的中心点位置
    const double div = 1.0 / (double)D;//均匀分布时，搜索域中的每个点的权重大小
    const int tD = templeteWindowSize*templeteWindowSize;//邻域中的元素个数
    const double tdiv = 1.0 / (double)(tD);//均匀分布时，搜索域中的每个点的权重大小

    //扩充边界
    Mat boardSrc;
    copyMakeBorder(src, boardSrc, bb, bb, bb, bb, cv::BORDER_DEFAULT);

    //weight computation;
    vector<double> weight(256 * 256 * src.channels());
    double* w = &weight[0];
    const double gauss_sd = (sigma == 0.0) ? h : sigma;//高斯标准差
    double gauss_color_coeff = -(1.0 / (double)(src.channels())) * (1.0 / (h*h));//高斯颜色系数
    int emax=0;

    //w[i]保存方差，即邻域平均欧氏距离对应的高斯加权权重，供后面计算出欧式距离后调用
    for (int i = 0; i < weight.size(); i++){
        double v = std::exp(max(i - 2.0*gauss_sd*gauss_sd, 0.0)*gauss_color_coeff);
        w[i] = v;
        if (v<0.001){
            emax = i;
            break;
        }
    }

    for (int i = emax; i < weight.size(); i++)
        w[i] = 0.0;

    int height = src.rows;
    int width = src.cols;

    if (src.channels() == 3){
        const int cstep = (int)boardSrc.step - templeteWindowSize * 3;
        const int csstep = (int)boardSrc.step - searchWindowSize * 3;
#pragma omp parallel for
        for (int j = 0; j<height; j++){
            uchar* d = dst.ptr(j);
            int* ww = new int[D];//D 为搜索域中的元素数量，ww用于记录搜索域每个点的邻域方差
            double* nw = new double[D];//根据方差大小高斯加权归一化后的权重
            for (int i = 0; i<width; i++){
                double tweight = 0.0;
                //search loop
                uchar* tprt = boardSrc.data + boardSrc.step * (sr + j) + 3 * (sr + i);
                uchar* sptr2 = boardSrc.data + boardSrc.step * j + 3 * i;
                for (int l = searchWindowSize, count = D - 1; l--;){
                    uchar* sptr = sptr2 + boardSrc.step * (l);
                    for (int k = searchWindowSize; k--;){
                        //templete loop
                        int e = 0;
                        uchar* t = tprt;
                        uchar* s = sptr + 3 * k;
                        for (int n = templeteWindowSize; n--;){
                            for (int m = templeteWindowSize; m--;){
                                // computing color L2 norm
                                e += (s[0] - t[0])*(s[0] - t[0]) + (s[1] - t[1])*(s[1] - t[1]) + (s[2] - t[2])*(s[2] - t[2]);//L2 norm
                                s += 3;
                                t += 3;
                            }
                            t += cstep;
                            s += cstep;
                        }
                        const int ediv = e*tdiv;
                        ww[count--] = ediv;
                        //get weighted Euclidean distance
                        tweight += w[ediv];
                    }
                }

                //weight normalization
                if (tweight == 0.0){
                    for (int z = 0; z<D; z++) nw[z] = 0;
                    nw[H] = 1;
                }else{
                    double itweight = 1.0 / (double)tweight;
                    for (int z = 0; z<D; z++) nw[z] = w[ww[z]] * itweight;
                }
                double r = 0.0, g = 0.0, b = 0.0;
                uchar* s = boardSrc.ptr(j + tr); s += 3 * (tr + i);
                for (int l = searchWindowSize, count = 0; l--;){
                    for (int k = searchWindowSize; k--;)
                    {
                        r += s[0] * nw[count];
                        g += s[1] * nw[count];
                        b += s[2] * nw[count++];
                        s += 3;
                    }
                    s += csstep;
                }
                d[0] = saturate_cast<uchar>(r);
                d[1] = saturate_cast<uchar>(g);
                d[2] = saturate_cast<uchar>(b);
                d += 3;
            }//i
            delete[] ww;
            delete[] nw;
        }//j
    } else if (src.channels() == 1){
        const int cstep = (int)boardSrc.step - templeteWindowSize;//在邻域比较时，从邻域的上一行末尾跳至下一行开头
        const int csstep = (int)boardSrc.step - searchWindowSize;//搜索域循环中，从搜索域的上一行末尾跳至下一行开头
#pragma omp parallel for
        for (int j = 0; j<height; j++){
            uchar* d = dst.ptr(j);
            int* ww = new int[D];
            double* nw = new double[D];
            for (int i = 0; i<width; i++){
                double tweight = 0.0;
                uchar* tprt = boardSrc.data + boardSrc.step * (sr + j) + (sr + i);
                uchar* sptr2 = boardSrc.data + boardSrc.step * j + i;
                for (int l = searchWindowSize, count = D - 1; l--;){
                    uchar* sptr = sptr2 + boardSrc.step * (l);
                    for (int k = searchWindowSize; k--;){
                        int e = 0;
                        uchar* t = tprt;
                        uchar* s = sptr + k;
                        for (int n = templeteWindowSize; n--;){
                            for (int m = templeteWindowSize; m--;){
                                e += (*s - *t)*(*s - *t);
                                s++;
                                t++;
                            }
                            t += cstep;
                            s += cstep;
                        }
                        const int ediv = e*tdiv;
                        ww[count--] = ediv;
                        tweight += w[ediv];
                    }
                }

                if (tweight == 0.0){
                    for (int z = 0; z<D; z++) nw[z] = 0;
                    nw[H] = 1;
                }else{
                    double itweight = 1.0 / (double)tweight;
                    for (int z = 0; z<D; z++) nw[z] = w[ww[z]] * itweight;
                }
                double v = 0.0;
                uchar* s = boardSrc.ptr(j + tr); s += (tr + i);
                for (int l = searchWindowSize, count = 0; l--;){
                    for (int k = searchWindowSize; k--;){
                        v += *(s++)*nw[count++];
                    }
                    s += csstep;
                }
                *(d++) = saturate_cast<uchar>(v);
            }//i
            delete[] ww;
            delete[] nw;
        }//j
    }
}

int main() {
    Mat src = imread(".../girl.jpg");

    imshow("src", src);

    Mat result;

    Mat dst4;
    addSpeckleNoise(src,0.5,dst4);
    imshow("addSpeckleNoise", dst4);

    nonlocalMeansFilter(dst4, result, 3, 15, 40,40);
    imshow("removeSpeckleNoise", result);

    waitKey(0);
    return 0;
}

斑点噪声和使用非局部均值滤波后的效果.png

OpenCV 提供了非局部均值滤波算法，并对其进行了加速。

• fastNlMeansDenoising()：对单个灰度图像进行去噪。
• fastNlMeansDenoisingColored() ：对彩色图像进行去噪。
• fastNlMeansDenoisingMulti() ：用于连续相关灰度图像的快速去噪（例如视频中的连续灰度帧）。
• fastNlMeansDenoisingColoredMulti() ：用于连续相关彩色图像的快速去噪（例如视频中的连续彩色帧）。

int main() {
    Mat src = imread(.../girl.jpg");

    imshow("src", src);

    Mat result;

    Mat dst4;
    addSpeckleNoise(src,0.5,dst4);
    imshow("addSpeckleNoise", dst4);

    fastNlMeansDenoisingColored(dst4,result,40,40);
    imshow("removeSpeckleNoise2", result);

    waitKey(0);
    return 0;
}

3. 总结

非局部均值滤波能够有效地去除图像中的各种噪声，包括高斯噪声、椒盐噪声、纹理噪声等。非局部均值滤波能够较好地保留图像的细节，对噪声的类型和分布不敏感，具有较强的鲁棒性。

当然，非局部均值滤波的缺点也很明显：计算量大，容易受到图像边缘的影响等等。

非局部均值滤波的计算量大、速度慢是可以通过减少搜索窗口大小、采用快速相似度计算方法、利用图像的稀疏性、并行化计算、利用硬件加速等方法来加速。