非局部均值降噪算法（NLM）原理及实现

一、概述

在日常生活中，最常见的 CT 图像噪声是高斯白噪声。目前，针对高斯白噪声的处理方法，主要有空间域中的以平滑为基本思想的均值滤波、高斯滤波、局部滤波等，此外还有频率域中的维纳滤波和小波阈值收缩等。局部强度统计特征是衡量区域内像素间的平均相似性，但这一特征难以准确辨别边缘与其邻近点之间的差异，导致了滤波结果中边缘信息的模糊。小波阈值收缩算法虽能够很好地估计信号的噪声，但是阈值和阈值函数的选择上存在不确定性，会造成信号的失真。
2005 年Buades 等提出了去除图像加性噪声的非局部均值滤波算法（Non-Local Means，NLM），该算法是利用图像非局部结构的相似性来去除噪声，恢复图像的主要几何结构。论文原文：A non-local algorithm for image denoising，还有一篇2011年的论文：Non-Local Means Denoising。

二、算法原理

如下所示为一般均值滤波的原理，其中权重仅通过中心像素和其邻域像素进行计算。当滤波半径过大的时候，图像的细节信息很难被保留，即使像双边滤波引入了像素间的灰度信息也会出现细节丢失严重的情况。
如下所示为双边滤波的公式及示意图：

非局部均值滤波是考虑到图像的非局部统计自相似性质，利用图像包含的大量重复结构进行去噪。通过计算非局部像素之间的相似度确定权值，然后加权平均来恢复图像。NLM 的主要特点 ：将相似像素定义为具有相同邻域模式的像素，对像素周围固定大小的窗口内的信息进行比较，而不只是比较图像单个像素的信息，因此得到的相似性更加可靠。NLM基本原理如下 ：

通过观察NLM的计算公式不难发现，其相对于常见的局部均值滤波，其区别主要体现在权重的计算上。在计算权重的相似度时，其通过像素间的邻域进行计算。

注意：通常非局部均值的滤波核都设置得比较大，如上图中所示的滤波核大小为9x9，因此NLM能够获得更佳的平滑度，同时更好的保留细节。

三、算法流程

四、MATLAB实现

无优化的实现：

function [output]=nlm_filter(input,t,f,h)

% 输入: 待平滑的图像
% t: 搜索窗口半径
% f: 相似性窗口半径
% h: 平滑参数
% NLmeans(ima,5,2,sigma);

% 图像大小
[m n]=size(input);
% 输出
output=zeros(m,n);
input2 = padarray(input,[f+t f+t],'symmetric');%边界作对称处理

% 高斯核
ksize = f*2 + 1;
kernel = fspecial('gaussian',[ksize,ksize],1);
kernel = kernel / sum(sum(kernel));
kernel(:) = 1;
h=h*h;

for i=1:m
    for j=1:n
        i1 = i+ f+t;%原始图像的像素位置 （中心像素）
        j1 = j+ f+t;
        
        W1= input2(i1-f:i1+f , j1-f:j1+f);%小窗口
        
        wmax=0;
        average=0;
        sweight=0;
        
        %rmin = max(i1-t,f+1);
        %rmax = min(i1+t,m+f);
        %smin = max(j1-t,f+1);
        %smax = min(j1+t,n+f);
        rmin=i1-t;
        rmax=i1+t;
        smin=j1-t;
        smax=j1+t;
        
        for r=rmin:1:rmax %大窗口
            for s=smin:1:smax
                
                if(r==i1 && s==j1)
                    continue;
                end
                
                W2= input2(r-f:r+f , s-f:s+f);    %大搜索窗口中的小相似性窗口
                d = sum(sum(kernel.*(W1-W2).*(W1-W2)));
                w=exp(-d/h); %权重
                
                if w>wmax
                    wmax=w;   %求最大权重
                end
                
                sweight = sweight + w;  %大窗口中的权重和
                average = average + w*input2(r,s);
            end
        end
        
        average = average + wmax*input2(i1,j1);
        sweight = sweight + wmax;
        
        if sweight > 0
            output(i,j) = average / sweight;
        else
            output(i,j) = input(i,j);
        end
    end
end

使用积分图加速（这部分还没仔细看其实现原理）：

clc;
clear;
close all;

src=imread('lena.jpg');
src=rgb2gray(src);
src=double(src);
figure,imshow(src,[]),title('src image')

[m,n]=size(src);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(src,[ds+Ds,ds+Ds],'symmetric','both');% 扩展图像，便于处理

sumimage=zeros(m,n);
sumweight=zeros(m,n);
maxweight=zeros(m,n);
image=PaddedImg(1+Ds:Ds+m+ds,1+Ds:Ds+n+ds);
[M,N]=size(image);
for r=-Ds:Ds
    for s=-Ds:Ds
        
        %跳过当前点偏移
        if(r==0&&s==0)
            continue;
        end
        
        %求得差值积分图
        wimage=PaddedImg(1+Ds+r:Ds+m+ds+r,1+Ds+s:Ds+n+ds+s);
        diff=image-wimage;
        diff=diff.^2;
        J=cumsum(diff,1);
        J=cumsum(J,2);
        
        %计算距离
        distance=J(M-m+1:M,N-n+1:N)+J(1:m,1:n)-J(M-m+1:M,1:n)-J(1:m,N-n+1:N);
        distance=distance/((2*ds+1).^2);
        
        %计算权重并获得单个偏移下的加权图像
        weight=exp(-distance./(h*h));
        sumimage=sumimage+weight.*wimage(ds+1:ds+m,ds+1:ds+n);
        sumweight=sumweight+weight;
        maxweight=max(maxweight,weight);
    end
end
sumimage=sumimage+maxweight.*image(ds+1:ds+m,ds+1:ds+n);
sumweight=sumweight+maxweight;
dst=sumimage./sumweight;

figure,imshow(dst,[]),title('dst');

降噪前后的效果比较：

五、C++实现

//计算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];
    }
  }
}

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

//sigma越大越平滑
//halfKernelSize邻域窗
//halfSearchSize搜索窗
void NL_mean(Mat src, Mat &dst, double sigma, 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 sigma_square = sigma*sigma;

  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 / sigma_square);
          p += boardSrc_p[i + sc] * w;
          sumw += w;
        }
      }
      dst_p[i - boardSize] = saturate_cast<uchar>(p / sumw);
    }
  }
}

参考文献

[1] 非局部均值滤波算法（NL-means）
[2] https://developer.aliyun.com/article/1156618
[3] opencv源码修改与使用：fastNlMeansDenoisingMulti()对CV_16U的支持