详解非局部均值滤波原理以及用MATLAB源码实现

详解非局部均值滤波原理以及用MATLAB源码实现

---

序言

均值滤波、中值滤波、高斯滤波在滤除噪声的过程中，无可避免的使图像的边缘细节和纹理信息所被滤除。针对此问题，Buades[1]等人提出了非局部均值滤波，它使图像中的冗余信息得到了有效地利用，在滤除噪声的同时将图像的细节特征进行了最大程度的保留。

非局部均值滤波使图像中的冗余信息得到了有效地利用，在滤除噪声的同时将图像的细节特征进行了最大程度的保留。

非局部均值滤波利用了自然图像中的每个小块都存在关联性，与均值滤波是对邻域内的所有像素求和再平均的方法不同，它先在整幅图像中寻找相似的图像块，再根据图像块的相似度大小来赋予其不同的权值，以此来实现图像去噪。

原理介绍

1 借助图片

可用N（p）来表示以待处理像素点p为中心、大小为N×N的图像块；同理，N（q1）、N（q2）、N（q3）分别代表以q1、q2、q3为中心、大小为N×N的图像块。

由于N（q1）、N（q2）、N（q3）与N（p）的相似度均不相同，它们会通过与图像块N（p）的相似度比较后，根据其与N（p）各自的相似性被赋予不同的相似性权值，N（q1）、N（q2）、N（q3）的相似性权值分别可用w(p*,* q1)、w(p*,* q2)、w(p*,* q3)来表示。

值得说明的是，从上图可看出，N（q1）、N（q2）与N（p）具有相同的属性，因此相似性权值w(p, q1)、w(p, q2)会大一些；而N（q3）与N（p）之间存在较大差异，相似性权值w(p, q3)会比较小。

---

2 进一步用公式阐述非局部均值滤波的原理

假定一幅离散的噪声图像为g=g{g(i)|i∈Ω}，其大小为N×N，Ω代表图像邻域，i为像素索引，对噪声图像进行非局部均值滤波处理后可表达为：

这是一个加权平均的过程。上式中分母的存在是为了使其归一化，Ωi表示中心像素为i、大小为q×q的搜索区域；w(i,j)代表赋予噪声图像g(i)的权值，数学形式为：

h代表控制滤波平滑程度的滤波系数；N(i)和N(j)分别代表以像素中心点为i、像素中心点为j的图像块，图像块的大小都是p×p；d(i,j)表示图像块N(i)和N(j)之间的欧式加权距离，α表示高斯核的标准差，α>0。在非局部均值滤波算法中，搜索区域中将q的设置一般不会超过21，并且图像块中p的取值会比q的值小的多。

非局部均值滤波中的一些参数选择会直接影响着滤波效果的好坏，因此是至关重要的。传统的非局部均值滤波中，用来计算两个图像块之间的欧式距离一般采用高斯加权欧式距离，上面d(i,j)可进一步表达为：

上式中“⊗”用来表示两个矩阵之间的点乘；Gα是高斯矩阵。使用高斯加权欧式距离的一个好处是：噪声可能会影响两个图像块之间的相似性计算，采用高斯核函数卷积类似于先采用高斯滤波来滤除噪声，这样后期计算相似度会更准确。

• 滤波参数h大小说明

非局部均值滤波采用的是指数函数的减函数，那么滤波系数h就会对加权核函数的衰减作用有着决定性的作用。

如果h取值较大的话，求取的权值w(i,j)的值就会趋于一致使滤波作用会等效于均值滤波，去噪的效果由于边缘细节丢失会变得比较模糊；

若h取值较小，对噪声的滤除作用就不好从未导致图像中会有哦大量噪声存在，去噪效果不好。一般h的大小是与噪声标准差σ成正比关系的，通常将其取为噪声的标准方差，即h=σ。

• NL-means的搜索窗口和相似窗口说明

从理论上将，非局部均值滤波的搜索区域大小为整幅图像的大小，即以每个待处理点为中心的图像块需要和整幅图像的所有图像块都要做相似度计算来确定权值，这样计算量是巨大的。与待处理像点为中心的图像块需要进行计算相似度的图像块，我们将其叫做相似窗口。在实际应用中，为了在减少算法复杂度的同时又要尽可能找到比较多的相似度高的图像块，我们会对搜索区域的大小以及相似窗口的大小有固定的参数选择。

在Buades等人的实验中，将搜索区域的窗口大小设为21×21，将相似窗口的大小设为7×7。由于我们的噪声图像大小为N×N，这样将参数固定后的非局部均值滤波算法的时间复杂度是49×441×N^2。事实上，尽管将搜索窗口从整幅图像的大小锁定为21×21的大小，但由于处理图像中每个像素点时，都需要计算以此像素点为中心的图像块与搜索区域中每个图像块之间的欧式距离，用来确定权值。这个过程会存在许多重复的计算，运行此滤波的效率也是相当低的。

---

MATLAB源码实现

函数

function [output]=NLmeansfilter(input,t,f,h)
 %  input: image to be filtered
 %  t: radio of search window 搜索窗口大小，下面实际相似窗口大小为10
 %  f: radio of similarity window 相似窗口半径，下面实际相似窗口大小为3
 %  h: degree of filtering 滤波参数，控制滤波程度

 % Size of the image
 [m n] = size(input);
 
 
 % Memory for the output
 output = zeros(m,n); %输出内存

 % Replicate the boundaries of the input image，复制输入图像的边界
 input2 = padarray(input,[f f],'symmetric');  %padarray扩展图像的函数，扩展m+2*f行，n+2*f列
 
 % Used kernel,高斯核
 kernel = make_kernel(f);
 kernel = kernel / sum(sum(kernel));
 
 h=h*h;
 
 for i=1:m
    for j=1:n        
         i1 = i+ f;
         j1 = j+ f;
    
         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);
         
       for r=rmin:1:rmax
           for s=smin:1:smax
                                               
                if(r==i1 && s==j1)%第一次循环是r=3,s=3;满足if语句，continue；=>第二次循环r=3.s=4=>不满足，到下面W2
                    continue
                end
                                
                W2= input2(r-f:r+f , s-f:s+f);                
                 
                d = sum(sum(kernel.*(W1-W2).*(W1-W2)));%就是N(i)与N(j)的高斯加权距离
                                               
                w=exp(-d/h);%根据加权距离衡量相似度，确定W2窗口与W1窗口的权值w                 
                                 
                if w>wmax                
                    wmax=w;                   
                end
                
                sweight = sweight + w;
                average = average + w*input2(r,s); % 将W2的中心点像素值乘以权值w                               
         end 
         end
             
        average = average + wmax*input2(i1,j1); %最大权值wmax实际就是待处理像素点(i1,j1)所在的图像块，将其累计到average中
        sweight = sweight + wmax;
                   
        if sweight > 0
            output(i,j) = average / sweight;
        else
            output(i,j) = input(i,j);
        end                
 end
 end
 
%高斯核函数
function [kernel] = make_kernel(f)              
kernel=zeros(2*f+1,2*f+1);   
for d=1:f    
  value= 1 / (2*d+1)^2 ;    
  for i=-d:d
    for j=-d:d
    kernel(f+1-i,f+1-j)= kernel(f+1-i,f+1-j) + value ;
    end
  end
end
kernel = kernel ./ f;

        

调用

% demo
clear;
x0=double(imread('Cameraman256.png'));
% add noise
randn('state',0)
sigma=25; % variance of noise 噪声方差
x=x0+5*randn(size(x0));
figure(1), imshow(x,[]),colorbar
disp(['psnr of noisy image=' num2str(psnr(x,x0)) 'dB'])

%--------parameters-------------
% delta for convergence control
f=2;   % radius of square patch in kernelized ar regresssion
t=10;   %searching range in each direction

tic  %秒表计时器
[output]=NLmeansfilter(x,t,f,sigma);%调用函数
toc
figure(2), imshow(output,[])
disp(['psnr of denoising output=' num2str(psnr(output,x0)) 'dB'])

运行效果图

添加了均值为0、方差为25的噪声图

采用相似窗口t=2、搜索窗口f=10、h=25大小非局部均值滤波后效果图

本人对代码关键点的解释（以本代码参数为例）

• f=2，t=2，input：256×256大小； input2:260×260（将边缘对称折叠上去，扩展原input再进行相似度计算）；output：256×256

• 便于理解代码可设置断点，观察参数变化，现将第一次执行参数的值以及示意图如下展示：

  

• 值得注意的一点是：当i，j确定时循环会获得wmax，也就意味着将最大权值赋予待处理像素点所在的图像块，累计到average中得到最终待处理像素点的像素值。

参考文献

[1] Buades A, Coll B, Morel J, et al. A non-local algorithm for image denoising[C]. computer vision and pattern recognition,2005,2(2):60-65.