双边滤波器 保边滤波器 Bilateral Filters

一. 背景

1. 高斯滤波器根据邻域像素$(x_i, y_i)$和中心点像素$(x, y)$的空间距离作为权值. 为了方便描述, 假设中心点位于原点处$(0, 0)$, 一维高斯滤波器和二维高斯滤波器分别为:
   $G(x,y)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}}$
   $G(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$
2. 实际使用中, 只计算指数项. 所有值计算出来后, 归一化使得算子和为1. 也可以使用查表法求算子值.
3. 高斯滤波器去除背景噪声的同时会模糊边缘. $\sigma$过大边缘模糊, $\sigma$过小起不到去噪效果.
4. 高斯滤波器只使用了空间距离作为权值, 其他的如灰度距离也可以作为权值. 后者正是被双边滤波器利用来实现保边功能.
5. 双边滤波器使用空间距离算子(高斯)和灰度距离算子点乘的方式来生成滤波器算子. 显然不同像素的灰度距离算子是不同的.
   $B(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma_g^2}}e^{-\frac{d(x,y)^2}{2\sigma_r^2}}$
   $\sigma_g$是高斯算子的标准差, $\sigma_r$是灰度距离算子的标准差.
6. 双边滤波是非线性的. 不需要迭代计算.

二. 2维双边滤波器

Matlab Function:

function [ ImgOut ] = BilateralFilter2D( ImgIn, sigma_g, sigma_r, kernel_half_size )
% BILATERALFILTER2D 灰度图双边滤波
% ImgIn, 输入图像,灰度图
% sigma_g, 高斯算子标准差
% sigma_r, 灰度距离算子标准差
% kernel_half_size, 双边算子半长,算子size为2*kernel_half_size+1
% Detailed explanation goes here
n = kernel_half_size;
kernel_g = fspecial('gaussian',[2*n+1,2*n+1],sigma_g);
ImgIn = padarray(ImgIn,[n,n],'symmetric','both');
ImgOut = zeros(size(ImgIn));
for i = n+1:size(ImgIn,1)-nfor j = n+1:size(ImgIn,2)-nkernel_r = exp(-(ImgIn(i-n:i+n,j-n:j+n)-ImgIn(i,j)).^2/2/sigma_r^2);kernel_r = kernel_r/sum(sum(kernel_r));kernel = kernel_g.*kernel_r;kernel = kernel/sum(sum(kernel));ImgOut(i,j) = sum(sum(kernel.*ImgIn(i-n:i+n,j-n:j+n)));end
end
ImgOut = ImgOut(n+1:end-n,n+1:end-n);
end

测试Script:

clear;clc;close all;
ImgIn = imread('E:\图像处理\5.冈萨雷斯图片库\DIP3E_Original_Images_CH06\Fig0635(middle_row_left_chalk ).tif');
ImgIn = rgb2gray(ImgIn);
ImgOut = BilateralFilter2D(single(ImgIn),5,5,14);
figure;
subplot(121);imshow(ImgIn,[]);title('原图');
subplot(122);imshow(ImgOut,[]);title('双边滤波');

输出:

很细的纹理被滤掉, 比较明显的边缘都保留了下来.

三. 双边滤波的缺点

1. 去噪能力远逊于高斯滤波器. $\sigma$和算子直径要比较大才可以;
2. 每个像素都要计算滤波器算子, 计算速度慢.

四. 使用Look Up Table 方法加速滤波过程

function [ ImgOut ] = BilateralFilter2D_U8_C1R( ImgIn, sigma_g, sigma_r, kernel_half_size )
%BILATERALFILTER2D 灰度图双边滤波
% ImgIn, 输入图像,灰度图
% sigma_g, 高斯算子标准差
% sigma_r, 灰度距离算子标准差
% kernel_half_size, 双边算子半长,算子size为2*kernel_half_size+1
n = kernel_half_size;
kernel_g = fspecial('gaussian',[2*n+1,2*n+1],sigma_g);
[X,Y] = meshgrid(0:255,0:255);
kernel_LUT = exp(-(X-Y).^2/2/sigma_r^2);
ImgIn = padarray(ImgIn,[n,n],'symmetric','both');
ImgOut = zeros(size(ImgIn));
for i = n+1:size(ImgIn,1)-nfor j = n+1:size(ImgIn,2)-n% kernel_r = exp(-(ImgIn(i-n:i+n,j-n:j+n)-ImgIn(i,j)).^2/2/sigma_r^2);kernel_idx = ImgIn(i-n:i+n,j-n:j+n)+(ImgIn(i,j)-1)*256;kernel_r = kernel_LUT(kernel_idx);% kernel_r = kernel_LUT(ImgIn(i-n:i+n,j-n:j+n),ImgIn(i,j)*ones(2*n+1,2*n+1));kernel_r = kernel_r/sum(sum(kernel_r));kernel = kernel_g.*kernel_r;kernel = kernel/sum(sum(kernel));ImgOut(i,j) = sum(sum(kernel.*ImgIn(i-n:i+n,j-n:j+n)));end
end
ImgOut = ImgOut(n+1:end-n,n+1:end-n);
end

五. 其他优化

1. http://people.csail.mit.edu/sparis/bf/
2. http://people.csail.mit.edu/jiawen/#code