双边滤波及其matlab代码

参考:
https://blog.csdn.net/MoFMan/article/details/77482794
https://www.jianshu.com/p/8d11e26c9665
代码:https://blog.csdn.net/l_eop/article/details/81812277
高斯滤波:https://blog.csdn.net/nima1994/article/details/79776802

滤波

在介绍双边滤波以前,简要介绍一下滤波的过程.就像下面的动图所展示的,我们使用一个3×3的卷积核以步长为1对图像进行遍历,然后得到一幅新的经过滤波的图像.

其中每次滤波的时候只是在计算中心点的新的像素值.也就是说每次遍历移动一步其实只是在利用邻域(n×n)计算一个像素点的值,直到将图像中所有的点都卷积一次.

那么卷积的结果就是,在3×3的邻域内,权重矩阵和图像块之间各对应位置相乘然后将这些乘积相加得到红色中心点的新的像素值.

双边滤波

双边滤波是一种非线性的方法,同时考虑到了图像的空域信息和灰度相似性,以达到保边去噪的目的.双边滤波通过空域矩阵和值域矩阵形成一个新的权重矩阵,其中空域矩阵用来模糊去噪；值域矩阵用来保护边缘.
在介绍相关公式之前,我们定义(i,j)为中心点坐标,(k,l)为以(i,j)点为中心的邻域S内的任意一点.
空域矩阵中点(k,l)到点(i,j)的空间距离定义为: $$d(i,j,k,l)=exp-\frac{(i-k{)}^2+(j-l{)}^2}{2{\sigma }_d^2}(1)$$
值得注意的一旦${\sigma }_d$确定了,那么空域矩阵中各点的值是始终不变的.
而值域矩阵定义为:$$r(i,j,k,l)=exp-\frac{||f(k,l)-f(i,j)|{|}^2}{2{\sigma }_r^2}(2)$$
两者相乘后，就会产生依赖于数据的双边滤波权重矩阵w(i,j,k,l)：
$$w(i,j,k,l)=exp-\frac{(i-k{)}^2+(j-l{)}^2}{2{\sigma }_d^2}-\frac{||f(k,l)-f(i,j)|{|}^2}{2{\sigma }_r^2}(3)$$
最后计算得到点(i,j)新的像素值g(i,j).

d函数根据像素距离选择权重，距离越近权重越大,这一点和方框滤波,高斯滤波方式相同.而r函数则是根据像素值的差异来分配权值.
在平坦区域,像素差异较小,对应值域权重r(i,j,k,l)接近于1,此时空域权重d(i,j,k,l)起主要作用,相当于直接对此平坦区域进行高斯模糊.需要注意的是空域核是始终不变的,而值域核却因为区域的不同而不同.
边缘的特点就是相距近的点像素值差异大.那么滤波的同时是希望保留边缘的这种特性的.所以在边缘区域,因为像素之间差异较大,即f(k,l)-f(i,j)之间差值大,此时值域核权重变小,导致此处总权重值w(i,j,k,l)下降(w=r*d),当前像素(i,j)受到大差异点的影响就越小,从而保持了边缘的细节信息.

matlab代码

双边滤波的实现代码如下.

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % 双边滤波的程序
% % 
% % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

clear;close all;
% 读取图像
filename = './rgb.png';
% filename = './gray.png';
img = imread(filename);
 
% 设置参数
r = 3;    % 滤波半径
a = 3;    % 全局方差
b = 0.1;  % 局部方差

[~,~,ch] = size(img);
% 判断是灰度图还是彩色图像
if ch == 1   
    g = bfilt_gray(img,r,a,b);
else
    g = bfilt_rgb(img,r,a,b);
end
 
% 显示
figure;subplot(121);imshow(img);
subplot(122);imshow(g);
 
%%  灰度图双边滤波
function res = bfilt_gray(img,r,a,b)
% f灰度图;r滤波半径 ;a全局方差;b局部方差
[x,y] = meshgrid(-r:r);

% 空域核 把中心点当做原点那么各点与中心点的距离就为(i-k)^2 +(j-l)^2 =k^2+l^2.
% 其中(i,j)是中心点坐标,(k,l)是邻域各点坐标
w_spatial = exp(-( x.^2+y.^2 )/(2*a^2));   % 二维高斯函数为 G(x,y) = 1/(2πσ) *  exp(-(x^2 + y^2)/2*σ^2)
img = im2double(img); 
 
[m,n] = size(img);
f_temp = padarray(img,[r r],'symmetric');  % 边缘填充之后的图像
res = zeros(m,n);
count = 0; % 记录有多少个点计算了

for i = r+1:m+r
    for j = r+1:n+r
        
        count = count +1;
        temp = f_temp(i-r:i+r,j-r:j+r); % 一个局部块的像素值
       
        w_pixel = exp(  -( temp- img(i-r,j-r) ).^2/(2*b^2)); % 值域核
        w = w_spatial .* w_pixel;
        s = temp.*w;
    
        res(i-r,j-r) = sum(s(:)) / sum(w(:)); % 计算该点新的像素值
    end

end
fprintf('count = %d \n', count);

end
 
%% 彩色图双边滤波
% 彩色图像每个通道单独处理即可
function res = bfilt_rgb(img,r,a,b)
% f灰度图；r滤波半径；a全局方差；b局部方差

res = zeros(size(img)) ;
for ch = 1:3
     res(:,:,ch) = bfilt_gray(img(:,:,ch),r,a,b);
end

end

实验结果为