Cleaned up adaptive bilateral filtering, added support for gaussian…

Cleaned up adaptive bilateral filtering, added support for gaussian interpolation, updated sample and docs

modules/imgproc/doc/filtering.rst

@@ -416,24 +416,23 @@ adaptiveBilateralFilter
 -----------------------
 Applies the adaptive bilateral filter to an image.
 
-.. ocv:function:: void adaptiveBilateralFilter( InputArray src, OutputArray dst, Size ksize, double sigmaSpace, Point anchor=Point(-1, -1), int borderType=BORDER_DEFAULT )
+.. ocv:function:: void adaptiveBilateralFilter( InputArray src, OutputArray dst, Size ksize, double sigmaSpace, double maxSigmaColor = 20.0, Point anchor=Point(-1, -1), int borderType=BORDER_DEFAULT )
 
 .. ocv:pyfunction:: cv2.adaptiveBilateralFilter(src, ksize, sigmaSpace[, dst[, anchor[, borderType]]]) -> dst
 
-    :param src: Source 8-bit, 1-channel or 3-channel image.
-
-    :param dst: Destination image of the same size and type as  ``src`` .
+    :param src: The source image
 
-    :param ksize: filter kernel size.
+    :param dst: The destination image; will have the same size and the same type as src
 
-    :param sigmaSpace: Filter sigma in the coordinate space. It has similar meaning with ``sigmaSpace`` in ``bilateralFilter``.
+    :param ksize: The kernel size. This is the neighborhood where the local variance will be calculated, and where pixels will contribute (in a weighted manner).
 
-    :param anchor: anchor point; default value ``Point(-1,-1)`` means that the anchor is at the kernel center. Only default value is supported now.
+    :param sigmaSpace: Filter sigma in the coordinate space. Larger value of the parameter means that farther pixels will influence each other (as long as their colors are close enough; see sigmaColor). Then d>0, it specifies the neighborhood size regardless of sigmaSpace, otherwise d is proportional to sigmaSpace.
 
-    :param borderType: border mode used to extrapolate pixels outside of the image.
+    :param maxSigmaColor: Maximum allowed sigma color (will clamp the value calculated in the ksize neighborhood. Larger value of the parameter means that more dissimilar pixels will influence each other (as long as their colors are close enough; see sigmaColor). Then d>0, it specifies the neighborhood size regardless of sigmaSpace, otherwise d is proportional to sigmaSpace.
 
-The function applies adaptive bilateral filtering to the input image. This filter is similar to ``bilateralFilter``, in that dissimilarity from and distance to the center pixel is punished. Instead of using ``sigmaColor``, we employ the variance of pixel values in the neighbourhood.
+    :param borderType: Pixel extrapolation method.
 
+A main part of our strategy will be to load each raw pixel once, and reuse it to calculate all pixels in the output (filtered) image that need this pixel value. The math of the filter is that of the usual bilateral filter, except that the sigma color is calculated in the neighborhood, and clamped by the optional input value.
 
 
 blur

modules/imgproc/include/opencv2/imgproc/imgproc.hpp

@@ -400,7 +400,7 @@ CV_EXPORTS_W void bilateralFilter( InputArray src, OutputArray dst, int d,
                                    int borderType=BORDER_DEFAULT );
 //! smooths the image using adaptive bilateral filter
 CV_EXPORTS_W void adaptiveBilateralFilter( InputArray src, OutputArray dst, Size ksize,
-                                           double sigmaSpace, Point anchor=Point(-1, -1),
+                                           double sigmaSpace, double maxSigmaColor = 20.0, Point anchor=Point(-1, -1),
                                            int borderType=BORDER_DEFAULT );
 //! smooths the image using the box filter. Each pixel is processed in O(1) time
 CV_EXPORTS_W void boxFilter( InputArray src, OutputArray dst, int ddepth,

modules/imgproc/src/smooth.cpp

@@ -2279,15 +2279,24 @@ void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,
 
 namespace cv
 {
-#define CALCVAR 1
-#define FIXED_WEIGHT 0
+#ifndef ABF_CALCVAR
+#define ABF_CALCVAR 1
+#endif
+
+#ifndef ABF_FIXED_WEIGHT
+#define ABF_FIXED_WEIGHT 0
+#endif
+
+#ifndef ABF_GAUSSIAN
+#define ABF_GAUSSIAN 1
+#endif
 
 class adaptiveBilateralFilter_8u_Invoker :
     public ParallelLoopBody
 {
 public:
-    adaptiveBilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, Size _ksize, double _sigma_space, Point _anchor) :
-        temp(&_temp), dest(&_dest), ksize(_ksize), sigma_space(_sigma_space), anchor(_anchor)
+    adaptiveBilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, Size _ksize, double _sigma_space, double _maxSigmaColor, Point _anchor) :
+        temp(&_temp), dest(&_dest), ksize(_ksize), sigma_space(_sigma_space), maxSigma_Color(_maxSigmaColor), anchor(_anchor)
     {
         if( sigma_space <= 0 )
             sigma_space = 1;
@@ -2300,7 +2309,11 @@ public:
         for(int y=-h; y<=h; y++)
             for(int x=-w; x<=w; x++)
         {
+#if ABF_GAUSSIAN
+            space_weight[idx++] = (float)exp ( -0.5*(x * x + y * y)/sigma2);
+#else
             space_weight[idx++] = (float)(sigma2 / (sigma2 + x * x + y * y));
+#endif
         }
     }
     virtual void operator()(const Range& range) const
@@ -2336,7 +2349,7 @@ public:
                     int startLMJ = 0;
                     int endLMJ  = ksize.width  - 1;
                     int howManyAll = (anX *2 +1)*(ksize.width );
-#if CALCVAR
+#if ABF_CALCVAR
                     for(int x = startLMJ; x< endLMJ; x++)
                     {
                         tptr = temp->ptr(startY + x) +j;
@@ -2348,8 +2361,14 @@ public:
                         }
                     }
                     var = ( (sumValSqr * howManyAll)- sumVal * sumVal )  /  ( (float)(howManyAll*howManyAll));
+
+                    if(var < 0.01)
+                        var = 0.01f;
+                    else if(var > (float)(maxSigma_Color*maxSigma_Color) )
+                        var =  (float)(maxSigma_Color*maxSigma_Color) ;
+
 #else
-                    var = 900.0;
+                    var = maxSigmaColor*maxSigmaColor;
 #endif
                     startLMJ = 0;
                     endLMJ = ksize.width;
@@ -2360,13 +2379,18 @@ public:
                         tptr = temp->ptr(startY + x) +j;
                         for(int y=-anX; y<=anX; y++)
                         {
-#if FIXED_WEIGHT
+#if ABF_FIXED_WEIGHT
                             weight = 1.0;
 #else
                             currVal = tptr[cn*(y+anX)];
                             currWRTCenter = currVal - currValCenter;
 
-                            weight = var / ( var + (currWRTCenter * currWRTCenter) ) * space_weight[x*ksize.width+y+anX];;
+#if ABF_GAUSSIAN
+                            weight = exp ( -0.5f * currWRTCenter * currWRTCenter/var ) * space_weight[x*ksize.width+y+anX];
+#else
+                            weight = var / ( var + (currWRTCenter * currWRTCenter) ) * space_weight[x*ksize.width+y+anX];
+#endif
+
 #endif
                             tmpSum += ((float)tptr[cn*(y+anX)] * weight);
                             totalWeight += weight;
@@ -2401,7 +2425,8 @@ public:
                     int startLMJ = 0;
                     int endLMJ  = ksize.width - 1;
                     int howManyAll = (anX *2 +1)*(ksize.width);
-#if CALCVAR
+#if ABF_CALCVAR
+                    float max_var = (float)( maxSigma_Color*maxSigma_Color);
                     for(int x = startLMJ; x< endLMJ; x++)
                     {
                         tptr = temp->ptr(startY + x) +j;
@@ -2416,11 +2441,27 @@ public:
                             sumValSqr_r += (currVal_r *currVal_r);
                         }
                     }
-                    var_b = ( (sumValSqr_b * howManyAll)- sumVal_b * sumVal_b )  /  ( (float)(howManyAll*howManyAll));
-                    var_g = ( (sumValSqr_g * howManyAll)- sumVal_g * sumVal_g )  /  ( (float)(howManyAll*howManyAll));
-                    var_r = ( (sumValSqr_r * howManyAll)- sumVal_r * sumVal_r )  /  ( (float)(howManyAll*howManyAll));
+                    var_b =  ( (sumValSqr_b * howManyAll)- sumVal_b * sumVal_b )  /  ( (float)(howManyAll*howManyAll));
+                    var_g =  ( (sumValSqr_g * howManyAll)- sumVal_g * sumVal_g )  /  ( (float)(howManyAll*howManyAll));
+                    var_r =  ( (sumValSqr_r * howManyAll)- sumVal_r * sumVal_r )  /  ( (float)(howManyAll*howManyAll));
+
+                    if(var_b < 0.01)
+                        var_b = 0.01f;
+                    else if(var_b > max_var )
+                        var_b =  (float)(max_var) ;
+
+                    if(var_g < 0.01)
+                        var_g = 0.01f;
+                    else if(var_g > max_var )
+                        var_g =  (float)(max_var) ;
+
+                    if(var_r < 0.01)
+                        var_r = 0.01f;
+                    else if(var_r > max_var )
+                        var_r =  (float)(max_var) ;
+
 #else
-                    var_b = 900.0; var_g = 900.0;var_r = 900.0;
+                    var_b = maxSigma_Color*maxSigma_Color; var_g = maxSigma_Color*maxSigma_Color; var_r = maxSigma_Color*maxSigma_Color;
 #endif
                     startLMJ = 0;
                     endLMJ = ksize.width;
@@ -2431,7 +2472,7 @@ public:
                         tptr = temp->ptr(startY + x) +j;
                         for(int y=-anX; y<=anX; y++)
                         {
-#if FIXED_WEIGHT
+#if ABF_FIXED_WEIGHT
                             weight_b = 1.0;
                             weight_g = 1.0;
                             weight_r = 1.0;
@@ -2442,9 +2483,16 @@ public:
                             currWRTCenter_r = currVal_r - currValCenter_r;
 
                             float cur_spw = space_weight[x*ksize.width+y+anX];
+
+#if ABF_GAUSSIAN
+                            weight_b = exp( -0.5f * currWRTCenter_b * currWRTCenter_b/ var_b ) * cur_spw;
+                            weight_g = exp( -0.5f * currWRTCenter_g * currWRTCenter_g/ var_g ) * cur_spw;
+                            weight_r = exp( -0.5f * currWRTCenter_r * currWRTCenter_r/ var_r ) * cur_spw;
+#else
                             weight_b = var_b / ( var_b + (currWRTCenter_b * currWRTCenter_b) ) * cur_spw;
                             weight_g = var_g / ( var_g + (currWRTCenter_g * currWRTCenter_g) ) * cur_spw;
                             weight_r = var_r / ( var_r + (currWRTCenter_r * currWRTCenter_r) ) * cur_spw;
+#endif
 #endif
                             tmpSum_b += ((float)tptr[cn*(y+anX)]   * weight_b);
                             tmpSum_g += ((float)tptr[cn*(y+anX)+1] * weight_g);
@@ -2468,10 +2516,11 @@ private:
     Mat *dest;
     Size ksize;
     double sigma_space;
+    double maxSigma_Color;
     Point anchor;
     vector<float> space_weight;
 };
-static void adaptiveBilateralFilter_8u( const Mat& src, Mat& dst, Size ksize, double sigmaSpace, Point anchor, int borderType )
+static void adaptiveBilateralFilter_8u( const Mat& src, Mat& dst, Size ksize, double sigmaSpace, double maxSigmaColor, Point anchor, int borderType )
 {
     Size size = src.size();
 
@@ -2481,12 +2530,12 @@ static void adaptiveBilateralFilter_8u( const Mat& src, Mat& dst, Size ksize, do
     Mat temp;
     copyMakeBorder(src, temp, anchor.x, anchor.y, anchor.x, anchor.y, borderType);
 
-    adaptiveBilateralFilter_8u_Invoker body(dst, temp, ksize, sigmaSpace, anchor);
+    adaptiveBilateralFilter_8u_Invoker body(dst, temp, ksize, sigmaSpace, maxSigmaColor, anchor);
     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
 }
 }
 void cv::adaptiveBilateralFilter( InputArray _src, OutputArray _dst, Size ksize,
-                                  double sigmaSpace, Point anchor, int borderType )
+                                  double sigmaSpace, double maxSigmaColor, Point anchor, int borderType )
 {
     Mat src = _src.getMat();
     _dst.create(src.size(), src.type());
@@ -2496,7 +2545,7 @@ void cv::adaptiveBilateralFilter( InputArray _src, OutputArray _dst, Size ksize,
 
     anchor = normalizeAnchor(anchor,ksize);
     if( src.depth() == CV_8U )
-        adaptiveBilateralFilter_8u( src, dst, ksize, sigmaSpace, anchor, borderType );
+        adaptiveBilateralFilter_8u( src, dst, ksize, sigmaSpace, maxSigmaColor, anchor, borderType );
     else
         CV_Error( CV_StsUnsupportedFormat,
         "Adaptive Bilateral filtering is only implemented for 8u images" );

modules/ocl/doc/image_filtering.rst

@@ -497,23 +497,21 @@ ocl::adaptiveBilateralFilter
 --------------------------------
 Returns void
 
-.. ocv:function:: void ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize, double sigmaSpace, Point anchor = Point(-1, -1), int borderType=BORDER_DEFAULT)
+.. ocv:function:: void ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize, double sigmaSpace, double maxSigmaColor = 20.0, Point anchor = Point(-1, -1), int borderType=BORDER_DEFAULT)
 
     :param src: The source image
 
     :param dst: The destination image; will have the same size and the same type as src
 
-    :param ksize: The kernel size
+    :param ksize: The kernel size. This is the neighborhood where the local variance will be calculated, and where pixels will contribute (in a weighted manner).
 
     :param sigmaSpace: Filter sigma in the coordinate space. Larger value of the parameter means that farther pixels will influence each other (as long as their colors are close enough; see sigmaColor). Then d>0, it specifies the neighborhood size regardless of sigmaSpace, otherwise d is proportional to sigmaSpace.
 
-    :param borderType: Pixel extrapolation method.
-
-A main part of our strategy will be to load each raw pixel once, and reuse it to calculate all pixels in the output (filtered) image that need this pixel value.
+    :param maxSigmaColor: Maximum allowed sigma color (will clamp the value calculated in the ksize neighborhood. Larger value of the parameter means that more dissimilar pixels will influence each other (as long as their colors are close enough; see sigmaColor). Then d>0, it specifies the neighborhood size regardless of sigmaSpace, otherwise d is proportional to sigmaSpace.
 
-.. math::
+    :param borderType: Pixel extrapolation method.
 
-    \emph{O}_i = \frac{1}{W_i}\sum\limits_{j\in{N(i)}}{\frac{1}{1+\frac{(V_i-V_j)^2}{\sigma_{N{'}(i)}^2}}*\frac{1}{1+\frac{d(i,j)^2}{\sum^2}}}V_j
+A main part of our strategy will be to load each raw pixel once, and reuse it to calculate all pixels in the output (filtered) image that need this pixel value. The math of the filter is that of the usual bilateral filter, except that the sigma color is calculated in the neighborhood, and clamped by the optional input value.
 
 Local memory organization
 

modules/ocl/include/opencv2/ocl/ocl.hpp

@@ -553,11 +553,12 @@ namespace cv
         CV_EXPORTS void bilateralFilter(const oclMat& src, oclMat& dst, int d, double sigmaColor, double sigmaSpace, int borderType=BORDER_DEFAULT);
 
         //! Applies an adaptive bilateral filter to the input image
-        //  This is not truly a bilateral filter. Instead of using user provided fixed parameters,
-        //  the function calculates a constant at each window based on local standard deviation,
-        //  and use this constant to do filtering.
+        //  Unlike the usual bilateral filter that uses fixed value for sigmaColor,
+        //  the adaptive version calculates the local variance in he ksize neighborhood
+        //  and use this as sigmaColor, for the value filtering. However, the local standard deviation is
+        //  clamped to the maxSigmaColor.
         //  supports 8UC1, 8UC3
-        CV_EXPORTS void adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize, double sigmaSpace, Point anchor = Point(-1, -1), int borderType=BORDER_DEFAULT);
+        CV_EXPORTS void adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize, double sigmaSpace, double maxSigmaColor=20.0, Point anchor = Point(-1, -1), int borderType=BORDER_DEFAULT);
 
         //! computes exponent of each matrix element (dst = e**src)
         // supports only CV_32FC1, CV_64FC1 type

modules/ocl/src/filtering.cpp

@@ -20,6 +20,7 @@
 //    Zero Lin, Zero.Lin@amd.com
 //    Zhang Ying, zhangying913@gmail.com
 //    Yao Wang, bitwangyaoyao@gmail.com
+//    Harris Gasparakis, harris.gasparakis@amd.com
 //
 // Redistribution and use in source and binary forms, with or without modification,
 // are permitted provided that the following conditions are met:
@@ -1407,7 +1408,7 @@ void cv::ocl::GaussianBlur(const oclMat &src, oclMat &dst, Size ksize, double si
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Adaptive Bilateral Filter
 
-void cv::ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize, double sigmaSpace, Point anchor, int borderType)
+void cv::ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize, double sigmaSpace, double maxSigmaColor, Point anchor, int borderType)
 {
     CV_Assert((ksize.width & 1) && (ksize.height & 1));  // ksize must be odd
     CV_Assert(src.type() == CV_8UC1 || src.type() == CV_8UC3);  // source must be 8bit RGB image
@@ -1418,10 +1419,24 @@ void cv::ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize
     int idx = 0;
     int w = ksize.width / 2;
     int h = ksize.height / 2;
-    for(int y=-h; y<=h; y++)
-        for(int x=-w; x<=w; x++)
+
+    int ABF_GAUSSIAN_ocl = 1;
+
+    if(ABF_GAUSSIAN_ocl)
     {
-        lut.at<float>(idx++) = sigma2 / (sigma2 + x * x + y * y);
+        for(int y=-h; y<=h; y++)
+            for(int x=-w; x<=w; x++)
+        {
+            lut.at<float>(idx++) = expf( (float)(-0.5 * (x * x + y * y)/sigma2));
+        }
+    }
+    else
+    {
+        for(int y=-h; y<=h; y++)
+            for(int x=-w; x<=w; x++)
+        {
+            lut.at<float>(idx++) = (float) (sigma2 / (sigma2 + x * x + y * y));
+        }
     }
 
     oclMat dlut(lut);
@@ -1429,7 +1444,7 @@ void cv::ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize
     int cn = src.oclchannels();
 
     normalizeAnchor(anchor, ksize);
-    const static String kernelName = "edgeEnhancingFilter";
+    const static String kernelName = "adaptiveBilateralFilter";
 
     dst.create(src.size(), src.type());
 
@@ -1478,9 +1493,10 @@ void cv::ocl::adaptiveBilateralFilter(const oclMat& src, oclMat& dst, Size ksize
 
     //LDATATYPESIZE is sizeof local data store. This is to exemplify effect of LDS on kernel performance
     sprintf(build_options,
-        "-D VAR_PER_CHANNEL=1 -D CALCVAR=1 -D FIXED_WEIGHT=0 -D EXTRA=%d"
+        "-D VAR_PER_CHANNEL=1 -D CALCVAR=1 -D FIXED_WEIGHT=0 -D EXTRA=%d -D MAX_VAR_VAL=%f -D ABF_GAUSSIAN=%d"
         " -D THREADS=%d -D anX=%d -D anY=%d -D ksX=%d -D ksY=%d -D %s",
-        static_cast<int>(EXTRA), static_cast<int>(blockSizeX), anchor.x, anchor.y, ksize.width, ksize.height, btype);
+        static_cast<int>(EXTRA), static_cast<float>(maxSigmaColor*maxSigmaColor), static_cast<int>(ABF_GAUSSIAN_ocl),
+        static_cast<int>(blockSizeX), anchor.x, anchor.y, ksize.width, ksize.height, btype);
 
     std::vector<pair<size_t , const void *> > args;
     args.push_back(std::make_pair(sizeof(cl_mem), &src.data));

modules/ocl/src/opencl/filtering_adaptive_bilateral.cl

@@ -85,7 +85,7 @@
 #endif
 
 __kernel void
-edgeEnhancingFilter_C4_D0(
+adaptiveBilateralFilter_C4_D0(
     __global const uchar4 * restrict src,
     __global uchar4 *dst,
     float alpha,
@@ -173,14 +173,14 @@ edgeEnhancingFilter_C4_D0(
         //find variance of all data
         int startLMj;
         int endLMj ;
-#if CALCVAR
         // Top row: don't sum the very last element
         for(int extraCnt = 0; extraCnt <=EXTRA; extraCnt++)
         {
+#if CALCVAR
             startLMj = extraCnt;
             endLMj =  ksY+extraCnt-1;
-            sumVal =0;
-            sumValSqr=0;
+            sumVal = (int4)0;
+            sumValSqr= (int4)0;
             for(int j = startLMj; j < endLMj; j++)
                 for(int i=-anX; i<=anX; i++)
                 {
@@ -190,9 +190,10 @@ edgeEnhancingFilter_C4_D0(
                     sumValSqr += mul24(currVal, currVal);
                 }
 
-            var[extraCnt] = convert_float4( ( (sumValSqr * howManyAll)- mul24(sumVal , sumVal) ) ) /  ( (float)(howManyAll*howManyAll) ) ;
+            var[extraCnt] = clamp( convert_float4( ( (sumValSqr * howManyAll)- mul24(sumVal , sumVal) ) ) /  ( (float)(howManyAll*howManyAll) ), (float4)(0.1f, 0.1f, 0.1f, 0.1f), (float4)(MAX_VAR_VAL, MAX_VAR_VAL, MAX_VAR_VAL, MAX_VAR_VAL)) ;
+
 #else
-        var[extraCnt] = (float4)(900.0, 900.0, 900.0, 0.0);
+            var[extraCnt] = (float4)(MAX_VAR_VAL, MAX_VAR_VAL, MAX_VAR_VAL, MAX_VAR_VAL);
 #endif
         }
 
@@ -221,32 +222,48 @@ edgeEnhancingFilter_C4_D0(
 #else
                     weight = 1.0f;
 #endif
-#else
+#else // !FIXED_WEIGHT
                     currVal = convert_int4(data[j][col+anX+i]);
                     currWRTCenter = currVal-currValCenter;
 
+#if ABF_GAUSSIAN
+
 #if VAR_PER_CHANNEL
-                    weight = var[extraCnt] / (var[extraCnt] + convert_float4(currWRTCenter * currWRTCenter)) *
+                    weight = exp( (float4)(-0.5f, -0.5f, -0.5f, -0.5f) * convert_float4(currWRTCenter * currWRTCenter) / var[extraCnt] )*
                         (float4)(lut[lut_j*lut_step+anX+i]);
 #else
-                    weight = 1.0f/(1.0f+( mul24(currWRTCenter.x, currWRTCenter.x) + mul24(currWRTCenter.y, currWRTCenter.y) +
-                        mul24(currWRTCenter.z, currWRTCenter.z))/(var.x+var.y+var.z));
+                    weight = exp( -0.5f * (mul24(currWRTCenter.x, currWRTCenter.x) + mul24(currWRTCenter.y, currWRTCenter.y) +
+                        mul24(currWRTCenter.z, currWRTCenter.z) ) / (var[extraCnt].x+var[extraCnt].y+var[extraCnt].z) ) * lut[lut_j*lut_step+anX+i];
 #endif
+
+#else // !ABF_GAUSSIAN
+
+#if VAR_PER_CHANNEL
+                    weight = var[extraCnt] / (var[extraCnt] + convert_float4(currWRTCenter * currWRTCenter)) *
+                        (float4)(lut[lut_j*lut_step+anX+i]);
+#else
+                    weight = ((float)lut[lut_j*lut_step+anX+i]) /(1.0f+( mul24(currWRTCenter.x, currWRTCenter.x) + mul24(currWRTCenter.y, currWRTCenter.y) +
+                        mul24(currWRTCenter.z, currWRTCenter.z))/(var[extraCnt].x+var[extraCnt].y+var[extraCnt].z));
 #endif
+
+#endif //ABF_GAUSSIAN
+
+
+
+#endif  // FIXED_WEIGHT
+
                     tmp_sum[extraCnt] += convert_float4(data[j][col+anX+i]) * weight;
                     totalWeight += weight;
                 }
             }
 
-            tmp_sum[extraCnt] /= totalWeight;
-
             if(posX >= 0 && posX < dst_cols && (posY+extraCnt) >= 0 && (posY+extraCnt) < dst_rows)
-                dst[(dst_startY+extraCnt) * (dst_step>>2)+ dst_startX + col] = convert_uchar4(tmp_sum[extraCnt]);
+                dst[(dst_startY+extraCnt) * (dst_step>>2)+ dst_startX + col] = convert_uchar4_rtz( (tmp_sum[extraCnt] / (float4)totalWeight) + (float4)0.5f);
 
 #if VAR_PER_CHANNEL
             totalWeight = (float4)(0,0,0,0);
 #else
-            totalWeight = 0;
+            totalWeight = 0.0f;
 #endif
         }
     }
@@ -254,7 +271,7 @@ edgeEnhancingFilter_C4_D0(
 
 
 __kernel void
-edgeEnhancingFilter_C1_D0(
+adaptiveBilateralFilter_C1_D0(
     __global const uchar * restrict src,
     __global uchar *dst,
     float alpha,
@@ -343,10 +360,11 @@ edgeEnhancingFilter_C1_D0(
         //find variance of all data
         int startLMj;
         int endLMj;
-#if CALCVAR
+
         // Top row: don't sum the very last element
         for(int extraCnt=0; extraCnt<=EXTRA; extraCnt++)
         {
+#if CALCVAR
             startLMj = extraCnt;
             endLMj =  ksY+extraCnt-1;
             sumVal = 0;
@@ -361,9 +379,9 @@ edgeEnhancingFilter_C1_D0(
                     sumValSqr += mul24(currVal, currVal);
                 }
             }
-            var[extraCnt] = (float)( ( (sumValSqr * howManyAll)- mul24(sumVal , sumVal) ) ) /  ( (float)(howManyAll*howManyAll) ) ;
+            var[extraCnt] =  clamp((float)( ( (sumValSqr * howManyAll)- mul24(sumVal , sumVal) ) ) /  ( (float)(howManyAll*howManyAll) ) , 0.1f, (float)(MAX_VAR_VAL) );
 #else
-        var[extraCnt] = (float)(900.0);
+            var[extraCnt] = (float)(MAX_VAR_VAL);
 #endif
         }
 
@@ -389,19 +407,20 @@ edgeEnhancingFilter_C1_D0(
                     currVal	= (int)(data[j][col+anX+i])	;
                     currWRTCenter = currVal-currValCenter;
 
+#if ABF_GAUSSIAN
+                    weight = exp( -0.5f * (float)mul24(currWRTCenter,currWRTCenter)/var[extraCnt]) * lut[lut_j*lut_step+anX+i] ;
+#else
                     weight = var[extraCnt] / (var[extraCnt] + (float)mul24(currWRTCenter,currWRTCenter)) * lut[lut_j*lut_step+anX+i] ;
+#endif
 #endif
                     tmp_sum[extraCnt] += (float)(data[j][col+anX+i] * weight);
                     totalWeight += weight;
                 }
             }
 
-            tmp_sum[extraCnt] /= totalWeight;
-
-
             if(posX >= 0 && posX < dst_cols && (posY+extraCnt) >= 0 && (posY+extraCnt) < dst_rows)
             {
-                dst[(dst_startY+extraCnt) * (dst_step)+ dst_startX + col] = (uchar)(tmp_sum[extraCnt]);
+                dst[(dst_startY+extraCnt) * (dst_step)+ dst_startX + col] = convert_uchar_rtz(tmp_sum[extraCnt]/totalWeight+0.5f);
             }
 
             totalWeight = 0;

modules/ocl/test/test_filters.cpp

@@ -338,8 +338,8 @@ OCL_TEST_P(AdaptiveBilateral, Mat)
     {
         random_roi();
 
-        adaptiveBilateralFilter(src_roi, dst_roi, kernelSize, 5, Point(-1, -1), borderType); // TODO anchor
-        ocl::adaptiveBilateralFilter(gsrc_roi, gdst_roi, kernelSize, 5, Point(-1, -1), borderType);
+        adaptiveBilateralFilter(src_roi, dst_roi, kernelSize, 5, 1, Point(-1, -1), borderType); // TODO anchor
+        ocl::adaptiveBilateralFilter(gsrc_roi, gdst_roi, kernelSize, 5, 1, Point(-1, -1), borderType);
 
         Near();
     }

samples/ocl/adaptive_bilateral_filter.cpp

@@ -12,7 +12,9 @@ int main( int argc, const char** argv )
 {
     const char* keys =
         "{ i | input   |          | specify input image }"
-        "{ k | ksize   |     5    | specify kernel size }"
+        "{ k | ksize   |     11   | specify kernel size }"
+        "{ s | sSpace  |     3    | specify sigma space }"
+        "{ c | sColor  |     30   | specify max color }"
         "{ h | help    | false    | print help message }";
 
     CommandLineParser cmd(argc, argv, keys);
@@ -26,27 +28,36 @@ int main( int argc, const char** argv )
 
     string src_path = cmd.get<string>("i");
     int ks = cmd.get<int>("k");
-    const char * winName[] = {"input", "adaptive bilateral CPU", "adaptive bilateral OpenCL", "bilateralFilter OpenCL"};
+    const char * winName[] = {"input", "ABF OpenCL", "BF OpenCL"};
 
-    Mat src = imread(src_path), abFilterCPU;
+    Mat src = imread(src_path);
     if (src.empty())
     {
         cout << "error read image: " << src_path << endl;
         return EXIT_FAILURE;
     }
 
-    ocl::oclMat dsrc(src), dABFilter, dBFilter;
+    double sigmaSpace = cmd.get<int>("s");
+
+    // sigma for checking pixel values. This is used as is in the "normal" bilateral filter,
+    // and it is used as an upper clamp on the adaptive case.
+    double sigmacolor = cmd.get<int>("c");
 
+    ocl::oclMat dsrc(src), dABFilter, dBFilter;
     Size ksize(ks, ks);
-    adaptiveBilateralFilter(src,abFilterCPU, ksize, 10);
-    ocl::adaptiveBilateralFilter(dsrc, dABFilter, ksize, 10);
-    ocl::bilateralFilter(dsrc, dBFilter, ks, 30, 9);
 
+    // ksize is the total width/height of neighborhood used to calculate local variance.
+    // sigmaSpace is not a priori related to ksize/2.
+    ocl::adaptiveBilateralFilter(dsrc, dABFilter, ksize, sigmaSpace, sigmacolor);
+    ocl::bilateralFilter(dsrc, dBFilter, ks, sigmacolor, sigmaSpace);
     Mat abFilter = dABFilter, bFilter = dBFilter;
+
+    ocl::finish();
+
     imshow(winName[0], src);
-    imshow(winName[1], abFilterCPU);
-    imshow(winName[2], abFilter);
-    imshow(winName[3], bFilter);
+    imshow(winName[1], abFilter);
+    imshow(winName[2], bFilter);
+
     waitKey();
 
     return EXIT_SUCCESS;