非局部均值滤波(NL-means)算法的原理与C++实现

常见的空间滤波算法有均值滤波、高斯滤波、中值滤波等。其中，均值滤波的核心思路是取每一个像素点邻域的矩形窗口，计算矩形窗口内所有像素点的像素平均值，作为该点滤波之后的像素值。高斯滤波与均值滤波类似，都是计算矩形窗口内所有像素点的像素值加权和，只不过其权重与均值滤波不一样，高斯滤波的权重服从二维正态分布，越靠近窗口中心点（也即当前滤波点），权重越大。

本文我们主要讲非局部均值（NL-means）滤波算法的原理与实现。其核心思路与高斯滤波很相似：计算矩形窗口内所有像素点的像素值加权和，权重服从高斯分布。区别在于：高斯滤波使用当前滤波点与矩形窗口内其它点的空间欧式距离来计算权重，距离越近权重越大；而非局部均值滤波则使用当前滤波点的邻域块与矩形窗口内其它点的邻域块的相似度来计算权重，相似度越大则权重越大。

高斯滤波使用的点与点的空间距离d如下图所示：

假设A的坐标为(x1,y1)，B的坐标为(x2,y2)，那么d可按照下式计算：

非局部均值滤波使用的点与点的邻域块相似度如下图所示：

要衡量两个邻域块的相似度，有多种指标，均方误差（MSE）是最常用的相似度衡量指标之一。非局部均值滤波算法就是使用MSE来计算两个邻域块的相似度。邻域块的行列相同，假设均为m行n列，那么MSE的计算如下式，其中A(i,j)为点A邻域块中的点(i,j)的像素值，B(i,j)为点B邻域块中的相同位置点的像素值：

由以上可知，非局部均值滤波算法有两个重要的参数：ksize和ssize。假设点A为当前待滤波点，邻域块分别取以点A和点B为中心的ksize*ksize区域，矩形窗口（也称为搜索窗口）取以点A为中心ssize*ssize大小。由于ksize和ssize都是奇数，通常传入它们的半值：half_ksize和half_ssize。所以，邻域块的大小为(2*half_ksize+1)*(2*half_ksize+1)，搜索窗口的大小为(2*half_ssize+1)*(2*half_ssize+1)。如下图所示：

点A的滤波值由搜索窗口内所有点的像素值加权平均得到：

上式中，w(A,B)为点A与搜索窗口中任意一点B的高斯权重，由两点各自邻域块的MSE相似度计算而得。如下式，w(A,B)需要除以搜索窗口内所有点的高斯系数之和，以实现归一化的目的。其中h也是一个重要的参数，h越大去噪效果越好，但是图像越模糊，反之h越小去噪效果越差，但去噪之后的失真度越小。

由于原图像中每一个待滤波的点都需要取完整的一个搜索窗口和多个邻域块，因此非局部均值滤波之前，首先要对图像进行边缘扩充，上、下、左、右分别扩充half_ssize+half_ksize的行、行、列、列：

好了，原理就讲到这里，下面上代码。

1. 查找表的计算

由于图像的像素值范围是确定的0~255，因此计算邻域块的MSE时可以提前计算好查找表，以节省一点时间。结合以下两个查找表可以快速计算两个数的差值平方。

//计算0~255的平方查找表
float table1[256];
static void cal_lookup_table1(void)
{
  for (int i = 0; i < 256; i++)
  {
    table1[i] = (float)(i*i);
  }
}


//计算两个0~255的数的绝对差值的查找表
uchar table2[256][256];
static void cal_lookup_table2(void)
{
  for (int i = 0; i < 256; i++)
  {
    for (int j = i; j < 256; j++)
    {
      table2[i][j] = abs(i - j);
      table2[j][i] = table2[i][j];
    }
  }
}

2. 计算两个邻域块的MSE代码

float MSE_block(Mat m1, Mat m2)
{
  float sum = 0.0;
  for (int j = 0; j < m1.rows; j++)
  {
    uchar *data1 = m1.ptr<uchar>(j);
    uchar *data2 = m2.ptr<uchar>(j);
    for (int i = 0; i < m1.cols; i++)
    {
      sum += table1[table2[data1[i]][data2[i]]];
    }
  }
  sum = sum / (m1.rows*m2.cols);
  return sum;
}

3. 主体函数

//h越大越平滑
//halfKernelSize小框
//halfSearchSize大框
void NL_mean(Mat src, Mat &dst, double h, int halfKernelSize, int halfSearchSize)
{
  Mat boardSrc;
  dst.create(src.rows, src.cols, CV_8UC1);
  int boardSize = halfKernelSize + halfSearchSize;
  copyMakeBorder(src, boardSrc, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT);   //边界扩展
  double h2 = h*h;


  int rows = src.rows;
  int cols = src.cols;


  cal_lookup_table1();
  cal_lookup_table2();


  for (int j = boardSize; j < boardSize + rows; j++)
  {
    uchar *dst_p = dst.ptr<uchar>(j - boardSize);
    for (int i = boardSize; i < boardSize + cols; i++)
    {
      Mat patchA = boardSrc(Range(j - halfKernelSize, j + halfKernelSize), Range(i - halfKernelSize, i + halfKernelSize));
      double w = 0;
      double p = 0;
      double sumw = 0;


      for (int sr = -halfSearchSize; sr <= halfSearchSize; sr++)   //在搜索框内滑动
      {
        uchar *boardSrc_p = boardSrc.ptr<uchar>(j + sr);
        for (int sc = -halfSearchSize; sc <= halfSearchSize; sc++)
        {
          Mat patchB = boardSrc(Range(j + sr - halfKernelSize, j + sr + halfKernelSize), Range(i + sc - halfKernelSize, i + sc + halfKernelSize));
          float d2 = MSE_block(patchA, patchB);


          w = exp(-d2 / h2);
          p += boardSrc_p[i + sc] * w;
          sumw += w;
        }
      }


      dst_p[i - boardSize] = saturate_cast<uchar>(p / sumw);


    }
  }


}

4. 测试代码

分别使用上述实现的代码，以及Opencv实现的均值滤波、高斯滤波、中值滤波函数，对有很大噪声的Lene图像进行滤波。

void nlmean_test(void)
{
  Mat img = imread("lena_noise.jpg", CV_LOAD_IMAGE_GRAYSCALE);
  
  Mat out;
  nlmean3(img, out, 20, 3, 15);   //NL-means
  
  Mat out1;
  blur(img, out1, Size(11, 11));   //均值滤波


  Mat out2;
  GaussianBlur(img, out2, Size(11, 11), 2.5);   //高斯滤波


  Mat out3;
  medianBlur(img, out3, 9);   //中值滤波


  imshow("img", img);
  imshow("非局部", out);
  imshow("均值", out1);
  imshow("高斯", out2);
  imshow("中值", out3);
  waitKey();
}

5. 滤波结果对比

噪声图

均值滤波结果

高斯滤波结果

中值滤波结果

非局部均值滤波结果

由上图可知，非局部均值滤波的效果在几种滤波算法中是最好的，但是非局部均值滤波的耗时也是最大的，Lena图像的大小为496*472，这个尺寸并不算大，但是却耗时约20秒。耗时太长很影响其实际应用，因此在接下来的篇章中我们将继续讲解该算法的加速优化，敬请期待~