双边滤波

1．概述

    双边滤波，Bilateral filter。是一种可以保边去噪的滤波器。之所以可以达到此去噪效果，是因为滤波器是由两个函数构成。一个函数是由几何空间距离决定滤波器系数。另一个由像素差值决定滤波器系数。
双边滤波器的好处是可以做边缘保存edge preserving，一般过去用的维纳滤波或者高斯滤波去降噪，都会较明显的模糊边缘，对于高频细节的保护效果并不明显。双边滤波器顾名思义比高斯滤波多了一个高斯方差sigma－d，它是基于空间分布的高斯滤波函数，所以在边缘附近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保存。

 1.1 双边滤波

     滤波算法中，目标点上的像素值通常是由其所在位置上的周围的一个小局部邻居像素的值所决定。在2D高斯滤波中的具体实现就是对周围的一定范围内的像素值分别赋以不同的高斯权重值，并在加权平均后得到当前点的最终结果。而这里的高斯权重因子是利用两个像素之间的空间距离（在图像中为2D）关系来生成。通过高斯分布的曲线可以发现，离目标像素越近的点对最终结果的贡献越大，反之则越小。其公式化的描述一般如下所述：

 

其中的c即为基于空间距离的高斯权重，而 用来对结果进行单位化。

高斯滤波在低通滤波算法中有不错的表现，但是其却有另外一个问题，那就是只考虑了像素间的空间位置上的关系，因此滤波的结果会丢失边缘的信息。这里的边缘主要是指图像中主要的不同颜色区域（比如蓝色的天空，黑色的头发等），而Bilateral就是在Gaussian blur中加入了另外的一个权重分部来解决这一问题。Bilateral滤波中对于边缘的保持通过下述表达式来实现：

 

 

其中的s为基于像素间相似程度的高斯权重， 同样用来对结果进行单位化。对两者进行结合即可以得到基于空间距离、相似程度综合考量的Bilateral滤波：

 

 其中k为归一化系数，其表达式为：

 

 

上式中的单位化分部 综合了两种高斯权重于一起而得到，其中的c与s计算可以详细描述如下：

 

2 调试

2.1 matlab调试结果及代码

 

窗口选择较大时，可以看出相对于高斯模糊处理，图像经双边滤波处理后，边缘保持较好，但是窗口较大，处理速度会大大延长。而窗口较小时，滤波不明显。

 

窗口大小：10*10   sigma = [3 0.1]

 

窗口大小3*3   sigma = [3 0.1]

 

 

窗口大小：3*3   sigma = [10 300] 测试结果

 

 

MATLAB代码

 

I=imread('D:\\Administrator\\My Pictures\\Lenagray.bmp'); 

I=double(I)/255;  

w  = 5;       % bilateral filter half-width 

sigma = [3 0.1]; % bilateral filter standard deviations   

I1=bfltGray(I,w,sigma(1),sigma(2));   

subplot(1,2,1); 

imshow(I); 

subplot(1,2,2); 

imshow(I1) 

 

 

% bilateral filtering for grayscale images. 

function B = bfltGray(A,w,sigma_d,sigma_r) 

[X,Y] = meshgrid(-w:w,-w:w); 

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. 

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);%ÌáÈ¡¸ÃÇøÓòµÄԴͼÏñÖµ¸³¸ø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); 

 

2.2 ccs调试结果及代码

2.2.1 双边滤波

 

CCS代码：

 

#include

#include

 

#define IMAGEWIDTH  256

#define IMAGEHEIGHT 256

#define Uint8       unsigned char

 

 

void ReadImage(char *cFileName);

void bmpDataPart(FILE* fpbmp);

void BilateralFilter(double dSigma1, double dSigma2, int nWin);

 

 

 

unsigned char grey[IMAGEHEIGHT][IMAGEWIDTH];

 

void main()

{

 

  

   ReadImage("D:\\Administrator\\My Pictures\\Lenagray.bmp");

  

 

   BilateralFilter(10,100, 1); 

 

  while (1);

                                              

}                                         

 

 

 

void ReadImage(char *cFileName)

{

 

 

        FILE *fp;

 

        if ( fp=fopen(cFileName,"rb" ) )

        {

          

             bmpDataPart(fp);

             fclose(fp);

        }

}

 

 

 

void bmpDataPart(FILE* fpbmp)

{

  int i, j=0;

 

  unsigned char* pix=NULL;

  fseek(fpbmp, 55L, SEEK_SET);

 

  pix=(unsigned char*)malloc(256);

 

  for(j=0;j

   {

    

     fread(pix, 1, 256, fpbmp);

     for(i=0;i

        {

                           

           grey[IMAGEHEIGHT-1-j][i]   =pix[i];

 

        }

 

    }

 

}

 

 

 

 

void BilateralFilter(double dSigma1, double dSigma2, int nWin) 

{ 

    int nNumX,nNumY,nX,nY;

       int nWindows;                                            //窗口边长(2*nWin+1)   大小(2*nWin+1) * (2*nWin+1)

       int i=0,j=0,m=0;

    double dis[9];                                           // 计算距离中间点的几何距离 

    double difGray[9];                                         // 定义中心点到当前点的灰度差                              

    double dData = 0.0,dTotal = 0.0,dGray,dGray1;                   // 用于进行归一化 

                  

 

     nWindows=2*nWin+1;

        

     for(m=0; m 

       { 

             nNumX = m/nWindows;               //行索引 

             nNumY = m%nWindows;               //列索引 

            

                    

             dis[m] = ((nWindows-1)/2-nNumX)*((nWindows-1)/2-nNumX)+     

                         ((nWindows-1)/2-nNumY)*((nWindows-1)/2-nNumY); 

             dis[m] = exp(-0.5*dis[m]/dSigma1/dSigma1); 

              }

     

    //求灰度值的差 

    for(i=50; i 

    { 

        for(j=50; j 

        { 

               

            dGray = grey[i][j];             //当前点的灰度值  

            dData = 0.0;                      

            dTotal = 0.0;                      //用于进行归一化

             

            for(m=0; m 

             { 

                nNumX = m/nWindows;               //行索引 

                nNumY = m%nWindows;               //列索引 

                nX = i-(nWindows-1)/2+nNumX; 

                nY = j-(nWindows-1)/2+nNumY; 

 

             if( (nY>=0) && (nY=0) && (nX

                        {

                // 求S参数

                    dGray1 = grey[nX][nY];

                    difGray[m] = fabs(dGray-dGray1);

                    difGray[m] = exp(-0.5*difGray[m]*difGray[m]/dSigma2/dSigma2);

                   

                //    求C和S参数综合       加权系数求和

                   

                  

                    dData  += dis[m] *dGray1*  difGray[m]; //             

                    dTotal += difGray[m] * dis[m] ;                    

                  

                }    

               

              }

          

             dData /=dTotal;

              if(dData>255)

               dData=255;

             grey[i][j]= dData;

           

        } 

    } 

 

}