1/*M///////////////////////////////////////////////////////////////////////////////////////
2//
3// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4//
5// By downloading, copying, installing or using the software you agree to this license.
6// If you do not agree to this license, do not download, install,
7// copy or use the software.
8//
9//
10// License Agreement
11// For Open Source Computer Vision Library
12//
13// Copyright (C) 2000-2008, 2018, Intel Corporation, all rights reserved.
14// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
15// Copyright (C) 2014-2015, Itseez Inc., all rights reserved.
16// Third party copyrights are property of their respective owners.
17//
18// Redistribution and use in source and binary forms, with or without modification,
19// are permitted provided that the following conditions are met:
20//
21// * Redistribution's of source code must retain the above copyright notice,
22// this list of conditions and the following disclaimer.
23//
24// * Redistribution's in binary form must reproduce the above copyright notice,
25// this list of conditions and the following disclaimer in the documentation
26// and/or other materials provided with the distribution.
27//
28// * The name of the copyright holders may not be used to endorse or promote products
29// derived from this software without specific prior written permission.
30//
31// This software is provided by the copyright holders and contributors "as is" and
32// any express or implied warranties, including, but not limited to, the implied
33// warranties of merchantability and fitness for a particular purpose are disclaimed.
34// In no event shall the Intel Corporation or contributors be liable for any direct,
35// indirect, incidental, special, exemplary, or consequential damages
36// (including, but not limited to, procurement of substitute goods or services;
37// loss of use, data, or profits; or business interruption) however caused
38// and on any theory of liability, whether in contract, strict liability,
39// or tort (including negligence or otherwise) arising in any way out of
40// the use of this software, even if advised of the possibility of such damage.
41//
42//M*/
43
44#include "precomp.hpp"
45
46#include <vector>
47
48#include "opencv2/core/hal/intrin.hpp"
49#include "opencl_kernels_imgproc.hpp"
50
51#include "bilateral_filter.simd.hpp"
52#include "bilateral_filter.simd_declarations.hpp" // defines CV_CPU_DISPATCH_MODES_ALL=AVX2,...,BASELINE based on CMakeLists.txt content
53
54/****************************************************************************************\
55 Bilateral Filtering
56\****************************************************************************************/
57
58namespace cv {
59
60#ifdef HAVE_OPENCL
61
62static bool ocl_bilateralFilter_8u(InputArray _src, OutputArray _dst, int d,
63 double sigma_color, double sigma_space,
64 int borderType)
65{
66 CV_INSTRUMENT_REGION();
67#ifdef __ANDROID__
68 if (ocl::Device::getDefault().isNVidia())
69 return false;
70#endif
71
72 int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
73 int i, j, maxk, radius;
74
75 if (depth != CV_8U || cn > 4)
76 return false;
77
78 constexpr double eps = 1e-6;
79 if( sigma_color <= eps || sigma_space <= eps )
80 {
81 _src.copyTo(arr: _dst);
82 return true;
83 }
84
85 double gauss_color_coeff = -0.5 / (sigma_color * sigma_color);
86 double gauss_space_coeff = -0.5 / (sigma_space * sigma_space);
87
88 if ( d <= 0 )
89 radius = cvRound(value: sigma_space * 1.5);
90 else
91 radius = d / 2;
92 radius = MAX(radius, 1);
93 d = radius * 2 + 1;
94
95 UMat src = _src.getUMat(), dst = _dst.getUMat(), temp;
96 if (src.u == dst.u)
97 return false;
98
99 copyMakeBorder(src, dst: temp, top: radius, bottom: radius, left: radius, right: radius, borderType);
100 std::vector<float> _space_weight(d * d);
101 std::vector<int> _space_ofs(d * d);
102 float * const space_weight = &_space_weight[0];
103 int * const space_ofs = &_space_ofs[0];
104
105 // initialize space-related bilateral filter coefficients
106 for( i = -radius, maxk = 0; i <= radius; i++ )
107 for( j = -radius; j <= radius; j++ )
108 {
109 double r = std::sqrt(x: (double)i * i + (double)j * j);
110 if ( r > radius )
111 continue;
112 space_weight[maxk] = (float)std::exp(x: r * r * gauss_space_coeff);
113 space_ofs[maxk++] = (int)(i * temp.step + j * cn);
114 }
115
116 char cvt[3][50];
117 String cnstr = cn > 1 ? format(fmt: "%d", cn) : "";
118 String kernelName("bilateral");
119 size_t sizeDiv = 1;
120 if ((ocl::Device::getDefault().isIntel()) &&
121 (ocl::Device::getDefault().type() == ocl::Device::TYPE_GPU))
122 {
123 //Intel GPU
124 if (dst.cols % 4 == 0 && cn == 1) // For single channel x4 sized images.
125 {
126 kernelName = "bilateral_float4";
127 sizeDiv = 4;
128 }
129 }
130 ocl::Kernel k(kernelName.c_str(), ocl::imgproc::bilateral_oclsrc,
131 format("-D radius=%d -D maxk=%d -D cn=%d -D int_t=%s -D uint_t=uint%s -D convert_int_t=%s"
132 " -D uchar_t=%s -D float_t=%s -D convert_float_t=%s -D convert_uchar_t=%s -D gauss_color_coeff=(float)%f",
133 radius, maxk, cn, ocl::typeToStr(CV_32SC(cn)), cnstr.c_str(),
134 ocl::convertTypeStr(CV_8U, CV_32S, cn, cvt[0], sizeof(cvt[0])),
135 ocl::typeToStr(type), ocl::typeToStr(CV_32FC(cn)),
136 ocl::convertTypeStr(CV_32S, CV_32F, cn, cvt[1], sizeof(cvt[1])),
137 ocl::convertTypeStr(CV_32F, CV_8U, cn, cvt[2], sizeof(cvt[2])), gauss_color_coeff));
138 if (k.empty())
139 return false;
140
141 Mat mspace_weight(1, d * d, CV_32FC1, space_weight);
142 Mat mspace_ofs(1, d * d, CV_32SC1, space_ofs);
143 UMat ucolor_weight, uspace_weight, uspace_ofs;
144
145 mspace_weight.copyTo(m: uspace_weight);
146 mspace_ofs.copyTo(m: uspace_ofs);
147
148 k.args(kernel_args: ocl::KernelArg::ReadOnlyNoSize(m: temp), kernel_args: ocl::KernelArg::WriteOnly(m: dst),
149 kernel_args: ocl::KernelArg::PtrReadOnly(m: uspace_weight),
150 kernel_args: ocl::KernelArg::PtrReadOnly(m: uspace_ofs));
151
152 size_t globalsize[2] = { (size_t)dst.cols / sizeDiv, (size_t)dst.rows };
153 return k.run(dims: 2, globalsize, NULL, sync: false);
154}
155#endif
156
157
158static void
159bilateralFilter_8u( const Mat& src, Mat& dst, int d,
160 double sigma_color, double sigma_space,
161 int borderType )
162{
163 CV_INSTRUMENT_REGION();
164
165 int cn = src.channels();
166 int i, j, maxk, radius;
167
168 CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) && src.data != dst.data );
169
170 constexpr double eps = 1e-6;
171 if( sigma_color <= eps || sigma_space <= eps )
172 {
173 src.copyTo(m: dst);
174 return;
175 }
176
177 double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
178 double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
179
180 if( d <= 0 )
181 radius = cvRound(value: sigma_space*1.5);
182 else
183 radius = d/2;
184 radius = MAX(radius, 1);
185 d = radius*2 + 1;
186
187 Mat temp;
188 copyMakeBorder( src, dst: temp, top: radius, bottom: radius, left: radius, right: radius, borderType );
189
190 std::vector<float> _color_weight(cn*256);
191 std::vector<float> _space_weight(d*d);
192 std::vector<int> _space_ofs(d*d);
193 float* color_weight = &_color_weight[0];
194 float* space_weight = &_space_weight[0];
195 int* space_ofs = &_space_ofs[0];
196
197 // initialize color-related bilateral filter coefficients
198
199 for( i = 0; i < 256*cn; i++ )
200 color_weight[i] = (float)std::exp(x: i*i*gauss_color_coeff);
201
202 // initialize space-related bilateral filter coefficients
203 for( i = -radius, maxk = 0; i <= radius; i++ )
204 {
205 j = -radius;
206
207 for( ; j <= radius; j++ )
208 {
209 double r = std::sqrt(x: (double)i*i + (double)j*j);
210 if( r > radius )
211 continue;
212 space_weight[maxk] = (float)std::exp(x: r*r*gauss_space_coeff);
213 space_ofs[maxk++] = (int)(i*temp.step + j*cn);
214 }
215 }
216
217 CV_CPU_DISPATCH(bilateralFilterInvoker_8u, (dst, temp, radius, maxk, space_ofs, space_weight, color_weight),
218 CV_CPU_DISPATCH_MODES_ALL);
219}
220
221
222static void
223bilateralFilter_32f( const Mat& src, Mat& dst, int d,
224 double sigma_color, double sigma_space,
225 int borderType )
226{
227 CV_INSTRUMENT_REGION();
228
229 int cn = src.channels();
230 int i, j, maxk, radius;
231 double minValSrc=-1, maxValSrc=1;
232 const int kExpNumBinsPerChannel = 1 << 12;
233 int kExpNumBins = 0;
234 float lastExpVal = 1.f;
235 float len, scale_index;
236
237 CV_Assert( (src.type() == CV_32FC1 || src.type() == CV_32FC3) && src.data != dst.data );
238
239 constexpr double eps = 1e-6;
240 if( sigma_color <= eps || sigma_space <= eps )
241 {
242 src.copyTo(m: dst);
243 return;
244 }
245
246 double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
247 double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
248
249 if( d <= 0 )
250 radius = cvRound(value: sigma_space*1.5);
251 else
252 radius = d/2;
253 radius = MAX(radius, 1);
254 d = radius*2 + 1;
255 // compute the min/max range for the input image (even if multichannel)
256
257 minMaxLoc( src: src.reshape(cn: 1), minVal: &minValSrc, maxVal: &maxValSrc );
258 if(std::abs(x: minValSrc - maxValSrc) < FLT_EPSILON)
259 {
260 src.copyTo(m: dst);
261 return;
262 }
263
264 // temporary copy of the image with borders for easy processing
265 Mat temp;
266 copyMakeBorder( src, dst: temp, top: radius, bottom: radius, left: radius, right: radius, borderType );
267
268 // allocate lookup tables
269 std::vector<float> _space_weight(d*d);
270 std::vector<int> _space_ofs(d*d);
271 float* space_weight = &_space_weight[0];
272 int* space_ofs = &_space_ofs[0];
273
274 // assign a length which is slightly more than needed
275 len = (float)(maxValSrc - minValSrc) * cn;
276 kExpNumBins = kExpNumBinsPerChannel * cn;
277 std::vector<float> _expLUT(kExpNumBins+2);
278 float* expLUT = &_expLUT[0];
279
280 scale_index = kExpNumBins/len;
281
282 // initialize the exp LUT
283 for( i = 0; i < kExpNumBins+2; i++ )
284 {
285 if( lastExpVal > 0.f )
286 {
287 double val = i / scale_index;
288 expLUT[i] = (float)std::exp(x: val * val * gauss_color_coeff);
289 lastExpVal = expLUT[i];
290 }
291 else
292 expLUT[i] = 0.f;
293 }
294
295 // initialize space-related bilateral filter coefficients
296 for( i = -radius, maxk = 0; i <= radius; i++ )
297 for( j = -radius; j <= radius; j++ )
298 {
299 double r = std::sqrt(x: (double)i*i + (double)j*j);
300 if( r > radius || ( i == 0 && j == 0 ) )
301 continue;
302 space_weight[maxk] = (float)std::exp(x: r*r*gauss_space_coeff);
303 space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
304 }
305
306 // parallel_for usage
307 CV_CPU_DISPATCH(bilateralFilterInvoker_32f, (cn, radius, maxk, space_ofs, temp, dst, scale_index, space_weight, expLUT),
308 CV_CPU_DISPATCH_MODES_ALL);
309}
310
311#ifdef HAVE_IPP
312#define IPP_BILATERAL_PARALLEL 1
313
314#ifdef HAVE_IPP_IW
315class ipp_bilateralFilterParallel: public ParallelLoopBody
316{
317public:
318 ipp_bilateralFilterParallel(::ipp::IwiImage &_src, ::ipp::IwiImage &_dst, int _radius, Ipp32f _valSquareSigma, Ipp32f _posSquareSigma, ::ipp::IwiBorderType _borderType, bool *_ok):
319 src(_src), dst(_dst)
320 {
321 pOk = _ok;
322
323 radius = _radius;
324 valSquareSigma = _valSquareSigma;
325 posSquareSigma = _posSquareSigma;
326 borderType = _borderType;
327
328 *pOk = true;
329 }
330 ~ipp_bilateralFilterParallel() {}
331
332 virtual void operator() (const Range& range) const CV_OVERRIDE
333 {
334 if(*pOk == false)
335 return;
336
337 try
338 {
339 ::ipp::IwiTile tile = ::ipp::IwiRoi(0, range.start, dst.m_size.width, range.end - range.start);
340 CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBilateral, src, dst, radius, valSquareSigma, posSquareSigma, ::ipp::IwDefault(), borderType, tile);
341 }
342 catch(const ::ipp::IwException &)
343 {
344 *pOk = false;
345 return;
346 }
347 }
348private:
349 ::ipp::IwiImage &src;
350 ::ipp::IwiImage &dst;
351
352 int radius;
353 Ipp32f valSquareSigma;
354 Ipp32f posSquareSigma;
355 ::ipp::IwiBorderType borderType;
356
357 bool *pOk;
358 const ipp_bilateralFilterParallel& operator= (const ipp_bilateralFilterParallel&);
359};
360#endif
361
362static bool ipp_bilateralFilter(Mat &src, Mat &dst, int d, double sigmaColor, double sigmaSpace, int borderType)
363{
364#ifdef HAVE_IPP_IW
365 CV_INSTRUMENT_REGION_IPP();
366
367 constexpr double eps = 1e-6;
368 if( sigmaColor <= eps || sigmaSpace <= eps )
369 {
370 src.copyTo(m: dst);
371 return true;
372 }
373
374 int radius = IPP_MAX(((d <= 0)?cvRound(sigmaSpace*1.5):d/2), 1);
375 Ipp32f valSquareSigma = (Ipp32f)(sigmaColor*sigmaColor);
376 Ipp32f posSquareSigma = (Ipp32f)(sigmaSpace*sigmaSpace);
377
378 // Acquire data and begin processing
379 try
380 {
381 ::ipp::IwiImage iwSrc = ippiGetImage(src);
382 ::ipp::IwiImage iwDst = ippiGetImage(src: dst);
383 ::ipp::IwiBorderSize borderSize(radius);
384 ::ipp::IwiBorderType ippBorder(ippiGetBorder(image&: iwSrc, ocvBorderType: borderType, borderSize));
385 if(!ippBorder)
386 return false;
387
388 const int threads = ippiSuggestThreadsNum(image: iwDst, multiplier: 2);
389 if(IPP_BILATERAL_PARALLEL && threads > 1) {
390 bool ok = true;
391 Range range(0, (int)iwDst.m_size.height);
392 ipp_bilateralFilterParallel invoker(iwSrc, iwDst, radius, valSquareSigma, posSquareSigma, ippBorder, &ok);
393 if(!ok)
394 return false;
395
396 parallel_for_(range, body: invoker, nstripes: threads*4);
397
398 if(!ok)
399 return false;
400 } else {
401 CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBilateral, iwSrc, iwDst, radius, valSquareSigma, posSquareSigma, ::ipp::IwDefault(), ippBorder);
402 }
403 }
404 catch (const ::ipp::IwException &)
405 {
406 return false;
407 }
408 return true;
409#else
410 CV_UNUSED(src); CV_UNUSED(dst); CV_UNUSED(d); CV_UNUSED(sigmaColor); CV_UNUSED(sigmaSpace); CV_UNUSED(borderType);
411 return false;
412#endif
413}
414#endif
415
416void bilateralFilter( InputArray _src, OutputArray _dst, int d,
417 double sigmaColor, double sigmaSpace,
418 int borderType )
419{
420 CV_INSTRUMENT_REGION();
421
422 CV_Assert(!_src.empty());
423
424 _dst.create( sz: _src.size(), type: _src.type() );
425
426 CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
427 ocl_bilateralFilter_8u(_src, _dst, d, sigma_color: sigmaColor, sigma_space: sigmaSpace, borderType))
428
429 Mat src = _src.getMat(), dst = _dst.getMat();
430
431 CALL_HAL(bilateralFilter, cv_hal_bilateralFilter, src.data, src.step, dst.data, dst.step, src.cols, src.rows, src.depth(),
432 src.channels(), d, sigmaColor, sigmaSpace, borderType);
433
434 CV_IPP_RUN_FAST(ipp_bilateralFilter(src, dst, d, sigmaColor, sigmaSpace, borderType));
435
436 if( src.depth() == CV_8U )
437 bilateralFilter_8u( src, dst, d, sigma_color: sigmaColor, sigma_space: sigmaSpace, borderType );
438 else if( src.depth() == CV_32F )
439 bilateralFilter_32f( src, dst, d, sigma_color: sigmaColor, sigma_space: sigmaSpace, borderType );
440 else
441 CV_Error( cv::Error::StsUnsupportedFormat,
442 "Bilateral filtering is only implemented for 8u and 32f images" );
443}
444
445} // namespace
446

source code of opencv/modules/imgproc/src/bilateral_filter.dispatch.cpp