diff --git a/modules/core/include/opencv2/core/parallel_tool.hpp b/modules/core/include/opencv2/core/parallel_tool.hpp new file mode 100644 index 000000000..08258d5c2 --- /dev/null +++ b/modules/core/include/opencv2/core/parallel_tool.hpp @@ -0,0 +1,108 @@ +/*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 +#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 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 inline + CV_EXPORTS void parallel_reduce_(const Range& range, Body& body) + { +#ifdef HAVE_TBB + tbb::parallel_reduce(tbb::blocked_range(range.start, range.end), body); +#else + body(range); +#endif // end HAVE_TBB + } + +} // namespace cv + +#endif // __cplusplus + +#endif // __OPENCV_PARALLEL_TOOL_HPP__ diff --git a/modules/core/src/parallel_tool.cpp b/modules/core/src/parallel_tool.cpp new file mode 100644 index 000000000..423d4787d --- /dev/null +++ b/modules/core/src/parallel_tool.cpp @@ -0,0 +1,112 @@ +/*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 +#elif defined HAVE_OPENMP +# include +#elif defined HAVE_GCD +# include +#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& 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(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(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(body), block_function); + +#else + + body(range); + +#endif // end HAVE_TBB + } + +} // namespace cv diff --git a/modules/core/src/precomp.hpp b/modules/core/src/precomp.hpp index 81b9d6e80..60429075a 100644 --- a/modules/core/src/precomp.hpp +++ b/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 #include diff --git a/modules/imgproc/perf/perf_bilateral.cpp b/modules/imgproc/perf/perf_bilateral.cpp new file mode 100644 index 000000000..85cfc7d0c --- /dev/null +++ b/modules/imgproc/perf/perf_bilateral.cpp @@ -0,0 +1,38 @@ +#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 > 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); +} diff --git a/modules/imgproc/src/precomp.hpp b/modules/imgproc/src/precomp.hpp index fef5f755b..998008ae2 100644 --- a/modules/imgproc/src/precomp.hpp +++ b/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 #include #include diff --git a/modules/imgproc/src/smooth.cpp b/modules/imgproc/src/smooth.cpp index faec530b8..1bc11c7fc 100644 --- a/modules/imgproc/src/smooth.cpp +++ b/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 _color_weight(cn*256); vector _space_weight(d*d); vector _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); } } diff --git a/modules/imgproc/test/test_bilateral_filter.cpp b/modules/imgproc/test/test_bilateral_filter.cpp new file mode 100644 index 000000000..034f9c363 --- /dev/null +++ b/modules/imgproc/test/test_bilateral_filter.cpp @@ -0,0 +1,290 @@ +/*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 _space_weight(d*d); + vector _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 _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