mobileliker
Opencv2.4.9源码分析——adaptiveBilateralFilter（转）

From：http://blog.csdn.net/zhaocj/article/details/39554747

 

上一篇文章我们介绍了双边滤波，它的公式为：

（1）

其中，，，，f(ξ)表示原图。

 

c(ξ,x)表示的是高斯距离的权值，σd值大则滤波结果会受到更远的像素影响；s(ξ,x)表示的是高斯相似度的权值，σr值大则意味着更无关的像素强度值（或颜色值）会影响滤波器结果。因此这两个值的选取会直接影响到滤波效果。

关于高斯距离的权值，还会受到滤波内核大小的影响，因此它的方差σd值对滤波结果的影响会受到一定的约束，但σr值的选取就难以把握，因此本算法的目的就是自适应的选取σr值的大小。

在opencv文档中没有说明该算法的出处，但从它的程序源码中可以分析得到，σr值是通过领域内的像素值得到，具体公式为：

         （2）

 

其中，n表示邻域内的像素个数，该邻域指的是滤波内核，I(i)表示的是像素值。

 

下面我们来分析一下具体的代码，该函数的原型为：

void adaptiveBilateralFilter(InputArraysrc, OutputArray dst, Size ksize, double sigmaSpace, double maxSigmaColor=20.0,Point anchor=Point(-1, -1), int borderType=BORDER_DEFAULT )

_src为输入原图像；_dst为滤波后的图像；ksize为滤波内核的大小；sigmaSpace为距离权值公式中的方差，即公式1中的σd；maxSigmaColor为相似度权值公式中的方差（σr）的最大值，自适应双边滤波的相似度方差是通过公式2计算得到，但如果计算的结果太大，超过了该值，则以该值为准；anchor为内核锚点；borderType表示用什么方式来处理加宽后的图像四周边界。

 

该函数的源码是在/sources/modules/imgproc/scr/smooth.cpp内：

 

[cpp] view plaincopy
 
  1. void cv::adaptiveBilateralFilter( InputArray _src, OutputArray _dst, Size ksize,  
  2.                                   double sigmaSpace, double maxSigmaColor, Point anchor, int borderType )  
  3. {  
  4.     //得到输入图像矩阵和与其尺寸类型一致的输出图像矩阵  
  5.     Mat src = _src.getMat();  
  6.     _dst.create(src.size(), src.type());  
  7.     Mat dst = _dst.getMat();  
  8.     //该算法只能处理8位二进制的灰度图像和三通道的彩色图像  
  9.     CV_Assert(src.type() == CV_8UC1 || src.type() == CV_8UC3);  
  10.     //得到滤波内核的锚点  
  11.     anchor = normalizeAnchor(anchor,ksize);  
  12.     if( src.depth() == CV_8U )  
  13.         adaptiveBilateralFilter_8u( src, dst, ksize, sigmaSpace, maxSigmaColor, anchor, borderType );  
  14.     else  
  15.         CV_Error( CV_StsUnsupportedFormat,  
  16.         "Adaptive Bilateral filtering is only implemented for 8u images" );  
  17. }  

 

[cpp] view plaincopy
 
  1. static void adaptiveBilateralFilter_8u( const Mat& src, Mat& dst, Size ksize, double sigmaSpace, double maxSigmaColor, Point anchor, int borderType )  
  2. {  
  3.     Size size = src.size();  
  4.     ////处理之前再次检查图像中的相关信息是否正确  
  5.     CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) &&  
  6.               src.type() == dst.type() && src.size() == dst.size() &&  
  7.               src.data != dst.data );  
  8.     //为了在图像边界处得到更好的处理效果，需要对图像四周边界做适当的处理  
  9.     //把原图的四周都加宽，而加宽部分的像素值由borderType值决定  
  10.     //待处理的图像由src换成了temp  
  11.     Mat temp;  
  12.     copyMakeBorder(src, temp, anchor.x, anchor.y, anchor.x, anchor.y, borderType);  
  13.     //通过实例化adaptiveBilateralFilter_8u_Invoker类计算得到自适应双边滤波的结果  
  14.     adaptiveBilateralFilter_8u_Invoker body(dst, temp, ksize, sigmaSpace, maxSigmaColor, anchor);  
  15.     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));  
  16. }  


我们先看一下adaptiveBilateralFilter_8u_Invoker类的构造函数：

 

 

[cpp] view plaincopy
 
  1. adaptiveBilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, Size _ksize, double _sigma_space, double _maxSigmaColor, Point _anchor) :  
  2.     temp(&_temp), dest(&_dest), ksize(_ksize), sigma_space(_sigma_space), maxSigma_Color(_maxSigmaColor), anchor(_anchor)  
  3. {  
  4.     if( sigma_space <= 0 )       //确保σd值不能小于零  
  5.         sigma_space = 1;  
  6.     CV_Assert((ksize.width & 1) && (ksize.height & 1)); //确保滤波内核的宽和高是奇数  
  7.     space_weight.resize(ksize.width * ksize.height);  
  8.     double sigma2 = sigma_space * sigma_space;  
  9.     int idx = 0;  
  10.     int w = ksize.width / 2;  
  11.     int h = ksize.height / 2;  
  12.     //遍历整个内核，计算高斯距离权值  
  13.     for(int y=-h; y<=h; y++)  
  14.         for(int x=-w; x<=w; x++)  
  15.     {  
  16. //在程序中定义了宏ABF_GAUSSIAN，因此高斯距离权值使用的是本文中给出的公式  
  17. #if ABF_GAUSSIAN  
  18.         space_weight[idx++] = (float)exp ( -0.5*(x * x + y * y)/sigma2);  
  19. #else  
  20.         space_weight[idx++] = (float)(sigma2 / (sigma2 + x * x + y * y));  
  21. #endif  
  22.     }  
  23. }  

下面再来介绍adaptiveBilateralFilter_8u_Invoker类中的operator()：

[cpp] view plaincopy
 
  1. virtual void operator()(const Range& range) const  
  2. {  
  3.     int cn = dest->channels();   //得到图像的通道数，即是灰度图像还是彩色图像  
  4.     int anX = anchor.x;  
  5.   
  6.     const uchar *tptr;  
  7.   
  8.     for(int i = range.start;i < range.end; i++)  
  9.     {  
  10.         int startY = i;  
  11.         if(cn == 1) //灰度图像处理方法  
  12.         {  
  13.             float var;      //方差  
  14.             int currVal;        //当前像素值  
  15.             int sumVal = 0;     //方差和  
  16.             int sumValSqr = 0;      //方差的平方和  
  17.             int currValCenter;      //当前内核中心的像素值，即待处理的像素值  
  18.             int currWRTCenter;      //像素的权值  
  19.             float weight;  
  20.             float totalWeight = 0.;  
  21.             float tmpSum = 0.;  
  22.   
  23.             for(int j = 0;j < dest->cols *cn; j+=cn)  
  24.             {  
  25.                 sumVal = 0;  
  26.                 sumValSqr= 0;  
  27.                 totalWeight = 0.;  
  28.                 tmpSum = 0.;  
  29.   
  30.                 // Top row: don't sum the very last element  
  31.                 int startLMJ = 0;  
  32.                 int endLMJ  = ksize.width  - 1;  
  33.                 int howManyAll = (anX *2 +1)*(ksize.width );    //内核像素个数总和  
  34. //在程序前面定义了宏ABF_CALCVAR，因此执行#if的内容  
  35. #if ABF_CALCVAR  
  36.                 //遍历整个内核，计算高斯相似度权值的方差  
  37.                 for(int x = startLMJ; x< endLMJ; x++)  
  38.                 {  
  39.                     //内核中某个像素的指针  
  40.                     tptr = temp->ptr(startY + x) +j;  
  41.                     for(int y=-anX; y<=anX; y++)  
  42.                     {  
  43.                         currVal = tptr[cn*(y+anX)];     //像素值  
  44.                         sumVal += currVal;          //像素之和  
  45.                         sumValSqr += (currVal *currVal);        //像素平方之和  
  46.                     }  
  47.                 }  
  48.                 //由公式2得到方差  
  49.                 var = ( (sumValSqr * howManyAll)- sumVal * sumVal )  /  ( (float)(howManyAll*howManyAll));  
  50.                 //如果计算得到的方差太小，则取值0.01  
  51.                 //如果计算得到的方差超过在调用该函数中所给定的方差，则以给定的方差为准  
  52.                 if(var < 0.01)  
  53.                     var = 0.01f;  
  54.                 else if(var > (float)(maxSigma_Color*maxSigma_Color) )  
  55.                     var =  (float)(maxSigma_Color*maxSigma_Color) ;  
  56.   
  57. #else  
  58.                 var = maxSigmaColor*maxSigmaColor;  
  59. #endif  
  60.                 startLMJ = 0;  
  61.                 endLMJ = ksize.width;  
  62.                 tptr = temp->ptr(startY + (startLMJ+ endLMJ)/2);  
  63.                 currValCenter =tptr[j+cn*anX];      //内核中心像素，即带处理的像素值  
  64.                 //再次遍历内核，这次是由公式1得到输出图像  
  65.                 for(int x = startLMJ; x< endLMJ; x++)  
  66.                 {  
  67.                     tptr = temp->ptr(startY + x) +j;  
  68.                     for(int y=-anX; y<=anX; y++)  
  69.                     {  
  70. #if ABF_FIXED_WEIGHT  
  71.                         weight = 1.0;  
  72. #else  
  73.                         currVal = tptr[cn*(y+anX)];  
  74.                         //内核领域内的像素与内核中心像素之差，  
  75.                         currWRTCenter = currVal - currValCenter;  
  76. //在程序前面定义了宏ABF_GAUSSIAN，因此利用高斯函数得到整个权值（距离权值和相似度权值）  
  77. #if ABF_GAUSSIAN  
  78.                         weight = exp ( -0.5f * currWRTCenter * currWRTCenter/var ) * space_weight[x*ksize.width+y+anX];  
  79. #else  
  80.                         weight = var / ( var + (currWRTCenter * currWRTCenter) ) * space_weight[x*ksize.width+y+anX];  
  81. #endif  
  82.   
  83. #endif  
  84.                         //得到公式1中的分子部分  
  85.                         tmpSum += ((float)tptr[cn*(y+anX)] * weight);  
  86.                         //得到公式1中的分母部分  
  87.                         totalWeight += weight;  
  88.                     }  
  89.                 }  
  90.                 tmpSum /= totalWeight;          //得到公式1的最终结果  
  91.                 //把结果赋值给输出图像  
  92.                 dest->at<uchar>(startY ,j)= static_cast<uchar>(tmpSum);  
  93.             }  
  94.         }  
  95.         else            //  处理彩色图像  
  96.         {  
  97.             assert(cn == 3);  
  98.             float var_b, var_g, var_r;  
  99.             int currVal_b, currVal_g, currVal_r;  
  100.             int sumVal_b= 0, sumVal_g= 0, sumVal_r= 0;  
  101.             int sumValSqr_b= 0, sumValSqr_g= 0, sumValSqr_r= 0;  
  102.             int currValCenter_b= 0, currValCenter_g= 0, currValCenter_r= 0;  
  103.             int currWRTCenter_b, currWRTCenter_g, currWRTCenter_r;  
  104.             float weight_b, weight_g, weight_r;  
  105.             float totalWeight_b= 0., totalWeight_g= 0., totalWeight_r= 0.;  
  106.             float tmpSum_b = 0., tmpSum_g= 0., tmpSum_r = 0.;  
  107.   
  108.             for(int j = 0;j < dest->cols *cn; j+=cn)  
  109.             {  
  110.                 sumVal_b= 0, sumVal_g= 0, sumVal_r= 0;  
  111.                 sumValSqr_b= 0, sumValSqr_g= 0, sumValSqr_r= 0;  
  112.                 totalWeight_b= 0., totalWeight_g= 0., totalWeight_r= 0.;  
  113.                 tmpSum_b = 0., tmpSum_g= 0., tmpSum_r = 0.;  
  114.   
  115.                 // Top row: don't sum the very last element  
  116.                 int startLMJ = 0;  
  117.                 int endLMJ  = ksize.width - 1;  
  118.                 int howManyAll = (anX *2 +1)*(ksize.width);  
  119. #if ABF_CALCVAR  
  120.                 float max_var = (float)( maxSigma_Color*maxSigma_Color);  
  121.                 //遍历内核，分别计算红、绿、蓝三个通道的相似度权值方差  
  122.                 for(int x = startLMJ; x< endLMJ; x++)  
  123.                 {  
  124.                     tptr = temp->ptr(startY + x) +j;  
  125.                     for(int y=-anX; y<=anX; y++)  
  126.                     {  
  127.                         currVal_b = tptr[cn*(y+anX)], currVal_g = tptr[cn*(y+anX)+1], currVal_r =tptr[cn*(y+anX)+2];  
  128.                         sumVal_b += currVal_b;  
  129.                         sumVal_g += currVal_g;  
  130.                         sumVal_r += currVal_r;  
  131.                         sumValSqr_b += (currVal_b *currVal_b);  
  132.                         sumValSqr_g += (currVal_g *currVal_g);  
  133.                         sumValSqr_r += (currVal_r *currVal_r);  
  134.                     }  
  135.                 }  
  136.                 var_b =  ( (sumValSqr_b * howManyAll)- sumVal_b * sumVal_b )  /  ( (float)(howManyAll*howManyAll));  
  137.                 var_g =  ( (sumValSqr_g * howManyAll)- sumVal_g * sumVal_g )  /  ( (float)(howManyAll*howManyAll));  
  138.                 var_r =  ( (sumValSqr_r * howManyAll)- sumVal_r * sumVal_r )  /  ( (float)(howManyAll*howManyAll));  
  139.   
  140.                 if(var_b < 0.01)  
  141.                     var_b = 0.01f;  
  142.                 else if(var_b > max_var )  
  143.                     var_b =  (float)(max_var) ;  
  144.   
  145.                 if(var_g < 0.01)  
  146.                     var_g = 0.01f;  
  147.                 else if(var_g > max_var )  
  148.                     var_g =  (float)(max_var) ;  
  149.   
  150.                 if(var_r < 0.01)  
  151.                     var_r = 0.01f;  
  152.                 else if(var_r > max_var )  
  153.                     var_r =  (float)(max_var) ;  
  154.   
  155. #else  
  156.                 var_b = maxSigma_Color*maxSigma_Color; var_g = maxSigma_Color*maxSigma_Color; var_r = maxSigma_Color*maxSigma_Color;  
  157. #endif  
  158.                 startLMJ = 0;  
  159.                 endLMJ = ksize.width;  
  160.                 tptr = temp->ptr(startY + (startLMJ+ endLMJ)/2) + j;  
  161.                 currValCenter_b =tptr[cn*anX], currValCenter_g =tptr[cn*anX+1], currValCenter_r =tptr[cn*anX+2];  
  162.                 //再次遍历内核，计算最终的结果  
  163.                 for(int x = startLMJ; x< endLMJ; x++)  
  164.                 {  
  165.                     tptr = temp->ptr(startY + x) +j;  
  166.                     for(int y=-anX; y<=anX; y++)  
  167.                     {  
  168. #if ABF_FIXED_WEIGHT  
  169.                         weight_b = 1.0;  
  170.                         weight_g = 1.0;  
  171.                         weight_r = 1.0;  
  172. #else  
  173.                         currVal_b = tptr[cn*(y+anX)];currVal_g=tptr[cn*(y+anX)+1];currVal_r=tptr[cn*(y+anX)+2];  
  174.                         currWRTCenter_b = currVal_b - currValCenter_b;  
  175.                         currWRTCenter_g = currVal_g - currValCenter_g;  
  176.                         currWRTCenter_r = currVal_r - currValCenter_r;  
  177.   
  178.                         float cur_spw = space_weight[x*ksize.width+y+anX];  
  179.   
  180. #if ABF_GAUSSIAN  
  181.                         weight_b = exp( -0.5f * currWRTCenter_b * currWRTCenter_b/ var_b ) * cur_spw;  
  182.                         weight_g = exp( -0.5f * currWRTCenter_g * currWRTCenter_g/ var_g ) * cur_spw;  
  183.                         weight_r = exp( -0.5f * currWRTCenter_r * currWRTCenter_r/ var_r ) * cur_spw;  
  184. #else  
  185.                         weight_b = var_b / ( var_b + (currWRTCenter_b * currWRTCenter_b) ) * cur_spw;  
  186.                         weight_g = var_g / ( var_g + (currWRTCenter_g * currWRTCenter_g) ) * cur_spw;  
  187.                         weight_r = var_r / ( var_r + (currWRTCenter_r * currWRTCenter_r) ) * cur_spw;  
  188. #endif  
  189. #endif  
  190.                         tmpSum_b += ((float)tptr[cn*(y+anX)]   * weight_b);  
  191.                         tmpSum_g += ((float)tptr[cn*(y+anX)+1] * weight_g);  
  192.                         tmpSum_r += ((float)tptr[cn*(y+anX)+2] * weight_r);  
  193.                         totalWeight_b += weight_b, totalWeight_g += weight_g, totalWeight_r += weight_r;  
  194.                     }  
  195.                 }  
  196.                 tmpSum_b /= totalWeight_b;  
  197.                 tmpSum_g /= totalWeight_g;  
  198.                 tmpSum_r /= totalWeight_r;  
  199.   
  200.                 dest->at<uchar>(startY,j  )= static_cast<uchar>(tmpSum_b);  
  201.                 dest->at<uchar>(startY,j+1)= static_cast<uchar>(tmpSum_g);  
  202.                 dest->at<uchar>(startY,j+2)= static_cast<uchar>(tmpSum_r);  
  203.             }  
  204.         }  
  205.     }  
  206. }  

下图是使用adaptiveBilateralFilter(src,dst,Size(7,7),75);的自适应双边滤波的结果。

 



这里还需要说明的是自适应双边滤波adaptiveBilateralFilter要比双边滤波bilareralFilter运行时间更长，而且从源码来看，明显感觉到两个函数不是一个人写的。更重要的是adaptiveBilateralFilter有bug，当滤波内核尺寸取得更大一些的话，比如Size(10,10)，会出现“The application has requested the Runtime toterminate in an unusual way，……”的错误对话框。

posted on 2015-02-02 16:46  mobileliker  阅读(975)  评论(0)  编辑  收藏  举报

Copyright © 2024 mobileliker
Powered by .NET 8.0 on Kubernetes Powered By博客园