parallel version of bilateral filter was implemented using parallel_for_

modules/core/include/opencv2/core/parallel_tool.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_PARALLEL_TOOL_HPP__
+#define __OPENCV_PARALLEL_TOOL_HPP__
+
+#ifdef HAVE_CVCONFIG_H
+# include <cvconfig.h>
+#endif // HAVE_CVCONFIG_H
+
+/*
+    HAVE_TBB - using TBB
+    HAVE_GCD - using GCD
+    HAVE_OPENMP - using OpenMP
+    HAVE_CONCURRENCY - using visual studio 2010 concurrency
+*/
+
+#ifdef HAVE_TBB
+#  include "tbb/tbb_stddef.h"
+#  if TBB_VERSION_MAJOR*100 + TBB_VERSION_MINOR >= 202
+#    include "tbb/tbb.h"
+#    include "tbb/task.h"
+#    undef min
+#    undef max
+#  else
+#    undef HAVE_TBB
+#   endif // end TBB version
+#endif // HAVE_TBB
+
+#ifdef __cplusplus
+
+namespace cv
+{
+    // a base body class
+    class CV_EXPORTS ParallelLoopBody
+    {
+    public:
+        virtual void operator() (const Range& range) const = 0;
+        virtual ~ParallelLoopBody();
+    };
+
+    CV_EXPORTS void parallel_for_(const Range& range, const ParallelLoopBody& body);
+
+    template <typename Iterator, typename Body> inline
+    CV_EXPORTS void parallel_do_(Iterator first, Iterator last, const Body& body)
+    {
+#ifdef HAVE_TBB
+        tbb::parallel_do(first, last, body);
+#else
+        for ( ; first != last; ++first)
+            body(*first);
+#endif // HAVE_TBB
+    }
+
+    template <typename Body> inline
+    CV_EXPORTS void parallel_reduce_(const Range& range, Body& body)
+    {
+#ifdef HAVE_TBB
+        tbb::parallel_reduce(tbb::blocked_range<int>(range.start, range.end), body);
+#else
+        body(range);
+#endif // end HAVE_TBB
+    }
+
+} // namespace cv
+
+#endif // __cplusplus
+
+#endif // __OPENCV_PARALLEL_TOOL_HPP__

modules/core/src/parallel_tool.cpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "precomp.hpp"
+
+#ifdef HAVE_CONCURRENCY
+#  include <ppl.h>
+#elif defined HAVE_OPENMP
+#  include <omp.h>
+#elif defined HAVE_GCD
+#  include <dispatch/dispatch.h>
+#endif // HAVE_CONCURRENCY
+
+namespace cv
+{
+    ParallelLoopBody::~ParallelLoopBody() { }
+
+#ifdef HAVE_TBB
+    class TbbProxyLoopBody
+    {
+    public:
+        TbbProxyLoopBody(const ParallelLoopBody& _body) :
+            body(&_body)
+        { }
+
+        void operator ()(const tbb::blocked_range<int>& range) const
+        {
+            body->operator()(Range(range.begin(), range.end()));
+        }
+
+    private:
+        const ParallelLoopBody* body;
+    };
+#endif // end HAVE_TBB
+
+#ifdef HAVE_GCD
+    static
+    void block_function(void* context, size_t index)
+    {
+        ParallelLoopBody* ptr_body = static_cast<ParallelLoopBody*>(context);
+        ptr_body->operator()(Range(index, index + 1));
+    }
+#endif // HAVE_GCD
+
+    void parallel_for_(const Range& range, const ParallelLoopBody& body)
+    {
+#ifdef HAVE_TBB
+
+        tbb::parallel_for(tbb::blocked_range<int>(range.start, range.end), TbbProxyLoopBody(body));
+
+#elif defined HAVE_CONCURRENCY
+
+        Concurrency::parallel_for(range.start, range.end, body);
+
+#elif defined HAVE_OPENMP
+
+#pragma omp parallel for schedule(dynamic)
+        for (int i = range.start; i < range.end; ++i)
+            body(Range(i, i + 1));
+
+#elif defined (HAVE_GCD)
+
+        dispatch_queue_t concurrent_queue = dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0);
+        dispatch_apply_f(range.end - range.start, concurrent_queue, &const_cast<ParallelLoopBody&>(body), block_function);
+
+#else
+
+        body(range);
+
+#endif // end HAVE_TBB
+    }
+
+} // namespace cv

modules/core/src/precomp.hpp

@@ -50,6 +50,7 @@
 #include "opencv2/core/core.hpp"
 #include "opencv2/core/core_c.h"
 #include "opencv2/core/internal.hpp"
+#include "opencv2/core/parallel_tool.hpp"
 
 #include <assert.h>
 #include <ctype.h>

modules/imgproc/perf/perf_bilateral.cpp

+#include "perf_precomp.hpp"
+
+using namespace std;
+using namespace cv;
+using namespace perf;
+using namespace testing;
+using std::tr1::make_tuple;
+using std::tr1::get;
+
+CV_ENUM(Mat_Type, CV_8UC1, CV_8UC3, CV_32FC1, CV_32FC3)
+
+typedef TestBaseWithParam< tr1::tuple<Size, int, Mat_Type> > TestBilateralFilter;
+
+PERF_TEST_P( TestBilateralFilter, BilateralFilter,
+             Combine(
+                Values( szVGA, sz1080p ), // image size
+                Values( 3, 5 ), // d
+                ValuesIn( Mat_Type::all() ) // image type
+             )
+)
+{
+    Size sz;
+    int d, type;
+    const double sigmaColor = 1., sigmaSpace = 1.;
+
+    sz         = get<0>(GetParam());
+    d          = get<1>(GetParam());
+    type       = get<2>(GetParam());
+
+    Mat src(sz, type);
+    Mat dst(sz, type);
+
+    declare.in(src, WARMUP_RNG).out(dst).time(20);
+
+    TEST_CYCLE() bilateralFilter(src, dst, d, sigmaColor, sigmaSpace, BORDER_DEFAULT);
+
+    SANITY_CHECK(dst);
+}

modules/imgproc/src/precomp.hpp

@@ -50,6 +50,7 @@
 #include "opencv2/imgproc/imgproc.hpp"
 #include "opencv2/imgproc/imgproc_c.h"
 #include "opencv2/core/internal.hpp"
+#include "opencv2/core/parallel_tool.hpp"
 #include <math.h>
 #include <assert.h>
 #include <string.h>

modules/imgproc/src/smooth.cpp

@@ -1288,48 +1288,119 @@ void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize )
 namespace cv
 {
 
+class BilateralFilter_8u_Invoker :
+    public ParallelLoopBody
+{
+public:
+    BilateralFilter_8u_Invoker(const Mat &_src, Mat& _dst, Mat _temp, int _radius, int _maxk,
+        int* _space_ofs, float *_space_weight, float *_color_weight) :
+        ParallelLoopBody(), src(_src), dst(_dst), temp(_temp), radius(_radius),
+        maxk(_maxk), space_ofs(_space_ofs), space_weight(_space_weight), color_weight(_color_weight)
+    {
+    }
+    
+    virtual void operator() (const Range& range) const
+    {
+        int i, j, cn = src.channels(), k;
+        Size size = src.size();
+        
+        for( i = range.start; i < range.end; i++ )
+        {
+            const uchar* sptr = temp.data + (i+radius)*temp.step + radius*cn;
+            uchar* dptr = dst.data + i*dst.step;
+            
+            if( cn == 1 )
+            {
+                for( j = 0; j < size.width; j++ )
+                {
+                    float sum = 0, wsum = 0;
+                    int val0 = sptr[j];
+                    for( k = 0; k < maxk; k++ )
+                    {
+                        int val = sptr[j + space_ofs[k]];
+                        float w = space_weight[k]*color_weight[std::abs(val - val0)];
+                        sum += val*w;
+                        wsum += w;
+                    }
+                    // overflow is not possible here => there is no need to use CV_CAST_8U
+                    dptr[j] = (uchar)cvRound(sum/wsum);
+                }
+            }
+            else
+            {
+                assert( cn == 3 );
+                for( j = 0; j < size.width*3; j += 3 )
+                {
+                    float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
+                    int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
+                    for( k = 0; k < maxk; k++ )
+                    {
+                        const uchar* sptr_k = sptr + j + space_ofs[k];
+                        int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
+                        float w = space_weight[k]*color_weight[std::abs(b - b0) +
+                                                               std::abs(g - g0) + std::abs(r - r0)];
+                        sum_b += b*w; sum_g += g*w; sum_r += r*w;
+                        wsum += w;
+                    }
+                    wsum = 1.f/wsum;
+                    b0 = cvRound(sum_b*wsum);
+                    g0 = cvRound(sum_g*wsum);
+                    r0 = cvRound(sum_r*wsum);
+                    dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
+                }
+            }
+        }
+    }
+    
+private:
+    const Mat& src;
+    Mat &dst, temp;
+    int radius, maxk, * space_ofs;
+    float *space_weight, *color_weight;
+};
+
 static void
 bilateralFilter_8u( const Mat& src, Mat& dst, int d,
-                    double sigma_color, double sigma_space,
-                    int borderType )
+    double sigma_color, double sigma_space,
+    int borderType )
 {
     int cn = src.channels();
-    int i, j, k, maxk, radius;
+    int i, j, maxk, radius;
     Size size = src.size();
-
+    
     CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
-        src.type() == dst.type() && src.size() == dst.size() &&
-        src.data != dst.data );
-
+              src.type() == dst.type() && src.size() == dst.size() &&
+              src.data != dst.data );
+    
     if( sigma_color <= 0 )
         sigma_color = 1;
     if( sigma_space <= 0 )
         sigma_space = 1;
-
+    
     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
-
+    
     if( d <= 0 )
         radius = cvRound(sigma_space*1.5);
     else
         radius = d/2;
     radius = MAX(radius, 1);
     d = radius*2 + 1;
-
+    
     Mat temp;
     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
-
+    
     vector<float> _color_weight(cn*256);
     vector<float> _space_weight(d*d);
     vector<int> _space_ofs(d*d);
     float* color_weight = &_color_weight[0];
     float* space_weight = &_space_weight[0];
     int* space_ofs = &_space_ofs[0];
-
+    
     // initialize color-related bilateral filter coefficients
     for( i = 0; i < 256*cn; i++ )
         color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);
-
+    
     // initialize space-related bilateral filter coefficients
     for( i = -radius, maxk = 0; i <= radius; i++ )
         for( j = -radius; j <= radius; j++ )
@@ -1340,55 +1411,89 @@ bilateralFilter_8u( const Mat& src, Mat& dst, int d,
             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
             space_ofs[maxk++] = (int)(i*temp.step + j*cn);
         }
+    
+    BilateralFilter_8u_Invoker body(src, dst, temp, radius, maxk, space_ofs, space_weight, color_weight);
+    parallel_for_(Range(0, size.height), body);
+}
 
-    for( i = 0; i < size.height; i++ )
+
+class BilateralFilter_32f_Invoker :
+    public ParallelLoopBody
+{
+public:
+
+    BilateralFilter_32f_Invoker(int _cn, int _radius, int _maxk, int *_space_ofs,
+        Mat _temp, Mat *_dest, Size _size,
+        float _scale_index, float *_space_weight, float *_expLUT) :
+        ParallelLoopBody(), cn(_cn), radius(_radius), maxk(_maxk), space_ofs(_space_ofs),
+        temp(_temp), dest(_dest), size(_size), scale_index(_scale_index), space_weight(_space_weight), expLUT(_expLUT)
     {
-        const uchar* sptr = temp.data + (i+radius)*temp.step + radius*cn;
-        uchar* dptr = dst.data + i*dst.step;
+    }
 
-        if( cn == 1 )
+    virtual void operator() (const Range& range) const
+    {
+        Mat& dst = *dest;
+        int i, j, k;
+
+        for( i = range.start; i < range.end; i++ )
         {
-            for( j = 0; j < size.width; j++ )
+            const float* sptr = (const float*)(temp.data + (i+radius)*temp.step) + radius*cn;
+            float* dptr = (float*)(dst.data + i*dst.step);
+
+            if( cn == 1 )
             {
-                float sum = 0, wsum = 0;
-                int val0 = sptr[j];
-                for( k = 0; k < maxk; k++ )
+                for( j = 0; j < size.width; j++ )
                 {
-                    int val = sptr[j + space_ofs[k]];
-                    float w = space_weight[k]*color_weight[std::abs(val - val0)];
-                    sum += val*w;
-                    wsum += w;
+                    float sum = 0, wsum = 0;
+                    float val0 = sptr[j];
+                    for( k = 0; k < maxk; k++ )
+                    {
+                        float val = sptr[j + space_ofs[k]];
+                        float alpha = (float)(std::abs(val - val0)*scale_index);
+                        int idx = cvFloor(alpha);
+                        alpha -= idx;
+                        float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
+                        sum += val*w;
+                        wsum += w;
+                    }
+                    dptr[j] = (float)(sum/wsum);
                 }
-                // overflow is not possible here => there is no need to use CV_CAST_8U
-                dptr[j] = (uchar)cvRound(sum/wsum);
             }
-        }
-        else
-        {
-            assert( cn == 3 );
-            for( j = 0; j < size.width*3; j += 3 )
+            else
             {
-                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
-                int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
-                for( k = 0; k < maxk; k++ )
+                assert( cn == 3 );
+                for( j = 0; j < size.width*3; j += 3 )
                 {
-                    const uchar* sptr_k = sptr + j + space_ofs[k];
-                    int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
-                    float w = space_weight[k]*color_weight[std::abs(b - b0) +
-                        std::abs(g - g0) + std::abs(r - r0)];
-                    sum_b += b*w; sum_g += g*w; sum_r += r*w;
-                    wsum += w;
+                    float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
+                    float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
+                    for( k = 0; k < maxk; k++ )
+                    {
+                        const float* sptr_k = sptr + j + space_ofs[k];
+                        float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
+                        float alpha = (float)((std::abs(b - b0) +
+                            std::abs(g - g0) + std::abs(r - r0))*scale_index);
+                        int idx = cvFloor(alpha);
+                        alpha -= idx;
+                        float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
+                        sum_b += b*w; sum_g += g*w; sum_r += r*w;
+                        wsum += w;
+                    }
+                    wsum = 1.f/wsum;
+                    b0 = sum_b*wsum;
+                    g0 = sum_g*wsum;
+                    r0 = sum_r*wsum;
+                    dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
                 }
-                wsum = 1.f/wsum;
-                b0 = cvRound(sum_b*wsum);
-                g0 = cvRound(sum_g*wsum);
-                r0 = cvRound(sum_r*wsum);
-                dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
             }
         }
     }
-}
 
+private:
+    int cn, radius, maxk, *space_ofs;
+    Mat temp, *dest;
+    Size size;
+    float scale_index, *space_weight, *expLUT;
+};
 
 static void
 bilateralFilter_32f( const Mat& src, Mat& dst, int d,
@@ -1396,7 +1501,7 @@ bilateralFilter_32f( const Mat& src, Mat& dst, int d,
                      int borderType )
 {
     int cn = src.channels();
-    int i, j, k, maxk, radius;
+    int i, j, maxk, radius;
     double minValSrc=-1, maxValSrc=1;
     const int kExpNumBinsPerChannel = 1 << 12;
     int kExpNumBins = 0;
@@ -1474,57 +1579,10 @@ bilateralFilter_32f( const Mat& src, Mat& dst, int d,
             space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
         }
 
-    for( i = 0; i < size.height; i++ )
-    {
-        const float* sptr = (const float*)(temp.data + (i+radius)*temp.step) + radius*cn;
-        float* dptr = (float*)(dst.data + i*dst.step);
+    // parallel_for usage
 
-        if( cn == 1 )
-        {
-            for( j = 0; j < size.width; j++ )
-            {
-                float sum = 0, wsum = 0;
-                float val0 = sptr[j];
-                for( k = 0; k < maxk; k++ )
-                {
-                    float val = sptr[j + space_ofs[k]];
-                    float alpha = (float)(std::abs(val - val0)*scale_index);
-                    int idx = cvFloor(alpha);
-                    alpha -= idx;
-                    float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
-                    sum += val*w;
-                    wsum += w;
-                }
-                dptr[j] = (float)(sum/wsum);
-            }
-        }
-        else
-        {
-            assert( cn == 3 );
-            for( j = 0; j < size.width*3; j += 3 )
-            {
-                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
-                float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
-                for( k = 0; k < maxk; k++ )
-                {
-                    const float* sptr_k = sptr + j + space_ofs[k];
-                    float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
-                    float alpha = (float)((std::abs(b - b0) +
-                        std::abs(g - g0) + std::abs(r - r0))*scale_index);
-                    int idx = cvFloor(alpha);
-                    alpha -= idx;
-                    float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
-                    sum_b += b*w; sum_g += g*w; sum_r += r*w;
-                    wsum += w;
-                }
-                wsum = 1.f/wsum;
-                b0 = sum_b*wsum;
-                g0 = sum_g*wsum;
-                r0 = sum_r*wsum;
-                dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
-            }
-        }
-    }
+    BilateralFilter_32f_Invoker body(cn, radius, maxk, space_ofs, temp, &dst, size, scale_index, space_weight, expLUT);
+    parallel_for_(Range(0, size.height), body);
 }
 
 }

modules/imgproc/test/test_bilateral_filter.cpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "test_precomp.hpp"
+
+using namespace cv;
+
+namespace cvtest 
+{ 
+    class CV_BilateralFilterTest : 
+        public cvtest::BaseTest
+    {
+    public:
+        enum 
+        {
+            MAX_WIDTH = 1920, MIN_WIDTH = 1,
+            MAX_HEIGHT = 1080, MIN_HEIGHT = 1 
+        };
+
+        CV_BilateralFilterTest();
+        ~CV_BilateralFilterTest();
+
+    protected:
+        virtual void run_func();
+        virtual int prepare_test_case(int test_case_index);
+        virtual int validate_test_results(int test_case_index);
+        
+    private:
+        void reference_bilateral_filter(const Mat& src, Mat& dst, int d, double sigma_color,
+            double sigma_space, int borderType = BORDER_DEFAULT);
+            
+        int getRandInt(RNG& rng, int min_value, int max_value) const;
+                       
+        double _sigma_color;
+        double _sigma_space;
+        
+        Mat _src;
+        Mat _parallel_dst;
+        int _d;
+    };
+
+    CV_BilateralFilterTest::CV_BilateralFilterTest() : 
+        cvtest::BaseTest(), _src(), _parallel_dst(), _d()
+    {
+        test_case_count = 1000;
+    }
+
+    CV_BilateralFilterTest::~CV_BilateralFilterTest()
+    {
+    }
+
+    int CV_BilateralFilterTest::getRandInt(RNG& rng, int min_value, int max_value) const
+    {
+        double rand_value = rng.uniform(log(min_value), log(max_value + 1));
+        return cvRound(exp(rand_value));
+    }
+
+    void CV_BilateralFilterTest::reference_bilateral_filter(const Mat &src, Mat &dst, int d, 
+        double sigma_color, double sigma_space, int borderType)
+    {
+        int cn = src.channels();
+        int i, j, k, maxk, radius;
+        double minValSrc = -1, maxValSrc = 1;
+        const int kExpNumBinsPerChannel = 1 << 12;
+        int kExpNumBins = 0;
+        float lastExpVal = 1.f;
+        float len, scale_index;
+        Size size = src.size();
+
+        dst.create(size, src.type());
+
+        CV_Assert( (src.type() == CV_32FC1 || src.type() == CV_32FC3) &&
+            src.type() == dst.type() && src.size() == dst.size() &&
+            src.data != dst.data );
+
+        if( sigma_color <= 0 )
+            sigma_color = 1;
+        if( sigma_space <= 0 )
+            sigma_space = 1;
+
+        double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
+        double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
+
+        if( d <= 0 )
+            radius = cvRound(sigma_space*1.5);
+        else
+            radius = d/2;
+        radius = MAX(radius, 1);
+        d = radius*2 + 1;
+        // compute the min/max range for the input image (even if multichannel)
+
+        minMaxLoc( src.reshape(1), &minValSrc, &maxValSrc );
+        if(std::abs(minValSrc - maxValSrc) < FLT_EPSILON)
+        {
+            src.copyTo(dst);
+            return;
+        }
+
+        // temporary copy of the image with borders for easy processing
+        Mat temp;
+        copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
+        patchNaNs(temp);
+
+        // allocate lookup tables
+        vector<float> _space_weight(d*d);
+        vector<int> _space_ofs(d*d);
+        float* space_weight = &_space_weight[0];
+        int* space_ofs = &_space_ofs[0];
+
+        // assign a length which is slightly more than needed
+        len = (float)(maxValSrc - minValSrc) * cn;
+        kExpNumBins = kExpNumBinsPerChannel * cn;
+        vector<float> _expLUT(kExpNumBins+2);
+        float* expLUT = &_expLUT[0];
+
+        scale_index = kExpNumBins/len;
+
+        // initialize the exp LUT
+        for( i = 0; i < kExpNumBins+2; i++ )
+        {
+            if( lastExpVal > 0.f )
+            {
+                double val =  i / scale_index;
+                expLUT[i] = (float)std::exp(val * val * gauss_color_coeff);
+                lastExpVal = expLUT[i];
+            }
+            else
+                expLUT[i] = 0.f;
+        }
+
+        // initialize space-related bilateral filter coefficients
+        for( i = -radius, maxk = 0; i <= radius; i++ )
+            for( j = -radius; j <= radius; j++ )
+            {
+                double r = std::sqrt((double)i*i + (double)j*j);
+                if( r > radius )
+                    continue;
+                space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
+                space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
+            }
+
+        for( i = 0; i < size.height; i++ )
+        {
+            const float* sptr = (const float*)(temp.data + (i+radius)*temp.step) + radius*cn;
+            float* dptr = (float*)(dst.data + i*dst.step);
+
+            if( cn == 1 )
+            {
+                for( j = 0; j < size.width; j++ )
+                {
+                    float sum = 0, wsum = 0;
+                    float val0 = sptr[j];
+                    for( k = 0; k < maxk; k++ )
+                    {
+                        float val = sptr[j + space_ofs[k]];
+                        float alpha = (float)(std::abs(val - val0)*scale_index);
+                        int idx = cvFloor(alpha);
+                        alpha -= idx;
+                        float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
+                        sum += val*w;
+                        wsum += w;
+                    }
+                    dptr[j] = (float)(sum/wsum);
+                }
+            }
+            else
+            {
+                assert( cn == 3 );
+                for( j = 0; j < size.width*3; j += 3 )
+                {
+                    float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
+                    float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
+                    for( k = 0; k < maxk; k++ )
+                    {
+                        const float* sptr_k = sptr + j + space_ofs[k];
+                        float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
+                        float alpha = (float)((std::abs(b - b0) +
+                            std::abs(g - g0) + std::abs(r - r0))*scale_index);
+                        int idx = cvFloor(alpha);
+                        alpha -= idx;
+                        float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
+                        sum_b += b*w; sum_g += g*w; sum_r += r*w;
+                        wsum += w;
+                    }
+                    wsum = 1.f/wsum;
+                    b0 = sum_b*wsum;
+                    g0 = sum_g*wsum;
+                    r0 = sum_r*wsum;
+                    dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
+                }
+            }
+        }
+    }
+
+    int CV_BilateralFilterTest::prepare_test_case(int /* test_case_index */)
+    {
+        const static int types[] = { CV_32FC1, CV_32FC3, CV_8UC1, CV_8UC3 };
+        RNG& rng = ts->get_rng();
+        Size size(getRandInt(rng, MIN_WIDTH, MAX_WIDTH), getRandInt(rng, MIN_HEIGHT, MAX_HEIGHT));
+        int type = types[rng(sizeof(types) / sizeof(types[0]))];
+    
+        _d = rng.uniform(0., 1.) > 0.5 ? 5 : 3; 
+        
+        _src.create(size, type);
+        
+        rng.fill(_src, RNG::UNIFORM, 0, 256);
+        
+        _sigma_color = _sigma_space = 1.;
+        
+        return 1;
+    }
+
+    int CV_BilateralFilterTest::validate_test_results(int test_case_index)
+    {
+        static const double eps = 1;
+        
+        Mat reference_dst, reference_src;
+        if (_src.depth() == CV_32F)
+            reference_bilateral_filter(_src, reference_dst, _d, _sigma_color, _sigma_space);
+        else
+        {
+            int type = _src.type();
+            _src.convertTo(reference_src, CV_32F);
+            reference_bilateral_filter(reference_src, reference_dst, _d, _sigma_color, _sigma_space);
+            reference_dst.convertTo(reference_dst, type);
+        }
+        
+        double e = norm(reference_dst, _parallel_dst);
+        if (e > eps)
+        {
+            ts->printf(cvtest::TS::CONSOLE, "actual error: %g, expected: %g", e, eps);
+            ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
+        }
+        else
+            ts->set_failed_test_info(cvtest::TS::OK);
+        
+        return BaseTest::validate_test_results(test_case_index);
+    }
+
+    void CV_BilateralFilterTest::run_func()
+    {
+        bilateralFilter(_src, _parallel_dst, _d, _sigma_color, _sigma_space);
+    }
+
+    TEST(Imgproc_BilateralFilter, accuracy)
+    { 
+        CV_BilateralFilterTest test; 
+        test.safe_run(); 
+    }
+
+} // end of namespace cvtest