双边滤波器（一）

双边滤波其实是来源于高斯滤波，充分的利用了空域和值域的信息，从而得到和好的滤波效果。高斯滤波器

 

       双边滤波（Bilateral filter）是一种非线性的滤波方法，是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，达到保边去噪的目的。具有简单、非迭代、局部的特点。

 

       双边滤波器的好处是可以做边缘保存（edge preserving），一般过去用的维纳滤波或者高斯滤波去降噪，都会较明显地模糊边缘，对于高频细节的保护效果并不明显。双边滤波器顾名思义比高斯滤波多了一个高斯方差sigma－d，它是基于空间分布的高斯滤波函数，所以在边缘附近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保存。但是由于保存了过多的高频信息，对于彩色图像里的高频噪声，双边滤波器不能够干净的滤掉，只能够对于低频信息进行较好的滤波。对于脉冲噪声，双边滤波会把它当成边缘从而不能去除。

 

       在图像处理领域中不管是去雾还是HDR亦或是美容算法的最核心的其实就是这种边缘保持滤波器，从双边滤波的问世以来，人们对它进行了不断的改进，其它的改进算法将在以后的博客中给出。这里先介绍下双边滤波器的基本原理和MATLAB实现。原文链接Bilateral Filtering for Gray and Color Images

 

       C.Tomasi和R.Manduchi提出了一种非迭代的简单策略，用于边缘保持滤波，称之为双边滤波。双边滤波的内在想法是：在图像的值域（range）上做传统滤波器在空域（domain）上做的工作。空域滤波对空间上邻近的点进行加权平均，加权系数随着距离的增加而减少；值域滤波则是对像素值相近的点进行加权平均，加权系数随着值差的增大而减少。

       低通空域滤波器定义如下，f(x)为输入图像，h(x)为输出图像：

       在上式中，度量了邻域中心点x与邻近点的几何邻近度，Kd为归一化参数，其值与图像的内容无关，在相同的几何位置其值恒定。

       同理，值域滤波可定义如下：

      度量了邻域中心点x与邻近点像素的光度相似性，Kr为归一化参数，与Kd不同，Kr的值依赖于图像f。

 

      单单应用值域滤波是没有意义的，它仅仅改变了图像的颜色映射，因为在空间上远离领域中心x的点的取值不应该对x的最终值有所影响。恰当的做法是结合空域与值域滤波，也即同时考虑几何位置与光度大小。组合后的滤波即为双边滤波器：

      这种结合后的滤波器称之为双边滤波器。它用与x点空间邻近且光度相似的点的像素值平均来取代x点上原来的像素值。在平滑的区域，双边滤波器表现为标准的网域滤波器，通过平均过滤掉噪声。如图2（a）所示，存在一条锐利的分界线将图像分为暗和明的区域两个区域，当双边滤波器被定位在明区域上的一个像素点时，相似函数s对同一侧的像素值取近似于1的值，对暗区域的像素点取近似于0的值，如图2（b）所示。滤波操作的结果是，滤波器用中心点周围的明区域像素点的平均值来取代原来中心点的值，而忽略与其邻近的暗区域像素点，反之亦然。因为公式中网域要素的作用，过滤后的图像达到良好的平滑效果，同时因为值域要素的作用，明暗区域的分界线也被很好地保留了（见图2（c））。

 

图2

 

       双边滤波能够在保持强边缘的同时有效地对图像的细小变化进行平滑，所以它被用于噪声去除、细节分解、HDR压缩和图像抽象。Petschnigg和G., Agrawala等人在双边滤波器的基础上提出了联合双边滤波的概念。它与双边滤波不同的是，联合双边滤波器的掩膜权重不是基于输入图像而是基于导向图像进行计算的。这种方法适用于输入图像无法提供准确边缘信息的情况，例如闪光/无闪光去噪、图像上采样、图像去卷积。

 

      虽然双边滤波有上述的优点，但同时也存在着一些不足。其中一点是当它被用于细节分解和HDR压缩时，会引起梯度反转人造物（表现为光晕）的出现。这是因为当一个像素点周围的相似像素很少时，其高斯加权平均是不稳定的。另外，双边滤波的效率偏低，其时间复杂度达到O(Nr2)，当掩膜半径r过大时其计算时间是无法接受的。

 

          关于双边滤波的加速方法会在下次的博文中给出。

 

 

下面是MATLAB的实现代码

CopyRight：

% Douglas R. Lanman, Brown University, September 2006.
% [e-mail], http://mesh.brown.edu/dlanman

 

%简单地说:  %A为给定图像，归一化到[0,1]的矩阵  %W为双边滤波器（核）的边长/2  %定义域方差σd记为SIGMA(1),值域方差σr记为SIGMA(2)    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % Pre-process input and select appropriate filter.  function B = bfilter2(A,w,sigma)    % Verify that the input image exists and is valid.  if ~exist('A','var') || isempty(A)     error('Input image A is undefined or invalid.');  end  if ~isfloat(A) || ~sum([1,3] == size(A,3)) || ...        min(A(:)) < 0 || max(A(:)) > 1     error(['Input image A must be a double precision ',...            'matrix of size NxMx1 or NxMx3 on the closed ',...            'interval [0,1].']);        end    % Verify bilateral filter window size.  if ~exist('w','var') || isempty(w) || ...        numel(w) ~= 1 || w < 1     w = 5;  end  w = ceil(w);    % Verify bilateral filter standard deviations.  if ~exist('sigma','var') || isempty(sigma) || ...        numel(sigma) ~= 2 || sigma(1) <= 0 || sigma(2) <= 0     sigma = [3 0.1];  end    % Apply either grayscale or color bilateral filtering.  if size(A,3) == 1     B = bfltGray(A,w,sigma(1),sigma(2));  else     B = bfltColor(A,w,sigma(1),sigma(2));  end      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % Implements bilateral filtering for grayscale images.  function B = bfltGray(A,w,sigma_d,sigma_r)    % Pre-compute Gaussian distance weights.  [X,Y] = meshgrid(-w:w,-w:w);  %创建核距离矩阵，e.g.  %  [x,y]=meshgrid(-1:1,-1:1)  %   % x =  %   %     -1     0     1  %     -1     0     1  %     -1     0     1  %   %   % y =  %   %     -1    -1    -1  %      0     0     0  %      1     1     1  %计算定义域核  G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));    % Create waitbar.  h = waitbar(0,'Applying bilateral filter...');  set(h,'Name','Bilateral Filter Progress');    % Apply bilateral filter.  %计算值域核H 并与定义域核G 乘积得到双边权重函数F  dim = size(A);  B = zeros(dim);  for i = 1:dim(1)     for j = 1:dim(2)                   % Extract local region.           iMin = max(i-w,1);           iMax = min(i+w,dim(1));           jMin = max(j-w,1);           jMax = min(j+w,dim(2));           %定义当前核所作用的区域为(iMin:iMax,jMin:jMax)           I = A(iMin:iMax,jMin:jMax);%提取该区域的源图像值赋给I                   % Compute Gaussian intensity weights.           H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));                   % Calculate bilateral filter response.           F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);           B(i,j) = sum(F(:).*I(:))/sum(F(:));                      end     waitbar(i/dim(1));  end    % Close waitbar.  close(h);      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % Implements bilateral filter for color images.  function B = bfltColor(A,w,sigma_d,sigma_r)    % Convert input sRGB image to CIELab color space.  if exist('applycform','file')     A = applycform(A,makecform('srgb2lab'));  else     A = colorspace('Lab<-RGB',A);  end    % Pre-compute Gaussian domain weights.  [X,Y] = meshgrid(-w:w,-w:w);  G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));    % Rescale range variance (using maximum luminance).  sigma_r = 100*sigma_r;    % Create waitbar.  h = waitbar(0,'Applying bilateral filter...');  set(h,'Name','Bilateral Filter Progress');    % Apply bilateral filter.  dim = size(A);  B = zeros(dim);  for i = 1:dim(1)     for j = 1:dim(2)                   % Extract local region.           iMin = max(i-w,1);           iMax = min(i+w,dim(1));           jMin = max(j-w,1);           jMax = min(j+w,dim(2));           I = A(iMin:iMax,jMin:jMax,:);                   % Compute Gaussian range weights.           dL = I(:,:,1)-A(i,j,1);           da = I(:,:,2)-A(i,j,2);           db = I(:,:,3)-A(i,j,3);           H = exp(-(dL.^2+da.^2+db.^2)/(2*sigma_r^2));                   % Calculate bilateral filter response.           F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);           norm_F = sum(F(:));           B(i,j,1) = sum(sum(F.*I(:,:,1)))/norm_F;           B(i,j,2) = sum(sum(F.*I(:,:,2)))/norm_F;           B(i,j,3) = sum(sum(F.*I(:,:,3)))/norm_F;                       end     waitbar(i/dim(1));  end    % Convert filtered image back to sRGB color space.  if exist('applycform','file')     B = applycform(B,makecform('lab2srgb'));  else       B = colorspace('RGB<-Lab',B);  end    % Close waitbar.  close(h);