Non-Local Means滤波器

Halcom

Non-Local Means滤波器
D = Non_Local_Means(im,3,3,0.15);  % 应用Non-Local Means滤波图像

1. function DeNimg = Non_Local_Means(Nimg,PSH,WSH,Sigma)
   
2. % Non_Local_Means滤波器
   
3. %函数输入:
   
4. %         Nimg:  输入的图像矩阵 + 带噪声的
   
5. %         PSH:  扩展窗尺寸大小
   
6. %         WSH: 窗尺寸大小
   
7. %         Sigma：方差
   
8. %函数输出:
   
9. %         DeNimg: 重构滤波图像
   
10. 
    
11. if ~isa(Nimg,'double')
    
12.     Nimg = double(Nimg)/255;
    
13. end
    
14. 
    
15. % 图像维数
    
16. [Height,Width] = size(Nimg);
    
17. u = zeros(Height,Width); % 初始化去噪图像矩阵
    
18. M = u; % 初始化权值矩阵
    
19. Z = M; % 初始化叠加权值 accumlated weights
    
20. % 避免边界效应
    
21. PP = padarray(Nimg,[PSH,PSH],'symmetric','both');
    
22. PW = padarray(Nimg,[WSH,WSH],'symmetric','both');
    
23. % padarray使用
    
24. % A =
    
25. %      1     3     4
    
26. %      2     3     4
    
27. %      3     4     5
    
28. % B = padarray(A, 2 * [1 1], 0, 'both')
    
29. %      0     0     0     0     0     0     0
    
30. %      0     0     0     0     0     0     0
    
31. %      0     0     1     3     4     0     0
    
32. %      0     0     2     3     4     0     0
    
33. %      0     0     3     4     5     0     0
    
34. %      0     0     0     0     0     0     0
    
35. %      0     0     0     0     0     0     0
    
36. % 主循环
    
37. for dx = -WSH:WSH
    
38.     for dy = -WSH:WSH
    
39.         if dx ~= 0 || dy ~= 0
    
40.             Sd = integral_img(PP,dx,dy);  % 插值图像
    
41.             % 获取对应像素点的平方差矩阵
    
42.             SDist = Sd(PSH+1:end-PSH,PSH+1:end-PSH)+Sd(1:end-2*PSH,1:end-2*PSH)-Sd(1:end-2*PSH,PSH+1:end-PSH)-Sd(PSH+1:end-PSH,1:end-2*PSH);      
    
43.             % 计算每一个像素点的权值
    
44.             w = exp(-SDist/(2*Sigma^2));
    
45.             % 得到相应的噪声点
    
46.             v = PW((WSH+1+dx):(WSH+dx+Height),(WSH+1+dy):(WSH+dy+Width));
    
47.             % 更新去噪图像矩阵
    
48.             u = u+w.*v;
    
49.             % 更新权值去噪图像矩阵
    
50.             M = max(M,w);
    
51.             % 更新叠加权值 accumlated weights
    
52.             Z = Z+w;
    
53.         end
    
54.     end
    
55. end
    
56. % 重构图像        
    
57. f = 1;
    
58. u = u+f*M.*Nimg;
    
59. u = u./(Z+f*M);
    
60. DeNimg = u; % 重构去噪图像
    
61. 
    
62. function Sd = integral_img(v,dx,dy)
    
63. % 根据平方差，插值图像
    
64. % 变换计算：tx = vx+dx; ty = vy+dy
    
65. t = img_Shift(v,dx,dy);
    
66. % 平方差图像
    
67. diff = (v-t).^2;
    
68. % 沿行插值
    
69. Sd = cumsum(diff,1);  % 行叠加
    
70. % 沿列插值
    
71. Sd = cumsum(Sd,2);    % 列叠加
    
72. 
    
73. function t = img_Shift(v,dx,dy)
    
74. % 在xy坐标系下，进行图像变换操作
    
75. t = zeros(size(v));
    
76. type = (dx>0)*2+(dy>0);
    
77. switch type
    
78.     case 0 % dx<0,dy<0: 向右下方移动
    
79.         t(-dx+1:end,-dy+1:end) = v(1:end+dx,1:end+dy);
    
80.     case 1 % dx<0,dy>0: 向左下方移动
    
81.         t(-dx+1:end,1:end-dy) = v(1:end+dx,dy+1:end);
    
82.     case 2 % dx>0,dy<0: 向右上方移动
    
83.         t(1:end-dx,-dy+1:end) = v(dx+1:end,1:end+dy);
    
84.     case 3 % dx>0,dy>0: 向左上方移动
    
85.         t(1:end-dx,1:end-dy) = v(dx+1:end,dy+1:end);
    
86. end

复制代码