Bilateral双边滤波器

Halcom

Bilateral双边滤波器
[im1, PSNR] = bif_filter(im,3,0.2);

1. function [out, psn]=bif_filter(im,sigd,sigr)
   
2. % bilateral filter双边滤波器
   
3. % 函数输入：
   
4. %           im    输入的图像
   
5. %           sigd  空间内核的时域参数
   
6. %           sigr  内核参数强度变化范围
   
7. % 函数输出：
   
8. %          out  滤波图像 = output imagespatial kernel
   
9. 
   
10. w=(2*sigd)+1;
    
11. % sigr=(n*100)^2/(.003*(sigd^2));  % 自适应R值，n为高斯噪声强度,n=0.001
    
12. 
    
13. % 高斯滤波器
    
14. gw=zeros(w,w);       % 高斯权值矩阵初始化
    
15. c=ceil(w/2);         % 向前取整
    
16. c=[c c];             % 中心元素位置
    
17. 
    
18. for i=1:w   
    
19.     for j=1:w
    
20.         q=[i,j]; % 记录相连像素位置标识位
    
21.         gw(i,j)=norm(c-q); % 欧氏距离
    
22.     end
    
23. end
    
24. 
    
25. Gwd=(exp(-(gw.^2)/(2*(sigd^2)))); % 高斯函数
    
26. 
    
27. % Padding 扩展图像的边界，防止滑动窗口边界值溢出
    
28. proci=padarray(im,[sigd sigd],'replicate');
    
29. % A = [1 2; 3 4];
    
30. % B = padarray(A,[3 2],'replicate','post')
    
31. % B =
    
32. %      1     2     2     2
    
33. %      3     4     4     4
    
34. %      3     4     4     4
    
35. %      3     4     4     4
    
36. %      3     4     4     4
    
37. 
    
38. [row clm]=size(proci);    % Size of image
    
39. if ~isa(proci,'double')
    
40.     proci = double(proci)/255;   % 转换为double类型
    
41. end
    
42. 
    
43. K=sigd;
    
44. L=[-K:K];
    
45. c=K+1;   % 中心元素位置
    
46. iter=length(L); % 迭代次数
    
47. ind=1;
    
48. 
    
49. for r=(1+K):(row-K)          % 行   
    
50.     for s=(1+K):(clm-K)      % 列     
    
51.             for i=1:iter     % 窗口大小 行
    
52.                 for j=1:iter % 窗口大小 列                  
    
53.                     win(i,j)=proci((r+L(i)),(s+L(j))); % 获取窗口                  
    
54.                 end
    
55.             end
    
56.             I=win; % 灰度矩阵
    
57.             win=win(c,c)-win; % 相对中心点处的强度差异，中心点为参考灰度值
    
58.             win=sqrt(win.^2); % 保证win中的每一个元素为正
    
59.             Gwi=exp(-(win.^2)/(2*(sigr^2))); % 高斯函数      
    
60.             weights=(Gwi.*Gwd)/sum(sum(Gwi.*Gwd)); % 高斯权值
    
61.             Ii=sum(sum(weights.*I));               % 得到当前双边滤波值  
    
62.             proci(r,s)=Ii;                         % 替换当前灰度值
    
63.             win=[];
    
64.     end
    
65. end
    
66. 
    
67. % 移除边界扩展值
    
68. proci=rpadd(proci,K);  % 移除边界扩展值
    
69. out=im2uint8(proci);   % 类型转换
    
70. 
    
71. %% 滤波重建后，图像峰值信噪比计算
    
72. if ~isa(out,'double')
    
73.     dimg = double(out)/255;   % 转换为double类型
    
74. end
    
75. psn = PSN(im,dimg); % PSNR,峰值信噪比
    
76.         
    
77. end
    
78. 
    
79. function x=rpadd(R,K)
    
80. % 移除边界扩展值
    
81. % 函数输入：
    
82. %         R    输入的图像矩阵
    
83. %         K    窗口大小(2*K + 1)
    
84. % 函数输出：
    
85. %         x    移除边界扩展值后的原图像矩阵
    
86. for i=1:K
    
87.     R(1,:)=[];
    
88.     R(:,1)=[];
    
89.     [ro cl]= size(R);
    
90.     R(ro,:)=[];
    
91.     R(:,cl)=[];;
    
92. end
    
93. x=R;
    
94. end
    
95. 
    
96. function [out]=PSN(orgimg,mimg)
    
97. % 峰值信噪比计算
    
98. 
    
99. orgimg =im2double(orgimg);
    
100. mimg   =im2double(mimg);
     
101. Mse=sum(sum((orgimg-mimg).^2))/(numel(orgimg)); % Mean square Error均方差
     
102. out=10*log10(1/Mse);
     
103. end   

复制代码