Bilateral Filtering双边滤波

1.Introduction

双边滤波属于局部滤波的一种，本质是利用局部纹理的相似性，对依赖距离的权重滤波器（Gaussian，mean Filtering）等进行再加权，典型的模型如：

           

图像示例即：

                        

一个简单的matlab版本双边滤波demo如下，此处定义距离加权因子为：

                                                                                    

clc;
clear all;
close all;

image=imread('src1.bmp');
image=double(image);

%1.normal gaussian filtering
lpf=fspecial('gaussian',5,3);
dst=imfilter(image,lpf,'replicate');

%2.simple bilateral Filtering demo
[row,col]=size(image);
dst1=image;
fz=5;
for i=3:row-2
    for j=3:col-2
        src_block=image(i-2:i+2,j-2:j+2);
        dw=1./(1+abs(src_block-image(i,j)));
        w=dw.*lpf;
        w=w./sum(abs(w(:)));
        dst1(i,j)=sum(sum(w.*src_block));
    end
end

%display result
figure,subplot(2,2,1),imshow(image,[]);title('src');
hold on,subplot(2,2,2),imshow(dst,[]);title('dst-gaussian');
hold on,subplot(2,2,3),imshow(dst1,[]);title('dst1-bilateral');
hold on,subplot(2,2,4),imshow(dst-dst1,[]);title('diff=dst-dst1');

效果图如下：

通过diff图像，很明显双边滤波能很好的保留图像边缘，而平坦区域得到一定程度的平滑。

2.加速的双边滤波

由于无法进行行列分离且依赖局部信息定义卷积核，理想的双边滤波较为耗时，Kunal Narayan Chaudhury等在论文《Fast O(1) bilateral filtering using trigonometric range kernels》提出，基于升余弦Raised cosines函数来逼近高斯函数，利用一些特性把值域函数分解为一系列函数的叠加，以此，可以分解两步：（1）计算分解图像；（2）对分解进行空域高斯滤波（此处可以利用傅立叶变换在频域求解，对于尺寸为M的图像，单次计算时间复杂度，继而normalization。

具体方法如：

实现代码如：

clc;
clear all;
close all;

image=imread('src1.bmp');
image=double(image);

%1.init
deltaSpatial=0.4;
deltaRange=0.5;
T=max(image(:));
gamma=pi/(2*T);
rho=gamma*deltaRange;
N=3;
f=image;

%2.auxlliary image
for n=0:N
    temp=gamma*(2*n-N)*f;
    b(:,:,n+1)=exp(1i*temp/(rho*sqrt(N)));
    g(:,:,n+1)=b(:,:,n+1).*f;
    c=factorial(N)/(factorial(n)*factorial(N-n));
    d(:,:,n+1)=c*exp(-1i*temp/(rho*sqrt(N)))/(2.^N);
end

%3.gaussian filtering
size=15;
lpf=fspecial('gaussian',size,deltaSpatial);
for n=0:N
    b(:,:,n+1)=imfilter(b(:,:,n+1),lpf,'replicate');
    g(:,:,n+1)=imfilter(g(:,:,n+1),lpf,'replicate');
end

%4.display result
result=abs(sum(d.*g,3)./sum(d.*b,3));

figure,subplot(1,3,1),imshow(f,[]);title('src');
hold on,subplot(1,3,2),imshow(result,[]);title('dst');
hold on,subplot(1,3,3),imshow(f-result,[]);title('diff');