非局部均值滤波matlab代码理解

左手的ㄟ右手 2022-06-15 02:28 471阅读 0赞

以下代码转载自：http://blog.csdn.net/frankgoogle/article/details/52209901

本文在以上博主的分析下进一步理解，主要是对高斯加权距离的理解。非局部均值滤波的代码如下：

function \[output\]=NLmeansfilter(input,t,f,h)  
\[m n\]=size(input); %t:搜索框半径；f:相似框半径；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;

在运行过程中输入代码：

load('IMG\_Normal002.mat')  
imag=Img(:,:,50);  
imshow(imag,\[0 255\])  
fima=NLmeansfilter(imag,10,3,10);

![Center][]

也就是对上图进行非局部均值滤波（这是一张脑部的MRI图像，大小为448×448）

其中，搜索框半径设置为10，相似框半径设置为3，滤波频率设置为10.

那么，程序运行大致过程如下。先扩充原图像为545×454（各边以相似框半径扩充）；然后以像素点p（4，4）为中心，构成一个7×7的相似框W1，设置搜索框范围为(行：4-14，列：4-14），再以点p1（4，5）为中心，构成同样大小的相似框W2。然后就这两个相似框的高斯加权距离d = sum(sum(kernel.\*(W1-W2).\*(W1-W2)));

这里为什么叫高斯加权距离呢？是因为(W1-W2).\*(W1-W2）计算出两个相似框对应像素差值的平方W3，在用高斯核进行卷积，相当于把W3中的各像素点按高斯模板赋予权重（离中心点近的权重大些，远的就权重小些），然后再求和得到高斯加权距离。而普通的欧氏距离就是d = sum(sum((W1-W2).\*(W1-W2)));

同理，循环进行，最后得到结果如下

![Center 1][]

至于h我现在还不知道该怎么理解，如果大神知道，望解惑，谢谢！![微笑][smile.gif]

[Center]: /images/20220615/5511f0ab074d4f6f8090c855bcb6a575.png
[Center 1]: /images/20220615/3ee539a76f01417289aaeaa96980553b.png
[smile.gif]: /images/20220615/820b17df040d470c95c6db872747ecc0.png