bilateral filter

最新推荐文章于 2025-04-04 13:30:00 发布

chezhai

最新推荐文章于 2025-04-04 13:30:00 发布

阅读量652

点赞数

分类专栏： study_C++

study_C++ 专栏收录该内容

50 篇文章

订阅专栏

"Bilateral Filter"一词出自<Bilateral Filtering for Gray and Color Images>这篇文章当中，作者对其的描述是"smooths images while preserving edges"，即在对图像进行平滑滤波的同时又能够保持边缘清晰。这个特点来自于Bilateral Filter的定义，既考虑到了geometric closeness（几何邻近性，即像素点之间的距离），又考虑到了photometric similarity（光学相似性，即像素值之间的距离）。需要注意的是，对于彩色的图像，我们在计算photometric similarity的时候，应该将3个颜色通道同时考虑，否则，单独处理可能会导致颜色失真。

在引出Bilateral Filter这个名词之前，作者还介绍了两种滤波方式：

1.Domain Filtering：传统的滤波方法，在计算权重的时候只考虑周围像素点离中心像素点的距离。

2.Range Filtering：这是作者提出的一种新的滤波方式，在计算权重时只考虑图像的灰度或颜色。（单独作用可能使图像的color map失真）

下面用图简单解释一下这两种滤波方式的差异，假设左边的是原始图像，右边的是滤波后的图像。首先，我们考虑Domain Filtering，在计算滤波后图像的A点的像素值时，原图像中的各点的影响权重取决于它们到原图像中A点的几何距离。即：B和C点的权重是一样的，即使B和C点的像素值相差很大，也不会造成影响。而Range Filtering则不然，就算B和C点到A点的距离一致，哪个点的像素值与A点相差比较小，哪个点的影响权重就比较大。（这种方式相当于对原图像的色彩直方图进行了重新的加权平均，直接用文字表述可能不太直观，在后面我会给出处理结果进行对比。）

将两者结合在一起考虑，就得到了Bilateral Filter。从上面可以看出，Bilateral Filter的定义很简单，只要我们分别定义好两种距离(geometric closeness, photometric similarity)，就可以根据这两种距离得出相应的Bilateral Filter。

下面给出以上三种滤波器的数学表示：

1.Domain Filter

其中，k的作用是归一化，使各权重因子的和为1。x代表中心像素点，f是像素值，c是geometric closeness的度量函数。

2.Range Filter

其中，k的作用同样是归一化。s的photometric similarity的度量函数。

现在对比一下两种滤波方式，可以发现Domain Filter的度量函数刚好作用于图像函数f的domain，而Range Filter的度量函数作用于图像函数f的range。因此，用这两种方式为它们命名。

3.Combine Filter（其实就是Bilateral Filter）

就是简单将两种度量函数乘起来，综合考虑。

下面用一幅图简单讲讲Bilateral Filter，可以看到，输入是一张黑白对分，且有很多起伏噪点的图像。将箭头所指的该点作为中心点，考虑其他像素点对它的影响权重。空间上的权重分布图Gs以及值域上的权重分布图Gr分别计算得出，将两者乘在一起就得到了Bilateral的权重（注意，目前得出的仅仅是针对当前像素点的，其他像素点对其的影响权重）。所有点都用这样的方式计算以后，我们可以看到最终的输出有不错的保边效果。

通常，我们都采用高斯函数作为两种距离的度量函数，定义如下：

其中||.||表示欧式范数，即l2-norm，再简单地说就是“差的平方”。

文中其实还提到了Range Filter与Bilateral Filter的比较，从而得出了两种滤波器融合的必要性。不过那里的一些论证我个人觉得有点不太合理，也可能我还没搞懂，这里就不详细描述了，有兴趣的可以去看看原文。

说了这么多，下面附上我的代码：

实现平台：MATLAB

[plain] view plain copy

print ?

clear;close all;clc;
%Implement of Bilateral Filter (the Gauss Case)
%Also show the domain filter and range filter as comparasion
%% Set paremeter
sigma_d=3; %geometric spread
sigma_r=0.3; %photometric spread
radius=3;
patch=2*radius+1;
%Get the square of distance
dSquare=zeros(patch);
for i=1:patch
for j=1:patch
dSquare(i,j)=(i-radius-1)^2+(j-radius-1)^2;
end
end
%% For gray image
% %Ig=imread('lena.bmp');
% Ig=rgb2gray(imread('lena_small.bmp'));
% figure(1);imshow(Ig);title('Initial image');
% [row,col]=size(Ig);
% %Add noise
% Ig_n=imnoise(Ig,'gaussian',0,0.0005);
% figure(2);imshow(Ig_n);title('Noisy image');
% Ig_n=double(Ig_n)/255;
% Ig_d=zeros(row,col);
% Ig_r=zeros(row,col);
% Ig_b=zeros(row,col);
% for i=1:row
% for j=1:col
% temp=ones(patch);
% img_temp=zeros(patch);
% r=zeros(patch);
% %extract the patch from initial picture
% %deal with the boundary case (a bit complex)
% if(i<=radius)
% temp(1:(radius-i+1),:)=0;
% else if(i>=row-radius+1)
% temp(patch-(i-(row-radius+1)):patch,:)=0;
% end
% end
%
% if(j<=radius)
% temp(:,1:(radius-j+1))=0;
% else if(j>=col-radius+1)
% temp(:,patch-(j-(col-radius+1)):patch)=0;
% end
% end
% %extract
% for m=1:patch
% for n=1:patch
% if(temp(m,n)==1)
% img_temp(m,n)=Ig_n(i-(radius-m+1),j-(radius-n+1));
% end
% end
% end
% %Calculate geometric closeness (c)
% c=exp((-1/2)*(dSquare/sigma_d^2)).*temp;
% Ig_d(i,j)=sum(img_temp.*c)/sum(c);
% %Calculate photemotric similarity (whole) (s_w)
% delta_w=Ig_n-Ig_n(i,j);
% s_w=exp((-1/2)*(delta_w/sigma_r)^2);
% Ig_r(i,j)=sum(double(Ig_n).*s_w)/sum(s_w);
% %Calculate photometric similarity (local) (s)
% delta=img_temp-img_temp(radius+1,radius+1);
% % delta(delta==-img_temp(radius+1,radius+1))=0;
% s=exp((-1/2)*(delta/sigma_r)^2).*c;
% Ig_b(i,j)=sum(img_temp.*s)/sum(s);
% end
% end
%
% %Result of Domain Filter
% figure(3);imshow(uint8(Ig_d*255));title('Domain Filtered Result');
% %Result of Range Filter
% figure(4);imshow(uint8(Ig_r*255));title('Range Filtered Result');
% %Result of Bilateral Filter
% figure(5);imshow(uint8(Ig_b*255));title('Bilateral Filtered Result');
%% For color image
Ig=imread('lena128.bmp');
figure(1);imshow(Ig);title('Initial image');
[row,col,channel]=size(Ig);
%Add noise
Ig_n=imnoise(Ig,'gaussian',0,0.0005);
figure(2);imshow(Ig_n);title('Noisy image');
Ig_n=double(Ig_n)/255;
Ig_d=zeros(row,col,3);
Ig_r=zeros(row,col,3);
Ig_b=zeros(row,col,3);
for i=1:row
for j=1:col
temp=ones(patch);
img_temp=zeros(patch,patch,3);
r=zeros(patch);
%extract the patch from initial picture
%deal with the boundary case (a bit complex)
if(i<=radius)
temp(1:(radius-i+1),:)=0;
else if(i>=row-radius+1)
temp(patch-(i-(row-radius+1)):patch,:)=0;
end
end
if(j<=radius)
temp(:,1:(radius-j+1))=0;
else if(j>=col-radius+1)
temp(:,patch-(j-(col-radius+1)):patch)=0;
end
end
%extract
for m=1:patch
for n=1:patch
if(temp(m,n)==1)
img_temp(m,n,:)=Ig_n(i-(radius-m+1),j-(radius-n+1),:);
end
end
end
%Calculate geometric closeness (c)
c=exp((-1/2)*(dSquare/sigma_d^2)).*temp;
Ig_d(i,j,1)=sum(img_temp(:,:,1).*c)/sum(c);
Ig_d(i,j,2)=sum(img_temp(:,:,2).*c)/sum(c);
Ig_d(i,j,3)=sum(img_temp(:,:,3).*c)/sum(c);
%Calculate photemotric similarity (whole) (s_w)
delta_wr=Ig_n(:,:,1)-Ig_n(i,j,1);
delta_wg=Ig_n(:,:,2)-Ig_n(i,j,2);
delta_wb=Ig_n(:,:,3)-Ig_n(i,j,3);
s_w=exp((-1/2)*((delta_wr^2+delta_wg^2+delta_wb^2)/(sigma_r^2)));
Ig_r(i,j,1)=sum(Ig_n(:,:,1).*s_w)/sum(s_w);
Ig_r(i,j,2)=sum(Ig_n(:,:,2).*s_w)/sum(s_w);
Ig_r(i,j,3)=sum(Ig_n(:,:,3).*s_w)/sum(s_w);
%Calculate photometric similarity (local) (s)
delta_r=img_temp(:,:,1)-img_temp(radius+1,radius+1,1);
delta_g=img_temp(:,:,2)-img_temp(radius+1,radius+1,2);
delta_b=img_temp(:,:,3)-img_temp(radius+1,radius+1,3);
% delta(delta==-img_temp(radius+1,radius+1))=0;
s=exp((-1/2)*((delta_r^2+delta_g^2+delta_b^2)/(sigma_r^2))).*c;
Ig_b(i,j,1)=sum(img_temp(:,:,1).*s)/sum(s);
Ig_b(i,j,2)=sum(img_temp(:,:,2).*s)/sum(s);
Ig_b(i,j,3)=sum(img_temp(:,:,3).*s)/sum(s);
end
end
%Result of Domain Filter
figure(3);imshow(uint8(Ig_d*255));title('Domain Filtered Result');
%Result of Range Filter
figure(4);imshow(uint8(Ig_r*255));title('Range Filtered Result');
%Result of Bilateral Filter
figure(5);imshow(uint8(Ig_b*255));title('Bilateral Filtered Result');

clear;close all;clc;
%Implement of Bilateral Filter (the Gauss Case)
%Also show the domain filter and range filter as comparasion

%% Set paremeter
sigma_d=3;   %geometric spread
sigma_r=0.3;   %photometric spread
radius=3;    
patch=2*radius+1;

%Get the square of distance
dSquare=zeros(patch);
for i=1:patch
    for j=1:patch
        dSquare(i,j)=(i-radius-1)^2+(j-radius-1)^2;
    end
end

%% For gray image
% %Ig=imread('lena.bmp');
% Ig=rgb2gray(imread('lena_small.bmp'));
% figure(1);imshow(Ig);title('Initial image');
% [row,col]=size(Ig);
% %Add noise
% Ig_n=imnoise(Ig,'gaussian',0,0.0005);
% figure(2);imshow(Ig_n);title('Noisy image');
% Ig_n=double(Ig_n)/255;
% Ig_d=zeros(row,col);
% Ig_r=zeros(row,col);
% Ig_b=zeros(row,col);
% for i=1:row
%     for j=1:col
%         temp=ones(patch);
%         img_temp=zeros(patch);
%         r=zeros(patch);
%         %extract the patch from initial picture
%         %deal with the boundary case (a bit complex)
%         if(i<=radius)
%             temp(1:(radius-i+1),:)=0;
%         else if(i>=row-radius+1)
%                 temp(patch-(i-(row-radius+1)):patch,:)=0;
%             end
%         end
%         
%         if(j<=radius)
%             temp(:,1:(radius-j+1))=0;
%         else if(j>=col-radius+1)
%                 temp(:,patch-(j-(col-radius+1)):patch)=0;
%             end
%         end
%         %extract 
%         for m=1:patch
%             for n=1:patch
%                 if(temp(m,n)==1)
%                     img_temp(m,n)=Ig_n(i-(radius-m+1),j-(radius-n+1));
%                 end
%             end
%         end
%         %Calculate geometric closeness (c)
%         c=exp((-1/2)*(dSquare/sigma_d^2)).*temp;
%         Ig_d(i,j)=sum(img_temp.*c)/sum(c);
%         %Calculate photemotric similarity (whole) (s_w)
%         delta_w=Ig_n-Ig_n(i,j);
%         s_w=exp((-1/2)*(delta_w/sigma_r)^2);
%         Ig_r(i,j)=sum(double(Ig_n).*s_w)/sum(s_w);
%         %Calculate photometric similarity (local) (s)
%         delta=img_temp-img_temp(radius+1,radius+1);
% %        delta(delta==-img_temp(radius+1,radius+1))=0;
%         s=exp((-1/2)*(delta/sigma_r)^2).*c;
%         Ig_b(i,j)=sum(img_temp.*s)/sum(s);
%     end
% end
% 
% %Result of Domain Filter
% figure(3);imshow(uint8(Ig_d*255));title('Domain Filtered Result');
% %Result of Range Filter
% figure(4);imshow(uint8(Ig_r*255));title('Range Filtered Result');
% %Result of Bilateral Filter
% figure(5);imshow(uint8(Ig_b*255));title('Bilateral Filtered Result');

%% For color image
Ig=imread('lena128.bmp');
figure(1);imshow(Ig);title('Initial image');
[row,col,channel]=size(Ig);
%Add noise
Ig_n=imnoise(Ig,'gaussian',0,0.0005);
figure(2);imshow(Ig_n);title('Noisy image');
Ig_n=double(Ig_n)/255;
Ig_d=zeros(row,col,3);
Ig_r=zeros(row,col,3);
Ig_b=zeros(row,col,3);
for i=1:row
    for j=1:col
        temp=ones(patch);
        img_temp=zeros(patch,patch,3);
        r=zeros(patch);
        %extract the patch from initial picture
        %deal with the boundary case (a bit complex)
        if(i<=radius)
            temp(1:(radius-i+1),:)=0;
        else if(i>=row-radius+1)
                temp(patch-(i-(row-radius+1)):patch,:)=0;
            end
        end
        
        if(j<=radius)
            temp(:,1:(radius-j+1))=0;
        else if(j>=col-radius+1)
                temp(:,patch-(j-(col-radius+1)):patch)=0;
            end
        end
        %extract 
        for m=1:patch
            for n=1:patch
                if(temp(m,n)==1)
                    img_temp(m,n,:)=Ig_n(i-(radius-m+1),j-(radius-n+1),:);
                end
            end
        end
        %Calculate geometric closeness (c)
        c=exp((-1/2)*(dSquare/sigma_d^2)).*temp;
        Ig_d(i,j,1)=sum(img_temp(:,:,1).*c)/sum(c);
        Ig_d(i,j,2)=sum(img_temp(:,:,2).*c)/sum(c);
        Ig_d(i,j,3)=sum(img_temp(:,:,3).*c)/sum(c);
        %Calculate photemotric similarity (whole) (s_w)
        delta_wr=Ig_n(:,:,1)-Ig_n(i,j,1);
        delta_wg=Ig_n(:,:,2)-Ig_n(i,j,2);
        delta_wb=Ig_n(:,:,3)-Ig_n(i,j,3);
        s_w=exp((-1/2)*((delta_wr^2+delta_wg^2+delta_wb^2)/(sigma_r^2)));
        Ig_r(i,j,1)=sum(Ig_n(:,:,1).*s_w)/sum(s_w);
        Ig_r(i,j,2)=sum(Ig_n(:,:,2).*s_w)/sum(s_w);
        Ig_r(i,j,3)=sum(Ig_n(:,:,3).*s_w)/sum(s_w);
        %Calculate photometric similarity (local) (s)
        delta_r=img_temp(:,:,1)-img_temp(radius+1,radius+1,1);
        delta_g=img_temp(:,:,2)-img_temp(radius+1,radius+1,2);
        delta_b=img_temp(:,:,3)-img_temp(radius+1,radius+1,3);
%        delta(delta==-img_temp(radius+1,radius+1))=0;
        s=exp((-1/2)*((delta_r^2+delta_g^2+delta_b^2)/(sigma_r^2))).*c;
        Ig_b(i,j,1)=sum(img_temp(:,:,1).*s)/sum(s);
        Ig_b(i,j,2)=sum(img_temp(:,:,2).*s)/sum(s);
        Ig_b(i,j,3)=sum(img_temp(:,:,3).*s)/sum(s);
    end
end

%Result of Domain Filter
figure(3);imshow(uint8(Ig_d*255));title('Domain Filtered Result');
%Result of Range Filter
figure(4);imshow(uint8(Ig_r*255));title('Range Filtered Result');
%Result of Bilateral Filter
figure(5);imshow(uint8(Ig_b*255));title('Bilateral Filtered Result');

注：

1.我将灰度图以及彩色图的处理放在了同一个.m文件当中，只需要将不需要的部分注释掉就可以用了。代码写的时候还是写循环写得比较多，还没有进行优化，所以运行速度不会很快。尤其是对于Range Filter，因为不考虑几何距离的影响，我写的时候，每一点都要对全图进行一次遍历计算权重。可能有更快的方式吧，也希望大家能够告诉我。

2.Bilateral Filter对于参数的设置非常敏感，如果对于结果不理想，建议先调调参数。

结果对比：

1.灰度图，图像大小：128*128，sigma_c=3，sigma_s=0.3（后面的参数与此一致）