图像算法之十二：非局部均值滤波及其Matlab实现

保边去噪算法之二：

首先谈一下什么是非局部均值滤波。在此之前，我们先来看一下均值滤波的原理。

均值滤波

均值滤波的计算非常简单，将图像像素点灰度记录在数组中，然后设置方框半径的值，然后将方框中的所有点的像素求和取平均，得到的结果就是均值滤波后对应像素点的灰度值。 
优点： 
计算很快而且简单 
从算法可以看出，只是求了平均，并没有很复杂的计算 
缺点： 
得到的图像很模糊 
当方框的半径越大，得到的图像中那些变化较大的地方（边缘）计算后变化就越小，即边缘不明显，即模糊

非局部均值滤波

非局部均值滤波的基本原理与均值滤波类似，都是要取平均值，但是非局部均值滤波在计算中加入了每一个点的权重值，所以能够保证在相邻且相差很大的点在方框中求平均值时相互之间的影响减小，也就对图像边缘细节部分保留很多，这样图像看起来会更清晰。非局部均值滤波的算法我认为可以大致分为以下几个步骤： 
1. 首先在一个点A周围取一个大的框（搜索框），设边长为s，A在方框的中心，然后再在方框中取小的方框，即相似框，设边长为d 
2. 那么在A周围也有一个边长为d的方框，然后在大方框中找到所有边长为d的小方框的组合（就是一个小正方形在一个大正方形中到处移动，记录小正方形中心点的坐标就行了），设小方框的中心点为B，分别于A周围的相似框求减法，并且加入高斯核计算得到的加权值，这样可以计算出一个二维数组，里面存放着各个点的差值乘以权重后的值，加入高斯核主要是因为距离中心点距离不同对中心点的影响大小也不同，而且高斯核的权重和是1，所以就不用再归一化了。 
3. 然后将这个二维数组求和并平均，得到的值就是这个相似框的中心点B对于A的权重值。计算出A周围所有点的权重值，其实这个时候这个值和权重是成反比的，以A本身为例（以A为中心点的相似框），计算出来A对于A的所谓权重值是零。然后根据计算出来的值用一个指数减函数就得到了成正比的权重关系，具体的函数见下面的代码，w=exp(-d/h)，就是这个，其中d就是计算出来的值啦，代入后w就是成正比的权重关系啦，h是一个滤波百分比值。可以先固定为一个常数。 而且这个计算出来w就是一个(0,1)的值哦，自动归一化啦 
4. 然后就是根据得到的权重值以及各个点本身的灰度值计算出非局部均值滤波后A点的灰度值。 
5. 以此类推，可以计算出图中所有点经过非局部均值滤波后的值 
优点 
可以既去除噪声，又保留图像边缘细节 
当然去噪声指的一般是高斯白噪声，因为高斯白噪声的均值是0，所以求和取平均会比较有效果 
缺点 
计算起来很慢 
如果图像像素点比较多，而且计算的时候取的框还比较大的话，那么计算一般几分钟是要的了。 
下面是一段非局部均值滤波的代码：

function [output]=NLmeansfilter(input,t,f,h)  [m n]=size(input); %f:搜索框半径；t:相似框半径；h:滤波频率百分比 Output=zeros(m,n); input2 = padarray(input,[f f],‘symmetric‘); %将边缘对称折叠上去  % f :加宽的宽度值 kernel = make_kernel(f); %计算得到一个高斯核，用于后续的计算 kernel = kernel / sum(sum(kernel)); %sum：对矩阵k的每一列求和，k表示矩阵k的总和 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) 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 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; 
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
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

图像算法之十二：非局部均值滤波及其Matlab实现

原文：http://blog.csdn.net/soaringlee_fighting/article/details/54632426