modules/imgproc/src/smooth.cpp

   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, 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 #include "opencv2/core/hal/intrin.hpp"
  46 #include "opencl_kernels_imgproc.hpp"
  47
  48 #include "opencv2/core/openvx/ovx_defs.hpp"
  49
  50 /*
  51  * This file includes the code, contributed by Simon Perreault
  52  * (the function icvMedianBlur_8u_O1)
  53  *
  54  * Constant-time median filtering -- http://nomis80.org/ctmf.html
  55  * Copyright (C) 2006 Simon Perreault
  56  *
  57  * Contact:
  58  *  Laboratoire de vision et systemes numeriques
  59  *  Pavillon Adrien-Pouliot
  60  *  Universite Laval
  61  *  Sainte-Foy, Quebec, Canada
  62  *  G1K 7P4
  63  *
  64  *  perreaul@gel.ulaval.ca
  65  */
  66
  67 namespace cv
  68 {
  69
  70 /****************************************************************************************\
  71                                          Box Filter
  72 \****************************************************************************************/
  73
  74 template<typename T, typename ST>
  75 struct RowSum :
  76         public BaseRowFilter
  77 {
  78     RowSum( int _ksize, int _anchor ) :
  79         BaseRowFilter()
  80     {
  81         ksize = _ksize;
  82         anchor = _anchor;
  83     }
  84
  85     virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
  86     {
  87         const T* S = (const T*)src;
  88         ST* D = (ST*)dst;
  89         int i = 0, k, ksz_cn = ksize*cn;
  90
  91         width = (width - 1)*cn;
  92         if( ksize == 3 )
  93         {
  94             for( i = 0; i < width + cn; i++ )
  95             {
  96                 D[i] = (ST)S[i] + (ST)S[i+cn] + (ST)S[i+cn*2];
  97             }
  98         }
  99         else if( ksize == 5 )
 100         {
 101             for( i = 0; i < width + cn; i++ )
 102             {
 103                 D[i] = (ST)S[i] + (ST)S[i+cn] + (ST)S[i+cn*2] + (ST)S[i + cn*3] + (ST)S[i + cn*4];
 104             }
 105         }
 106         else if( cn == 1 )
 107         {
 108             ST s = 0;
 109             for( i = 0; i < ksz_cn; i++ )
 110                 s += (ST)S[i];
 111             D[0] = s;
 112             for( i = 0; i < width; i++ )
 113             {
 114                 s += (ST)S[i + ksz_cn] - (ST)S[i];
 115                 D[i+1] = s;
 116             }
 117         }
 118         else if( cn == 3 )
 119         {
 120             ST s0 = 0, s1 = 0, s2 = 0;
 121             for( i = 0; i < ksz_cn; i += 3 )
 122             {
 123                 s0 += (ST)S[i];
 124                 s1 += (ST)S[i+1];
 125                 s2 += (ST)S[i+2];
 126             }
 127             D[0] = s0;
 128             D[1] = s1;
 129             D[2] = s2;
 130             for( i = 0; i < width; i += 3 )
 131             {
 132                 s0 += (ST)S[i + ksz_cn] - (ST)S[i];
 133                 s1 += (ST)S[i + ksz_cn + 1] - (ST)S[i + 1];
 134                 s2 += (ST)S[i + ksz_cn + 2] - (ST)S[i + 2];
 135                 D[i+3] = s0;
 136                 D[i+4] = s1;
 137                 D[i+5] = s2;
 138             }
 139         }
 140         else if( cn == 4 )
 141         {
 142             ST s0 = 0, s1 = 0, s2 = 0, s3 = 0;
 143             for( i = 0; i < ksz_cn; i += 4 )
 144             {
 145                 s0 += (ST)S[i];
 146                 s1 += (ST)S[i+1];
 147                 s2 += (ST)S[i+2];
 148                 s3 += (ST)S[i+3];
 149             }
 150             D[0] = s0;
 151             D[1] = s1;
 152             D[2] = s2;
 153             D[3] = s3;
 154             for( i = 0; i < width; i += 4 )
 155             {
 156                 s0 += (ST)S[i + ksz_cn] - (ST)S[i];
 157                 s1 += (ST)S[i + ksz_cn + 1] - (ST)S[i + 1];
 158                 s2 += (ST)S[i + ksz_cn + 2] - (ST)S[i + 2];
 159                 s3 += (ST)S[i + ksz_cn + 3] - (ST)S[i + 3];
 160                 D[i+4] = s0;
 161                 D[i+5] = s1;
 162                 D[i+6] = s2;
 163                 D[i+7] = s3;
 164             }
 165         }
 166         else
 167             for( k = 0; k < cn; k++, S++, D++ )
 168             {
 169                 ST s = 0;
 170                 for( i = 0; i < ksz_cn; i += cn )
 171                     s += (ST)S[i];
 172                 D[0] = s;
 173                 for( i = 0; i < width; i += cn )
 174                 {
 175                     s += (ST)S[i + ksz_cn] - (ST)S[i];
 176                     D[i+cn] = s;
 177                 }
 178             }
 179     }
 180 };
 181
 182
 183 template<typename ST, typename T>
 184 struct ColumnSum :
 185         public BaseColumnFilter
 186 {
 187     ColumnSum( int _ksize, int _anchor, double _scale ) :
 188         BaseColumnFilter()
 189     {
 190         ksize = _ksize;
 191         anchor = _anchor;
 192         scale = _scale;
 193         sumCount = 0;
 194     }
 195
 196     virtual void reset() { sumCount = 0; }
 197
 198     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 199     {
 200         int i;
 201         ST* SUM;
 202         bool haveScale = scale != 1;
 203         double _scale = scale;
 204
 205         if( width != (int)sum.size() )
 206         {
 207             sum.resize(width);
 208             sumCount = 0;
 209         }
 210
 211         SUM = &sum[0];
 212         if( sumCount == 0 )
 213         {
 214             memset((void*)SUM, 0, width*sizeof(ST));
 215
 216             for( ; sumCount < ksize - 1; sumCount++, src++ )
 217             {
 218                 const ST* Sp = (const ST*)src[0];
 219
 220                 for( i = 0; i < width; i++ )
 221                     SUM[i] += Sp[i];
 222             }
 223         }
 224         else
 225         {
 226             CV_Assert( sumCount == ksize-1 );
 227             src += ksize-1;
 228         }
 229
 230         for( ; count--; src++ )
 231         {
 232             const ST* Sp = (const ST*)src[0];
 233             const ST* Sm = (const ST*)src[1-ksize];
 234             T* D = (T*)dst;
 235             if( haveScale )
 236             {
 237                 for( i = 0; i <= width - 2; i += 2 )
 238                 {
 239                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 240                     D[i] = saturate_cast<T>(s0*_scale);
 241                     D[i+1] = saturate_cast<T>(s1*_scale);
 242                     s0 -= Sm[i]; s1 -= Sm[i+1];
 243                     SUM[i] = s0; SUM[i+1] = s1;
 244                 }
 245
 246                 for( ; i < width; i++ )
 247                 {
 248                     ST s0 = SUM[i] + Sp[i];
 249                     D[i] = saturate_cast<T>(s0*_scale);
 250                     SUM[i] = s0 - Sm[i];
 251                 }
 252             }
 253             else
 254             {
 255                 for( i = 0; i <= width - 2; i += 2 )
 256                 {
 257                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 258                     D[i] = saturate_cast<T>(s0);
 259                     D[i+1] = saturate_cast<T>(s1);
 260                     s0 -= Sm[i]; s1 -= Sm[i+1];
 261                     SUM[i] = s0; SUM[i+1] = s1;
 262                 }
 263
 264                 for( ; i < width; i++ )
 265                 {
 266                     ST s0 = SUM[i] + Sp[i];
 267                     D[i] = saturate_cast<T>(s0);
 268                     SUM[i] = s0 - Sm[i];
 269                 }
 270             }
 271             dst += dststep;
 272         }
 273     }
 274
 275     double scale;
 276     int sumCount;
 277     std::vector<ST> sum;
 278 };
 279
 280
 281 template<>
 282 struct ColumnSum<int, uchar> :
 283         public BaseColumnFilter
 284 {
 285     ColumnSum( int _ksize, int _anchor, double _scale ) :
 286         BaseColumnFilter()
 287     {
 288         ksize = _ksize;
 289         anchor = _anchor;
 290         scale = _scale;
 291         sumCount = 0;
 292     }
 293
 294     virtual void reset() { sumCount = 0; }
 295
 296     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 297     {
 298         int* SUM;
 299         bool haveScale = scale != 1;
 300         double _scale = scale;
 301
 302 #if CV_SIMD128
 303         bool haveSIMD128 = hasSIMD128();
 304 #endif
 305
 306         if( width != (int)sum.size() )
 307         {
 308             sum.resize(width);
 309             sumCount = 0;
 310         }
 311
 312         SUM = &sum[0];
 313         if( sumCount == 0 )
 314         {
 315             memset((void*)SUM, 0, width*sizeof(int));
 316             for( ; sumCount < ksize - 1; sumCount++, src++ )
 317             {
 318                 const int* Sp = (const int*)src[0];
 319                 int i = 0;
 320 #if CV_SIMD128
 321                 if( haveSIMD128 )
 322                 {
 323                     for (; i <= width - 4; i += 4)
 324                     {
 325                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 326                     }
 327                 }
 328 #endif
 329                 for( ; i < width; i++ )
 330                     SUM[i] += Sp[i];
 331             }
 332         }
 333         else
 334         {
 335             CV_Assert( sumCount == ksize-1 );
 336             src += ksize-1;
 337         }
 338
 339         for( ; count--; src++ )
 340         {
 341             const int* Sp = (const int*)src[0];
 342             const int* Sm = (const int*)src[1-ksize];
 343             uchar* D = (uchar*)dst;
 344             if( haveScale )
 345             {
 346                 int i = 0;
 347 #if CV_SIMD128
 348                 if( haveSIMD128 )
 349                 {
 350
 351                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 352                     for( ; i <= width-8; i+=8 )
 353                     {
 354                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 355                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 356
 357                         v_uint32x4 v_s0d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s0) * v_scale));
 358                         v_uint32x4 v_s01d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s01) * v_scale));
 359
 360                         v_uint16x8 v_dst = v_pack(v_s0d, v_s01d);
 361                         v_pack_store(D + i, v_dst);
 362
 363                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 364                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 365                     }
 366                 }
 367 #endif
 368                 for( ; i < width; i++ )
 369                 {
 370                     int s0 = SUM[i] + Sp[i];
 371                     D[i] = saturate_cast<uchar>(s0*_scale);
 372                     SUM[i] = s0 - Sm[i];
 373                 }
 374             }
 375             else
 376             {
 377                 int i = 0;
 378 #if CV_SIMD128
 379                 if( haveSIMD128 )
 380                 {
 381                     for( ; i <= width-8; i+=8 )
 382                     {
 383                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 384                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 385
 386                         v_uint16x8 v_dst = v_pack(v_reinterpret_as_u32(v_s0), v_reinterpret_as_u32(v_s01));
 387                         v_pack_store(D + i, v_dst);
 388
 389                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 390                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 391                     }
 392                 }
 393 #endif
 394
 395                 for( ; i < width; i++ )
 396                 {
 397                     int s0 = SUM[i] + Sp[i];
 398                     D[i] = saturate_cast<uchar>(s0);
 399                     SUM[i] = s0 - Sm[i];
 400                 }
 401             }
 402             dst += dststep;
 403         }
 404     }
 405
 406     double scale;
 407     int sumCount;
 408     std::vector<int> sum;
 409 };
 410
 411
 412 template<>
 413 struct ColumnSum<ushort, uchar> :
 414 public BaseColumnFilter
 415 {
 416     enum { SHIFT = 23 };
 417
 418     ColumnSum( int _ksize, int _anchor, double _scale ) :
 419     BaseColumnFilter()
 420     {
 421         ksize = _ksize;
 422         anchor = _anchor;
 423         scale = _scale;
 424         sumCount = 0;
 425         divDelta = 0;
 426         divScale = 1;
 427         if( scale != 1 )
 428         {
 429             int d = cvRound(1./scale);
 430             double scalef = ((double)(1 << SHIFT))/d;
 431             divScale = cvFloor(scalef);
 432             scalef -= divScale;
 433             divDelta = d/2;
 434             if( scalef < 0.5 )
 435                 divDelta++;
 436             else
 437                 divScale++;
 438         }
 439     }
 440
 441     virtual void reset() { sumCount = 0; }
 442
 443     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 444     {
 445         const int ds = divScale;
 446         const int dd = divDelta;
 447         ushort* SUM;
 448         const bool haveScale = scale != 1;
 449
 450 #if CV_SIMD128
 451         bool haveSIMD128 = hasSIMD128();
 452 #endif
 453
 454         if( width != (int)sum.size() )
 455         {
 456             sum.resize(width);
 457             sumCount = 0;
 458         }
 459
 460         SUM = &sum[0];
 461         if( sumCount == 0 )
 462         {
 463             memset((void*)SUM, 0, width*sizeof(SUM[0]));
 464             for( ; sumCount < ksize - 1; sumCount++, src++ )
 465             {
 466                 const ushort* Sp = (const ushort*)src[0];
 467                 int i = 0;
 468 #if CV_SIMD128
 469                 if( haveSIMD128 )
 470                 {
 471                     for( ; i <= width - 8; i += 8 )
 472                     {
 473                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 474                     }
 475                 }
 476 #endif
 477                 for( ; i < width; i++ )
 478                     SUM[i] += Sp[i];
 479             }
 480         }
 481         else
 482         {
 483             CV_Assert( sumCount == ksize-1 );
 484             src += ksize-1;
 485         }
 486
 487         for( ; count--; src++ )
 488         {
 489             const ushort* Sp = (const ushort*)src[0];
 490             const ushort* Sm = (const ushort*)src[1-ksize];
 491             uchar* D = (uchar*)dst;
 492             if( haveScale )
 493             {
 494                 int i = 0;
 495             #if CV_SIMD128
 496                 v_uint32x4 ds4 = v_setall_u32((unsigned)ds);
 497                 v_uint16x8 dd8 = v_setall_u16((ushort)dd);
 498
 499                 for( ; i <= width-16; i+=16 )
 500                 {
 501                     v_uint16x8 _sm0 = v_load(Sm + i);
 502                     v_uint16x8 _sm1 = v_load(Sm + i + 8);
 503
 504                     v_uint16x8 _s0 = v_add_wrap(v_load(SUM + i), v_load(Sp + i));
 505                     v_uint16x8 _s1 = v_add_wrap(v_load(SUM + i + 8), v_load(Sp + i + 8));
 506
 507                     v_uint32x4 _s00, _s01, _s10, _s11;
 508
 509                     v_expand(_s0 + dd8, _s00, _s01);
 510                     v_expand(_s1 + dd8, _s10, _s11);
 511
 512                     _s00 = v_shr<SHIFT>(_s00*ds4);
 513                     _s01 = v_shr<SHIFT>(_s01*ds4);
 514                     _s10 = v_shr<SHIFT>(_s10*ds4);
 515                     _s11 = v_shr<SHIFT>(_s11*ds4);
 516
 517                     v_int16x8 r0 = v_pack(v_reinterpret_as_s32(_s00), v_reinterpret_as_s32(_s01));
 518                     v_int16x8 r1 = v_pack(v_reinterpret_as_s32(_s10), v_reinterpret_as_s32(_s11));
 519
 520                     _s0 = v_sub_wrap(_s0, _sm0);
 521                     _s1 = v_sub_wrap(_s1, _sm1);
 522
 523                     v_store(D + i, v_pack_u(r0, r1));
 524                     v_store(SUM + i, _s0);
 525                     v_store(SUM + i + 8, _s1);
 526                 }
 527             #endif
 528                 for( ; i < width; i++ )
 529                 {
 530                     int s0 = SUM[i] + Sp[i];
 531                     D[i] = (uchar)((s0 + dd)*ds >> SHIFT);
 532                     SUM[i] = (ushort)(s0 - Sm[i]);
 533                 }
 534             }
 535             else
 536             {
 537                 int i = 0;
 538                 for( ; i < width; i++ )
 539                 {
 540                     int s0 = SUM[i] + Sp[i];
 541                     D[i] = saturate_cast<uchar>(s0);
 542                     SUM[i] = (ushort)(s0 - Sm[i]);
 543                 }
 544             }
 545             dst += dststep;
 546         }
 547     }
 548
 549     double scale;
 550     int sumCount;
 551     int divDelta;
 552     int divScale;
 553     std::vector<ushort> sum;
 554 };
 555
 556
 557 template<>
 558 struct ColumnSum<int, short> :
 559         public BaseColumnFilter
 560 {
 561     ColumnSum( int _ksize, int _anchor, double _scale ) :
 562         BaseColumnFilter()
 563     {
 564         ksize = _ksize;
 565         anchor = _anchor;
 566         scale = _scale;
 567         sumCount = 0;
 568     }
 569
 570     virtual void reset() { sumCount = 0; }
 571
 572     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 573     {
 574         int i;
 575         int* SUM;
 576         bool haveScale = scale != 1;
 577         double _scale = scale;
 578
 579 #if CV_SIMD128
 580         bool haveSIMD128 = hasSIMD128();
 581 #endif
 582
 583         if( width != (int)sum.size() )
 584         {
 585             sum.resize(width);
 586             sumCount = 0;
 587         }
 588
 589         SUM = &sum[0];
 590         if( sumCount == 0 )
 591         {
 592             memset((void*)SUM, 0, width*sizeof(int));
 593             for( ; sumCount < ksize - 1; sumCount++, src++ )
 594             {
 595                 const int* Sp = (const int*)src[0];
 596                 i = 0;
 597 #if CV_SIMD128
 598                 if( haveSIMD128 )
 599                 {
 600                     for( ; i <= width - 4; i+=4 )
 601                     {
 602                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 603                     }
 604                 }
 605                 #endif
 606                 for( ; i < width; i++ )
 607                     SUM[i] += Sp[i];
 608             }
 609         }
 610         else
 611         {
 612             CV_Assert( sumCount == ksize-1 );
 613             src += ksize-1;
 614         }
 615
 616         for( ; count--; src++ )
 617         {
 618             const int* Sp = (const int*)src[0];
 619             const int* Sm = (const int*)src[1-ksize];
 620             short* D = (short*)dst;
 621             if( haveScale )
 622             {
 623                 i = 0;
 624 #if CV_SIMD128
 625                 if( haveSIMD128 )
 626                 {
 627                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 628                     for( ; i <= width-8; i+=8 )
 629                     {
 630                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 631                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 632
 633                         v_int32x4 v_s0d =  v_round(v_cvt_f32(v_s0) * v_scale);
 634                         v_int32x4 v_s01d = v_round(v_cvt_f32(v_s01) * v_scale);
 635                         v_store(D + i, v_pack(v_s0d, v_s01d));
 636
 637                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 638                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 639                     }
 640                 }
 641 #endif
 642                 for( ; i < width; i++ )
 643                 {
 644                     int s0 = SUM[i] + Sp[i];
 645                     D[i] = saturate_cast<short>(s0*_scale);
 646                     SUM[i] = s0 - Sm[i];
 647                 }
 648             }
 649             else
 650             {
 651                 i = 0;
 652 #if CV_SIMD128
 653                 if( haveSIMD128 )
 654                 {
 655                     for( ; i <= width-8; i+=8 )
 656                     {
 657                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 658                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 659
 660                         v_store(D + i, v_pack(v_s0, v_s01));
 661
 662                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 663                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 664                     }
 665                 }
 666 #endif
 667
 668                 for( ; i < width; i++ )
 669                 {
 670                     int s0 = SUM[i] + Sp[i];
 671                     D[i] = saturate_cast<short>(s0);
 672                     SUM[i] = s0 - Sm[i];
 673                 }
 674             }
 675             dst += dststep;
 676         }
 677     }
 678
 679     double scale;
 680     int sumCount;
 681     std::vector<int> sum;
 682 };
 683
 684
 685 template<>
 686 struct ColumnSum<int, ushort> :
 687         public BaseColumnFilter
 688 {
 689     ColumnSum( int _ksize, int _anchor, double _scale ) :
 690         BaseColumnFilter()
 691     {
 692         ksize = _ksize;
 693         anchor = _anchor;
 694         scale = _scale;
 695         sumCount = 0;
 696     }
 697
 698     virtual void reset() { sumCount = 0; }
 699
 700     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 701     {
 702         int* SUM;
 703         bool haveScale = scale != 1;
 704         double _scale = scale;
 705
 706 #if CV_SIMD128
 707         bool haveSIMD128 = hasSIMD128();
 708 #endif
 709
 710         if( width != (int)sum.size() )
 711         {
 712             sum.resize(width);
 713             sumCount = 0;
 714         }
 715
 716         SUM = &sum[0];
 717         if( sumCount == 0 )
 718         {
 719             memset((void*)SUM, 0, width*sizeof(int));
 720             for( ; sumCount < ksize - 1; sumCount++, src++ )
 721             {
 722                 const int* Sp = (const int*)src[0];
 723                 int i = 0;
 724 #if CV_SIMD128
 725                 if( haveSIMD128 )
 726                 {
 727                     for (; i <= width - 4; i += 4)
 728                     {
 729                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 730                     }
 731                 }
 732 #endif
 733                 for( ; i < width; i++ )
 734                     SUM[i] += Sp[i];
 735             }
 736         }
 737         else
 738         {
 739             CV_Assert( sumCount == ksize-1 );
 740             src += ksize-1;
 741         }
 742
 743         for( ; count--; src++ )
 744         {
 745             const int* Sp = (const int*)src[0];
 746             const int* Sm = (const int*)src[1-ksize];
 747             ushort* D = (ushort*)dst;
 748             if( haveScale )
 749             {
 750                 int i = 0;
 751 #if CV_SIMD128
 752                 if( haveSIMD128 )
 753                 {
 754                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 755                     for( ; i <= width-8; i+=8 )
 756                     {
 757                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 758                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 759
 760                         v_uint32x4 v_s0d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s0) * v_scale));
 761                         v_uint32x4 v_s01d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s01) * v_scale));
 762                         v_store(D + i, v_pack(v_s0d, v_s01d));
 763
 764                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 765                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 766                     }
 767                 }
 768 #endif
 769                 for( ; i < width; i++ )
 770                 {
 771                     int s0 = SUM[i] + Sp[i];
 772                     D[i] = saturate_cast<ushort>(s0*_scale);
 773                     SUM[i] = s0 - Sm[i];
 774                 }
 775             }
 776             else
 777             {
 778                 int i = 0;
 779 #if CV_SIMD128
 780                 if( haveSIMD128 )
 781                 {
 782                     for( ; i <= width-8; i+=8 )
 783                     {
 784                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 785                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 786
 787                         v_store(D + i, v_pack(v_reinterpret_as_u32(v_s0), v_reinterpret_as_u32(v_s01)));
 788
 789                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 790                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 791                     }
 792                 }
 793 #endif
 794                 for( ; i < width; i++ )
 795                 {
 796                     int s0 = SUM[i] + Sp[i];
 797                     D[i] = saturate_cast<ushort>(s0);
 798                     SUM[i] = s0 - Sm[i];
 799                 }
 800             }
 801             dst += dststep;
 802         }
 803     }
 804
 805     double scale;
 806     int sumCount;
 807     std::vector<int> sum;
 808 };
 809
 810 template<>
 811 struct ColumnSum<int, int> :
 812         public BaseColumnFilter
 813 {
 814     ColumnSum( int _ksize, int _anchor, double _scale ) :
 815         BaseColumnFilter()
 816     {
 817         ksize = _ksize;
 818         anchor = _anchor;
 819         scale = _scale;
 820         sumCount = 0;
 821     }
 822
 823     virtual void reset() { sumCount = 0; }
 824
 825     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 826     {
 827         int* SUM;
 828         bool haveScale = scale != 1;
 829         double _scale = scale;
 830
 831 #if CV_SIMD128
 832         bool haveSIMD128 = hasSIMD128();
 833 #endif
 834
 835         if( width != (int)sum.size() )
 836         {
 837             sum.resize(width);
 838             sumCount = 0;
 839         }
 840
 841         SUM = &sum[0];
 842         if( sumCount == 0 )
 843         {
 844             memset((void*)SUM, 0, width*sizeof(int));
 845             for( ; sumCount < ksize - 1; sumCount++, src++ )
 846             {
 847                 const int* Sp = (const int*)src[0];
 848                 int i = 0;
 849 #if CV_SIMD128
 850                 if( haveSIMD128 )
 851                 {
 852                     for( ; i <= width - 4; i+=4 )
 853                     {
 854                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 855                     }
 856                 }
 857 #endif
 858                 for( ; i < width; i++ )
 859                     SUM[i] += Sp[i];
 860             }
 861         }
 862         else
 863         {
 864             CV_Assert( sumCount == ksize-1 );
 865             src += ksize-1;
 866         }
 867
 868         for( ; count--; src++ )
 869         {
 870             const int* Sp = (const int*)src[0];
 871             const int* Sm = (const int*)src[1-ksize];
 872             int* D = (int*)dst;
 873             if( haveScale )
 874             {
 875                 int i = 0;
 876 #if CV_SIMD128
 877                 if( haveSIMD128 )
 878                 {
 879                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 880                     for( ; i <= width-4; i+=4 )
 881                     {
 882                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 883                         v_int32x4 v_s0d = v_round(v_cvt_f32(v_s0) * v_scale);
 884
 885                         v_store(D + i, v_s0d);
 886                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 887                     }
 888                 }
 889 #endif
 890                 for( ; i < width; i++ )
 891                 {
 892                     int s0 = SUM[i] + Sp[i];
 893                     D[i] = saturate_cast<int>(s0*_scale);
 894                     SUM[i] = s0 - Sm[i];
 895                 }
 896             }
 897             else
 898             {
 899                 int i = 0;
 900 #if CV_SIMD128
 901                 if( haveSIMD128 )
 902                 {
 903                     for( ; i <= width-4; i+=4 )
 904                     {
 905                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 906
 907                         v_store(D + i, v_s0);
 908                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 909                     }
 910                 }
 911 #endif
 912                 for( ; i < width; i++ )
 913                 {
 914                     int s0 = SUM[i] + Sp[i];
 915                     D[i] = s0;
 916                     SUM[i] = s0 - Sm[i];
 917                 }
 918             }
 919             dst += dststep;
 920         }
 921     }
 922
 923     double scale;
 924     int sumCount;
 925     std::vector<int> sum;
 926 };
 927
 928
 929 template<>
 930 struct ColumnSum<int, float> :
 931         public BaseColumnFilter
 932 {
 933     ColumnSum( int _ksize, int _anchor, double _scale ) :
 934         BaseColumnFilter()
 935     {
 936         ksize = _ksize;
 937         anchor = _anchor;
 938         scale = _scale;
 939         sumCount = 0;
 940     }
 941
 942     virtual void reset() { sumCount = 0; }
 943
 944     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 945     {
 946         int* SUM;
 947         bool haveScale = scale != 1;
 948         double _scale = scale;
 949
 950 #if CV_SIMD128
 951         bool haveSIMD128 = hasSIMD128();
 952 #endif
 953
 954         if( width != (int)sum.size() )
 955         {
 956             sum.resize(width);
 957             sumCount = 0;
 958         }
 959
 960         SUM = &sum[0];
 961         if( sumCount == 0 )
 962         {
 963             memset((void*)SUM, 0, width*sizeof(int));
 964             for( ; sumCount < ksize - 1; sumCount++, src++ )
 965             {
 966                 const int* Sp = (const int*)src[0];
 967                 int i = 0;
 968 #if CV_SIMD128
 969                 if( haveSIMD128 )
 970                 {
 971                     for( ; i <= width - 4; i+=4 )
 972                     {
 973                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 974                     }
 975                 }
 976 #endif
 977
 978                 for( ; i < width; i++ )
 979                     SUM[i] += Sp[i];
 980             }
 981         }
 982         else
 983         {
 984             CV_Assert( sumCount == ksize-1 );
 985             src += ksize-1;
 986         }
 987
 988         for( ; count--; src++ )
 989         {
 990             const int * Sp = (const int*)src[0];
 991             const int * Sm = (const int*)src[1-ksize];
 992             float* D = (float*)dst;
 993             if( haveScale )
 994             {
 995                 int i = 0;
 996
 997 #if CV_SIMD128
 998                 if( haveSIMD128 )
 999                 {
1000                     v_float32x4 v_scale = v_setall_f32((float)_scale);
1001                     for (; i <= width - 8; i += 8)
1002                     {
1003                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
1004                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
1005
1006                         v_store(D + i, v_cvt_f32(v_s0) * v_scale);
1007                         v_store(D + i + 4, v_cvt_f32(v_s01) * v_scale);
1008
1009                         v_store(SUM + i, v_s0 - v_load(Sm + i));
1010                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
1011                     }
1012                 }
1013 #endif
1014                 for( ; i < width; i++ )
1015                 {
1016                     int s0 = SUM[i] + Sp[i];
1017                     D[i] = (float)(s0*_scale);
1018                     SUM[i] = s0 - Sm[i];
1019                 }
1020             }
1021             else
1022             {
1023                 int i = 0;
1024
1025 #if CV_SIMD128
1026                 if( haveSIMD128 )
1027                 {
1028                     for( ; i <= width-8; i+=8 )
1029                     {
1030                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
1031                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
1032
1033                         v_store(D + i, v_cvt_f32(v_s0));
1034                         v_store(D + i + 4, v_cvt_f32(v_s01));
1035
1036                         v_store(SUM + i, v_s0 - v_load(Sm + i));
1037                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
1038                     }
1039                 }
1040 #endif
1041                 for( ; i < width; i++ )
1042                 {
1043                     int s0 = SUM[i] + Sp[i];
1044                     D[i] = (float)(s0);
1045                     SUM[i] = s0 - Sm[i];
1046                 }
1047             }
1048             dst += dststep;
1049         }
1050     }
1051
1052     double scale;
1053     int sumCount;
1054     std::vector<int> sum;
1055 };
1056
1057 #ifdef HAVE_OPENCL
1058
1059 static bool ocl_boxFilter3x3_8UC1( InputArray _src, OutputArray _dst, int ddepth,
1060                                    Size ksize, Point anchor, int borderType, bool normalize )
1061 {
1062     const ocl::Device & dev = ocl::Device::getDefault();
1063     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
1064
1065     if (ddepth < 0)
1066         ddepth = sdepth;
1067
1068     if (anchor.x < 0)
1069         anchor.x = ksize.width / 2;
1070     if (anchor.y < 0)
1071         anchor.y = ksize.height / 2;
1072
1073     if ( !(dev.isIntel() && (type == CV_8UC1) &&
1074          (_src.offset() == 0) && (_src.step() % 4 == 0) &&
1075          (_src.cols() % 16 == 0) && (_src.rows() % 2 == 0) &&
1076          (anchor.x == 1) && (anchor.y == 1) &&
1077          (ksize.width == 3) && (ksize.height == 3)) )
1078         return false;
1079
1080     float alpha = 1.0f / (ksize.height * ksize.width);
1081     Size size = _src.size();
1082     size_t globalsize[2] = { 0, 0 };
1083     size_t localsize[2] = { 0, 0 };
1084     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
1085
1086     globalsize[0] = size.width / 16;
1087     globalsize[1] = size.height / 2;
1088
1089     char build_opts[1024];
1090     sprintf(build_opts, "-D %s %s", borderMap[borderType], normalize ? "-D NORMALIZE" : "");
1091
1092     ocl::Kernel kernel("boxFilter3x3_8UC1_cols16_rows2", cv::ocl::imgproc::boxFilter3x3_oclsrc, build_opts);
1093     if (kernel.empty())
1094         return false;
1095
1096     UMat src = _src.getUMat();
1097     _dst.create(size, CV_MAKETYPE(ddepth, cn));
1098     if (!(_dst.offset() == 0 && _dst.step() % 4 == 0))
1099         return false;
1100     UMat dst = _dst.getUMat();
1101
1102     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
1103     idxArg = kernel.set(idxArg, (int)src.step);
1104     idxArg = kernel.set(idxArg, ocl::KernelArg::PtrWriteOnly(dst));
1105     idxArg = kernel.set(idxArg, (int)dst.step);
1106     idxArg = kernel.set(idxArg, (int)dst.rows);
1107     idxArg = kernel.set(idxArg, (int)dst.cols);
1108     if (normalize)
1109         idxArg = kernel.set(idxArg, (float)alpha);
1110
1111     return kernel.run(2, globalsize, (localsize[0] == 0) ? NULL : localsize, false);
1112 }
1113
1114 #define DIVUP(total, grain) ((total + grain - 1) / (grain))
1115 #define ROUNDUP(sz, n)      ((sz) + (n) - 1 - (((sz) + (n) - 1) % (n)))
1116
1117 static bool ocl_boxFilter( InputArray _src, OutputArray _dst, int ddepth,
1118                            Size ksize, Point anchor, int borderType, bool normalize, bool sqr = false )
1119 {
1120     const ocl::Device & dev = ocl::Device::getDefault();
1121     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type), esz = CV_ELEM_SIZE(type);
1122     bool doubleSupport = dev.doubleFPConfig() > 0;
1123
1124     if (ddepth < 0)
1125         ddepth = sdepth;
1126
1127     if (cn > 4 || (!doubleSupport && (sdepth == CV_64F || ddepth == CV_64F)) ||
1128         _src.offset() % esz != 0 || _src.step() % esz != 0)
1129         return false;
1130
1131     if (anchor.x < 0)
1132         anchor.x = ksize.width / 2;
1133     if (anchor.y < 0)
1134         anchor.y = ksize.height / 2;
1135
1136     int computeUnits = ocl::Device::getDefault().maxComputeUnits();
1137     float alpha = 1.0f / (ksize.height * ksize.width);
1138     Size size = _src.size(), wholeSize;
1139     bool isolated = (borderType & BORDER_ISOLATED) != 0;
1140     borderType &= ~BORDER_ISOLATED;
1141     int wdepth = std::max(CV_32F, std::max(ddepth, sdepth)),
1142         wtype = CV_MAKE_TYPE(wdepth, cn), dtype = CV_MAKE_TYPE(ddepth, cn);
1143
1144     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
1145     size_t globalsize[2] = { (size_t)size.width, (size_t)size.height };
1146     size_t localsize_general[2] = { 0, 1 }, * localsize = NULL;
1147
1148     UMat src = _src.getUMat();
1149     if (!isolated)
1150     {
1151         Point ofs;
1152         src.locateROI(wholeSize, ofs);
1153     }
1154
1155     int h = isolated ? size.height : wholeSize.height;
1156     int w = isolated ? size.width : wholeSize.width;
1157
1158     size_t maxWorkItemSizes[32];
1159     ocl::Device::getDefault().maxWorkItemSizes(maxWorkItemSizes);
1160     int tryWorkItems = (int)maxWorkItemSizes[0];
1161
1162     ocl::Kernel kernel;
1163
1164     if (dev.isIntel() && !(dev.type() & ocl::Device::TYPE_CPU) &&
1165         ((ksize.width < 5 && ksize.height < 5 && esz <= 4) ||
1166          (ksize.width == 5 && ksize.height == 5 && cn == 1)))
1167     {
1168         if (w < ksize.width || h < ksize.height)
1169             return false;
1170
1171         // Figure out what vector size to use for loading the pixels.
1172         int pxLoadNumPixels = cn != 1 || size.width % 4 ? 1 : 4;
1173         int pxLoadVecSize = cn * pxLoadNumPixels;
1174
1175         // Figure out how many pixels per work item to compute in X and Y
1176         // directions.  Too many and we run out of registers.
1177         int pxPerWorkItemX = 1, pxPerWorkItemY = 1;
1178         if (cn <= 2 && ksize.width <= 4 && ksize.height <= 4)
1179         {
1180             pxPerWorkItemX = size.width % 8 ? size.width % 4 ? size.width % 2 ? 1 : 2 : 4 : 8;
1181             pxPerWorkItemY = size.height % 2 ? 1 : 2;
1182         }
1183         else if (cn < 4 || (ksize.width <= 4 && ksize.height <= 4))
1184         {
1185             pxPerWorkItemX = size.width % 2 ? 1 : 2;
1186             pxPerWorkItemY = size.height % 2 ? 1 : 2;
1187         }
1188         globalsize[0] = size.width / pxPerWorkItemX;
1189         globalsize[1] = size.height / pxPerWorkItemY;
1190
1191         // Need some padding in the private array for pixels
1192         int privDataWidth = ROUNDUP(pxPerWorkItemX + ksize.width - 1, pxLoadNumPixels);
1193
1194         // Make the global size a nice round number so the runtime can pick
1195         // from reasonable choices for the workgroup size
1196         const int wgRound = 256;
1197         globalsize[0] = ROUNDUP(globalsize[0], wgRound);
1198
1199         char build_options[1024], cvt[2][40];
1200         sprintf(build_options, "-D cn=%d "
1201                 "-D ANCHOR_X=%d -D ANCHOR_Y=%d -D KERNEL_SIZE_X=%d -D KERNEL_SIZE_Y=%d "
1202                 "-D PX_LOAD_VEC_SIZE=%d -D PX_LOAD_NUM_PX=%d "
1203                 "-D PX_PER_WI_X=%d -D PX_PER_WI_Y=%d -D PRIV_DATA_WIDTH=%d -D %s -D %s "
1204                 "-D PX_LOAD_X_ITERATIONS=%d -D PX_LOAD_Y_ITERATIONS=%d "
1205                 "-D srcT=%s -D srcT1=%s -D dstT=%s -D dstT1=%s -D WT=%s -D WT1=%s "
1206                 "-D convertToWT=%s -D convertToDstT=%s%s%s -D PX_LOAD_FLOAT_VEC_CONV=convert_%s -D OP_BOX_FILTER",
1207                 cn, anchor.x, anchor.y, ksize.width, ksize.height,
1208                 pxLoadVecSize, pxLoadNumPixels,
1209                 pxPerWorkItemX, pxPerWorkItemY, privDataWidth, borderMap[borderType],
1210                 isolated ? "BORDER_ISOLATED" : "NO_BORDER_ISOLATED",
1211                 privDataWidth / pxLoadNumPixels, pxPerWorkItemY + ksize.height - 1,
1212                 ocl::typeToStr(type), ocl::typeToStr(sdepth), ocl::typeToStr(dtype),
1213                 ocl::typeToStr(ddepth), ocl::typeToStr(wtype), ocl::typeToStr(wdepth),
1214                 ocl::convertTypeStr(sdepth, wdepth, cn, cvt[0]),
1215                 ocl::convertTypeStr(wdepth, ddepth, cn, cvt[1]),
1216                 normalize ? " -D NORMALIZE" : "", sqr ? " -D SQR" : "",
1217                 ocl::typeToStr(CV_MAKE_TYPE(wdepth, pxLoadVecSize)) //PX_LOAD_FLOAT_VEC_CONV
1218                 );
1219
1220
1221         if (!kernel.create("filterSmall", cv::ocl::imgproc::filterSmall_oclsrc, build_options))
1222             return false;
1223     }
1224     else
1225     {
1226         localsize = localsize_general;
1227         for ( ; ; )
1228         {
1229             int BLOCK_SIZE_X = tryWorkItems, BLOCK_SIZE_Y = std::min(ksize.height * 10, size.height);
1230
1231             while (BLOCK_SIZE_X > 32 && BLOCK_SIZE_X >= ksize.width * 2 && BLOCK_SIZE_X > size.width * 2)
1232                 BLOCK_SIZE_X /= 2;
1233             while (BLOCK_SIZE_Y < BLOCK_SIZE_X / 8 && BLOCK_SIZE_Y * computeUnits * 32 < size.height)
1234                 BLOCK_SIZE_Y *= 2;
1235
1236             if (ksize.width > BLOCK_SIZE_X || w < ksize.width || h < ksize.height)
1237                 return false;
1238
1239             char cvt[2][50];
1240             String opts = format("-D LOCAL_SIZE_X=%d -D BLOCK_SIZE_Y=%d -D ST=%s -D DT=%s -D WT=%s -D convertToDT=%s -D convertToWT=%s"
1241                                  " -D ANCHOR_X=%d -D ANCHOR_Y=%d -D KERNEL_SIZE_X=%d -D KERNEL_SIZE_Y=%d -D %s%s%s%s%s"
1242                                  " -D ST1=%s -D DT1=%s -D cn=%d",
1243                                  BLOCK_SIZE_X, BLOCK_SIZE_Y, ocl::typeToStr(type), ocl::typeToStr(CV_MAKE_TYPE(ddepth, cn)),
1244                                  ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)),
1245                                  ocl::convertTypeStr(wdepth, ddepth, cn, cvt[0]),
1246                                  ocl::convertTypeStr(sdepth, wdepth, cn, cvt[1]),
1247                                  anchor.x, anchor.y, ksize.width, ksize.height, borderMap[borderType],
1248                                  isolated ? " -D BORDER_ISOLATED" : "", doubleSupport ? " -D DOUBLE_SUPPORT" : "",
1249                                  normalize ? " -D NORMALIZE" : "", sqr ? " -D SQR" : "",
1250                                  ocl::typeToStr(sdepth), ocl::typeToStr(ddepth), cn);
1251
1252             localsize[0] = BLOCK_SIZE_X;
1253             globalsize[0] = DIVUP(size.width, BLOCK_SIZE_X - (ksize.width - 1)) * BLOCK_SIZE_X;
1254             globalsize[1] = DIVUP(size.height, BLOCK_SIZE_Y);
1255
1256             kernel.create("boxFilter", cv::ocl::imgproc::boxFilter_oclsrc, opts);
1257             if (kernel.empty())
1258                 return false;
1259
1260             size_t kernelWorkGroupSize = kernel.workGroupSize();
1261             if (localsize[0] <= kernelWorkGroupSize)
1262                 break;
1263             if (BLOCK_SIZE_X < (int)kernelWorkGroupSize)
1264                 return false;
1265
1266             tryWorkItems = (int)kernelWorkGroupSize;
1267         }
1268     }
1269
1270     _dst.create(size, CV_MAKETYPE(ddepth, cn));
1271     UMat dst = _dst.getUMat();
1272
1273     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
1274     idxArg = kernel.set(idxArg, (int)src.step);
1275     int srcOffsetX = (int)((src.offset % src.step) / src.elemSize());
1276     int srcOffsetY = (int)(src.offset / src.step);
1277     int srcEndX = isolated ? srcOffsetX + size.width : wholeSize.width;
1278     int srcEndY = isolated ? srcOffsetY + size.height : wholeSize.height;
1279     idxArg = kernel.set(idxArg, srcOffsetX);
1280     idxArg = kernel.set(idxArg, srcOffsetY);
1281     idxArg = kernel.set(idxArg, srcEndX);
1282     idxArg = kernel.set(idxArg, srcEndY);
1283     idxArg = kernel.set(idxArg, ocl::KernelArg::WriteOnly(dst));
1284     if (normalize)
1285         idxArg = kernel.set(idxArg, (float)alpha);
1286
1287     return kernel.run(2, globalsize, localsize, false);
1288 }
1289
1290 #undef ROUNDUP
1291
1292 #endif
1293
1294 }
1295
1296
1297 cv::Ptr<cv::BaseRowFilter> cv::getRowSumFilter(int srcType, int sumType, int ksize, int anchor)
1298 {
1299     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
1300     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
1301
1302     if( anchor < 0 )
1303         anchor = ksize/2;
1304
1305     if( sdepth == CV_8U && ddepth == CV_32S )
1306         return makePtr<RowSum<uchar, int> >(ksize, anchor);
1307     if( sdepth == CV_8U && ddepth == CV_16U )
1308         return makePtr<RowSum<uchar, ushort> >(ksize, anchor);
1309     if( sdepth == CV_8U && ddepth == CV_64F )
1310         return makePtr<RowSum<uchar, double> >(ksize, anchor);
1311     if( sdepth == CV_16U && ddepth == CV_32S )
1312         return makePtr<RowSum<ushort, int> >(ksize, anchor);
1313     if( sdepth == CV_16U && ddepth == CV_64F )
1314         return makePtr<RowSum<ushort, double> >(ksize, anchor);
1315     if( sdepth == CV_16S && ddepth == CV_32S )
1316         return makePtr<RowSum<short, int> >(ksize, anchor);
1317     if( sdepth == CV_32S && ddepth == CV_32S )
1318         return makePtr<RowSum<int, int> >(ksize, anchor);
1319     if( sdepth == CV_16S && ddepth == CV_64F )
1320         return makePtr<RowSum<short, double> >(ksize, anchor);
1321     if( sdepth == CV_32F && ddepth == CV_64F )
1322         return makePtr<RowSum<float, double> >(ksize, anchor);
1323     if( sdepth == CV_64F && ddepth == CV_64F )
1324         return makePtr<RowSum<double, double> >(ksize, anchor);
1325
1326     CV_Error_( CV_StsNotImplemented,
1327         ("Unsupported combination of source format (=%d), and buffer format (=%d)",
1328         srcType, sumType));
1329
1330     return Ptr<BaseRowFilter>();
1331 }
1332
1333
1334 cv::Ptr<cv::BaseColumnFilter> cv::getColumnSumFilter(int sumType, int dstType, int ksize,
1335                                                      int anchor, double scale)
1336 {
1337     int sdepth = CV_MAT_DEPTH(sumType), ddepth = CV_MAT_DEPTH(dstType);
1338     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(dstType) );
1339
1340     if( anchor < 0 )
1341         anchor = ksize/2;
1342
1343     if( ddepth == CV_8U && sdepth == CV_32S )
1344         return makePtr<ColumnSum<int, uchar> >(ksize, anchor, scale);
1345     if( ddepth == CV_8U && sdepth == CV_16U )
1346         return makePtr<ColumnSum<ushort, uchar> >(ksize, anchor, scale);
1347     if( ddepth == CV_8U && sdepth == CV_64F )
1348         return makePtr<ColumnSum<double, uchar> >(ksize, anchor, scale);
1349     if( ddepth == CV_16U && sdepth == CV_32S )
1350         return makePtr<ColumnSum<int, ushort> >(ksize, anchor, scale);
1351     if( ddepth == CV_16U && sdepth == CV_64F )
1352         return makePtr<ColumnSum<double, ushort> >(ksize, anchor, scale);
1353     if( ddepth == CV_16S && sdepth == CV_32S )
1354         return makePtr<ColumnSum<int, short> >(ksize, anchor, scale);
1355     if( ddepth == CV_16S && sdepth == CV_64F )
1356         return makePtr<ColumnSum<double, short> >(ksize, anchor, scale);
1357     if( ddepth == CV_32S && sdepth == CV_32S )
1358         return makePtr<ColumnSum<int, int> >(ksize, anchor, scale);
1359     if( ddepth == CV_32F && sdepth == CV_32S )
1360         return makePtr<ColumnSum<int, float> >(ksize, anchor, scale);
1361     if( ddepth == CV_32F && sdepth == CV_64F )
1362         return makePtr<ColumnSum<double, float> >(ksize, anchor, scale);
1363     if( ddepth == CV_64F && sdepth == CV_32S )
1364         return makePtr<ColumnSum<int, double> >(ksize, anchor, scale);
1365     if( ddepth == CV_64F && sdepth == CV_64F )
1366         return makePtr<ColumnSum<double, double> >(ksize, anchor, scale);
1367
1368     CV_Error_( CV_StsNotImplemented,
1369         ("Unsupported combination of sum format (=%d), and destination format (=%d)",
1370         sumType, dstType));
1371
1372     return Ptr<BaseColumnFilter>();
1373 }
1374
1375
1376 cv::Ptr<cv::FilterEngine> cv::createBoxFilter( int srcType, int dstType, Size ksize,
1377                     Point anchor, bool normalize, int borderType )
1378 {
1379     int sdepth = CV_MAT_DEPTH(srcType);
1380     int cn = CV_MAT_CN(srcType), sumType = CV_64F;
1381     if( sdepth == CV_8U && CV_MAT_DEPTH(dstType) == CV_8U &&
1382         ksize.width*ksize.height <= 256 )
1383         sumType = CV_16U;
1384     else if( sdepth <= CV_32S && (!normalize ||
1385         ksize.width*ksize.height <= (sdepth == CV_8U ? (1<<23) :
1386             sdepth == CV_16U ? (1 << 15) : (1 << 16))) )
1387         sumType = CV_32S;
1388     sumType = CV_MAKETYPE( sumType, cn );
1389
1390     Ptr<BaseRowFilter> rowFilter = getRowSumFilter(srcType, sumType, ksize.width, anchor.x );
1391     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
1392         dstType, ksize.height, anchor.y, normalize ? 1./(ksize.width*ksize.height) : 1);
1393
1394     return makePtr<FilterEngine>(Ptr<BaseFilter>(), rowFilter, columnFilter,
1395            srcType, dstType, sumType, borderType );
1396 }
1397
1398 #ifdef HAVE_OPENVX
1399 namespace cv
1400 {
1401     namespace ovx {
1402         template <> inline bool skipSmallImages<VX_KERNEL_BOX_3x3>(int w, int h) { return w*h < 640 * 480; }
1403     }
1404     static bool openvx_boxfilter(InputArray _src, OutputArray _dst, int ddepth,
1405                                  Size ksize, Point anchor,
1406                                  bool normalize, int borderType)
1407     {
1408         if (ddepth < 0)
1409             ddepth = CV_8UC1;
1410         if (_src.type() != CV_8UC1 || ddepth != CV_8U || !normalize ||
1411             _src.cols() < 3 || _src.rows() < 3 ||
1412             ksize.width != 3 || ksize.height != 3 ||
1413             (anchor.x >= 0 && anchor.x != 1) ||
1414             (anchor.y >= 0 && anchor.y != 1) ||
1415             ovx::skipSmallImages<VX_KERNEL_BOX_3x3>(_src.cols(), _src.rows()))
1416             return false;
1417
1418         Mat src = _src.getMat();
1419
1420         if ((borderType & BORDER_ISOLATED) == 0 && src.isSubmatrix())
1421             return false; //Process isolated borders only
1422         vx_enum border;
1423         switch (borderType & ~BORDER_ISOLATED)
1424         {
1425         case BORDER_CONSTANT:
1426             border = VX_BORDER_CONSTANT;
1427             break;
1428         case BORDER_REPLICATE:
1429             border = VX_BORDER_REPLICATE;
1430             break;
1431         default:
1432             return false;
1433         }
1434
1435         _dst.create(src.size(), CV_8UC1);
1436         Mat dst = _dst.getMat();
1437
1438         try
1439         {
1440             ivx::Context ctx = ovx::getOpenVXContext();
1441
1442             Mat a;
1443             if (dst.data != src.data)
1444                 a = src;
1445             else
1446                 src.copyTo(a);
1447
1448             ivx::Image
1449                 ia = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
1450                                                   ivx::Image::createAddressing(a.cols, a.rows, 1, (vx_int32)(a.step)), a.data),
1451                 ib = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
1452                                                   ivx::Image::createAddressing(dst.cols, dst.rows, 1, (vx_int32)(dst.step)), dst.data);
1453
1454             //ATTENTION: VX_CONTEXT_IMMEDIATE_BORDER attribute change could lead to strange issues in multi-threaded environments
1455             //since OpenVX standart says nothing about thread-safety for now
1456             ivx::border_t prevBorder = ctx.immediateBorder();
1457             ctx.setImmediateBorder(border, (vx_uint8)(0));
1458             ivx::IVX_CHECK_STATUS(vxuBox3x3(ctx, ia, ib));
1459             ctx.setImmediateBorder(prevBorder);
1460         }
1461         catch (ivx::RuntimeError & e)
1462         {
1463             VX_DbgThrow(e.what());
1464         }
1465         catch (ivx::WrapperError & e)
1466         {
1467             VX_DbgThrow(e.what());
1468         }
1469
1470         return true;
1471     }
1472 }
1473 #endif
1474
1475 #if defined(HAVE_IPP)
1476 namespace cv
1477 {
1478 static bool ipp_boxfilter(Mat &src, Mat &dst, Size ksize, Point anchor, bool normalize, int borderType)
1479 {
1480 #ifdef HAVE_IPP_IW
1481     CV_INSTRUMENT_REGION_IPP()
1482
1483 #if IPP_VERSION_X100 < 201801
1484     // Problem with SSE42 optimization for 16s and some 8u modes
1485     if(ipp::getIppTopFeatures() == ippCPUID_SSE42 && (((src.depth() == CV_16S || src.depth() == CV_16U) && (src.channels() == 3 || src.channels() == 4)) || (src.depth() == CV_8U && src.channels() == 3 && (ksize.width > 5 || ksize.height > 5))))
1486         return false;
1487
1488     // Other optimizations has some degradations too
1489     if((((src.depth() == CV_16S || src.depth() == CV_16U) && (src.channels() == 4)) || (src.depth() == CV_8U && src.channels() == 1 && (ksize.width > 5 || ksize.height > 5))))
1490         return false;
1491 #endif
1492
1493     if(!normalize)
1494         return false;
1495
1496     if(!ippiCheckAnchor(anchor, ksize))
1497         return false;
1498
1499     try
1500     {
1501         ::ipp::IwiImage       iwSrc      = ippiGetImage(src);
1502         ::ipp::IwiImage       iwDst      = ippiGetImage(dst);
1503         ::ipp::IwiSize        iwKSize    = ippiGetSize(ksize);
1504         ::ipp::IwiBorderSize  borderSize(iwKSize);
1505         ::ipp::IwiBorderType  ippBorder(ippiGetBorder(iwSrc, borderType, borderSize));
1506         if(!ippBorder)
1507             return false;
1508
1509         CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBox, iwSrc, iwDst, iwKSize, ::ipp::IwDefault(), ippBorder);
1510     }
1511     catch (::ipp::IwException)
1512     {
1513         return false;
1514     }
1515
1516     return true;
1517 #else
1518     CV_UNUSED(src); CV_UNUSED(dst); CV_UNUSED(ksize); CV_UNUSED(anchor); CV_UNUSED(normalize); CV_UNUSED(borderType);
1519     return false;
1520 #endif
1521 }
1522 }
1523 #endif
1524
1525
1526 void cv::boxFilter( InputArray _src, OutputArray _dst, int ddepth,
1527                 Size ksize, Point anchor,
1528                 bool normalize, int borderType )
1529 {
1530     CV_INSTRUMENT_REGION()
1531
1532     CV_OCL_RUN(_dst.isUMat() &&
1533                (borderType == BORDER_REPLICATE || borderType == BORDER_CONSTANT ||
1534                 borderType == BORDER_REFLECT || borderType == BORDER_REFLECT_101),
1535                ocl_boxFilter3x3_8UC1(_src, _dst, ddepth, ksize, anchor, borderType, normalize))
1536
1537     CV_OCL_RUN(_dst.isUMat(), ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize))
1538
1539     CV_OVX_RUN(true,
1540                openvx_boxfilter(_src, _dst, ddepth, ksize, anchor, normalize, borderType))
1541
1542     Mat src = _src.getMat();
1543     int stype = src.type(), sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype);
1544     if( ddepth < 0 )
1545         ddepth = sdepth;
1546     _dst.create( src.size(), CV_MAKETYPE(ddepth, cn) );
1547     Mat dst = _dst.getMat();
1548     if( borderType != BORDER_CONSTANT && normalize && (borderType & BORDER_ISOLATED) != 0 )
1549     {
1550         if( src.rows == 1 )
1551             ksize.height = 1;
1552         if( src.cols == 1 )
1553             ksize.width = 1;
1554     }
1555 #ifdef HAVE_TEGRA_OPTIMIZATION
1556     if ( tegra::useTegra() && tegra::box(src, dst, ksize, anchor, normalize, borderType) )
1557         return;
1558 #endif
1559
1560     CV_IPP_RUN_FAST(ipp_boxfilter(src, dst, ksize, anchor, normalize, borderType));
1561
1562     Point ofs;
1563     Size wsz(src.cols, src.rows);
1564     if(!(borderType&BORDER_ISOLATED))
1565         src.locateROI( wsz, ofs );
1566     borderType = (borderType&~BORDER_ISOLATED);
1567
1568     Ptr<FilterEngine> f = createBoxFilter( src.type(), dst.type(),
1569                         ksize, anchor, normalize, borderType );
1570
1571     f->apply( src, dst, wsz, ofs );
1572 }
1573
1574
1575 void cv::blur( InputArray src, OutputArray dst,
1576            Size ksize, Point anchor, int borderType )
1577 {
1578     CV_INSTRUMENT_REGION()
1579
1580     boxFilter( src, dst, -1, ksize, anchor, true, borderType );
1581 }
1582
1583
1584 /****************************************************************************************\
1585                                     Squared Box Filter
1586 \****************************************************************************************/
1587
1588 namespace cv
1589 {
1590
1591 template<typename T, typename ST>
1592 struct SqrRowSum :
1593         public BaseRowFilter
1594 {
1595     SqrRowSum( int _ksize, int _anchor ) :
1596         BaseRowFilter()
1597     {
1598         ksize = _ksize;
1599         anchor = _anchor;
1600     }
1601
1602     virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
1603     {
1604         const T* S = (const T*)src;
1605         ST* D = (ST*)dst;
1606         int i = 0, k, ksz_cn = ksize*cn;
1607
1608         width = (width - 1)*cn;
1609         for( k = 0; k < cn; k++, S++, D++ )
1610         {
1611             ST s = 0;
1612             for( i = 0; i < ksz_cn; i += cn )
1613             {
1614                 ST val = (ST)S[i];
1615                 s += val*val;
1616             }
1617             D[0] = s;
1618             for( i = 0; i < width; i += cn )
1619             {
1620                 ST val0 = (ST)S[i], val1 = (ST)S[i + ksz_cn];
1621                 s += val1*val1 - val0*val0;
1622                 D[i+cn] = s;
1623             }
1624         }
1625     }
1626 };
1627
1628 static Ptr<BaseRowFilter> getSqrRowSumFilter(int srcType, int sumType, int ksize, int anchor)
1629 {
1630     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
1631     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
1632
1633     if( anchor < 0 )
1634         anchor = ksize/2;
1635
1636     if( sdepth == CV_8U && ddepth == CV_32S )
1637         return makePtr<SqrRowSum<uchar, int> >(ksize, anchor);
1638     if( sdepth == CV_8U && ddepth == CV_64F )
1639         return makePtr<SqrRowSum<uchar, double> >(ksize, anchor);
1640     if( sdepth == CV_16U && ddepth == CV_64F )
1641         return makePtr<SqrRowSum<ushort, double> >(ksize, anchor);
1642     if( sdepth == CV_16S && ddepth == CV_64F )
1643         return makePtr<SqrRowSum<short, double> >(ksize, anchor);
1644     if( sdepth == CV_32F && ddepth == CV_64F )
1645         return makePtr<SqrRowSum<float, double> >(ksize, anchor);
1646     if( sdepth == CV_64F && ddepth == CV_64F )
1647         return makePtr<SqrRowSum<double, double> >(ksize, anchor);
1648
1649     CV_Error_( CV_StsNotImplemented,
1650               ("Unsupported combination of source format (=%d), and buffer format (=%d)",
1651                srcType, sumType));
1652
1653     return Ptr<BaseRowFilter>();
1654 }
1655
1656 }
1657
1658 void cv::sqrBoxFilter( InputArray _src, OutputArray _dst, int ddepth,
1659                        Size ksize, Point anchor,
1660                        bool normalize, int borderType )
1661 {
1662     CV_INSTRUMENT_REGION()
1663
1664     int srcType = _src.type(), sdepth = CV_MAT_DEPTH(srcType), cn = CV_MAT_CN(srcType);
1665     Size size = _src.size();
1666
1667     if( ddepth < 0 )
1668         ddepth = sdepth < CV_32F ? CV_32F : CV_64F;
1669
1670     if( borderType != BORDER_CONSTANT && normalize )
1671     {
1672         if( size.height == 1 )
1673             ksize.height = 1;
1674         if( size.width == 1 )
1675             ksize.width = 1;
1676     }
1677
1678     CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2,
1679                ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize, true))
1680
1681     int sumDepth = CV_64F;
1682     if( sdepth == CV_8U )
1683         sumDepth = CV_32S;
1684     int sumType = CV_MAKETYPE( sumDepth, cn ), dstType = CV_MAKETYPE(ddepth, cn);
1685
1686     Mat src = _src.getMat();
1687     _dst.create( size, dstType );
1688     Mat dst = _dst.getMat();
1689
1690     Ptr<BaseRowFilter> rowFilter = getSqrRowSumFilter(srcType, sumType, ksize.width, anchor.x );
1691     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
1692                                                             dstType, ksize.height, anchor.y,
1693                                                             normalize ? 1./(ksize.width*ksize.height) : 1);
1694
1695     Ptr<FilterEngine> f = makePtr<FilterEngine>(Ptr<BaseFilter>(), rowFilter, columnFilter,
1696                                                 srcType, dstType, sumType, borderType );
1697     Point ofs;
1698     Size wsz(src.cols, src.rows);
1699     src.locateROI( wsz, ofs );
1700
1701     f->apply( src, dst, wsz, ofs );
1702 }
1703
1704
1705 /****************************************************************************************\
1706                                      Gaussian Blur
1707 \****************************************************************************************/
1708
1709 cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )
1710 {
1711     const int SMALL_GAUSSIAN_SIZE = 7;
1712     static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
1713     {
1714         {1.f},
1715         {0.25f, 0.5f, 0.25f},
1716         {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
1717         {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
1718     };
1719
1720     const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
1721         small_gaussian_tab[n>>1] : 0;
1722
1723     CV_Assert( ktype == CV_32F || ktype == CV_64F );
1724     Mat kernel(n, 1, ktype);
1725     float* cf = kernel.ptr<float>();
1726     double* cd = kernel.ptr<double>();
1727
1728     double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
1729     double scale2X = -0.5/(sigmaX*sigmaX);
1730     double sum = 0;
1731
1732     int i;
1733     for( i = 0; i < n; i++ )
1734     {
1735         double x = i - (n-1)*0.5;
1736         double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);
1737         if( ktype == CV_32F )
1738         {
1739             cf[i] = (float)t;
1740             sum += cf[i];
1741         }
1742         else
1743         {
1744             cd[i] = t;
1745             sum += cd[i];
1746         }
1747     }
1748
1749     sum = 1./sum;
1750     for( i = 0; i < n; i++ )
1751     {
1752         if( ktype == CV_32F )
1753             cf[i] = (float)(cf[i]*sum);
1754         else
1755             cd[i] *= sum;
1756     }
1757
1758     return kernel;
1759 }
1760
1761 namespace cv {
1762
1763 static void createGaussianKernels( Mat & kx, Mat & ky, int type, Size ksize,
1764                                    double sigma1, double sigma2 )
1765 {
1766     int depth = CV_MAT_DEPTH(type);
1767     if( sigma2 <= 0 )
1768         sigma2 = sigma1;
1769
1770     // automatic detection of kernel size from sigma
1771     if( ksize.width <= 0 && sigma1 > 0 )
1772         ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
1773     if( ksize.height <= 0 && sigma2 > 0 )
1774         ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
1775
1776     CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 &&
1777         ksize.height > 0 && ksize.height % 2 == 1 );
1778
1779     sigma1 = std::max( sigma1, 0. );
1780     sigma2 = std::max( sigma2, 0. );
1781
1782     kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) );
1783     if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )
1784         ky = kx;
1785     else
1786         ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) );
1787 }
1788
1789 }
1790
1791 cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize,
1792                                         double sigma1, double sigma2,
1793                                         int borderType )
1794 {
1795     Mat kx, ky;
1796     createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
1797
1798     return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );
1799 }
1800
1801 namespace cv
1802 {
1803 #ifdef HAVE_OPENCL
1804
1805 static bool ocl_GaussianBlur_8UC1(InputArray _src, OutputArray _dst, Size ksize, int ddepth,
1806                                   InputArray _kernelX, InputArray _kernelY, int borderType)
1807 {
1808     const ocl::Device & dev = ocl::Device::getDefault();
1809     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
1810
1811     if ( !(dev.isIntel() && (type == CV_8UC1) &&
1812          (_src.offset() == 0) && (_src.step() % 4 == 0) &&
1813          ((ksize.width == 5 && (_src.cols() % 4 == 0)) ||
1814          (ksize.width == 3 && (_src.cols() % 16 == 0) && (_src.rows() % 2 == 0)))) )
1815         return false;
1816
1817     Mat kernelX = _kernelX.getMat().reshape(1, 1);
1818     if (kernelX.cols % 2 != 1)
1819         return false;
1820     Mat kernelY = _kernelY.getMat().reshape(1, 1);
1821     if (kernelY.cols % 2 != 1)
1822         return false;
1823
1824     if (ddepth < 0)
1825         ddepth = sdepth;
1826
1827     Size size = _src.size();
1828     size_t globalsize[2] = { 0, 0 };
1829     size_t localsize[2] = { 0, 0 };
1830
1831     if (ksize.width == 3)
1832     {
1833         globalsize[0] = size.width / 16;
1834         globalsize[1] = size.height / 2;
1835     }
1836     else if (ksize.width == 5)
1837     {
1838         globalsize[0] = size.width / 4;
1839         globalsize[1] = size.height / 1;
1840     }
1841
1842     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
1843     char build_opts[1024];
1844     sprintf(build_opts, "-D %s %s%s", borderMap[borderType],
1845             ocl::kernelToStr(kernelX, CV_32F, "KERNEL_MATRIX_X").c_str(),
1846             ocl::kernelToStr(kernelY, CV_32F, "KERNEL_MATRIX_Y").c_str());
1847
1848     ocl::Kernel kernel;
1849
1850     if (ksize.width == 3)
1851         kernel.create("gaussianBlur3x3_8UC1_cols16_rows2", cv::ocl::imgproc::gaussianBlur3x3_oclsrc, build_opts);
1852     else if (ksize.width == 5)
1853         kernel.create("gaussianBlur5x5_8UC1_cols4", cv::ocl::imgproc::gaussianBlur5x5_oclsrc, build_opts);
1854
1855     if (kernel.empty())
1856         return false;
1857
1858     UMat src = _src.getUMat();
1859     _dst.create(size, CV_MAKETYPE(ddepth, cn));
1860     if (!(_dst.offset() == 0 && _dst.step() % 4 == 0))
1861         return false;
1862     UMat dst = _dst.getUMat();
1863
1864     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
1865     idxArg = kernel.set(idxArg, (int)src.step);
1866     idxArg = kernel.set(idxArg, ocl::KernelArg::PtrWriteOnly(dst));
1867     idxArg = kernel.set(idxArg, (int)dst.step);
1868     idxArg = kernel.set(idxArg, (int)dst.rows);
1869     idxArg = kernel.set(idxArg, (int)dst.cols);
1870
1871     return kernel.run(2, globalsize, (localsize[0] == 0) ? NULL : localsize, false);
1872 }
1873
1874 #endif
1875
1876 #ifdef HAVE_OPENVX
1877
1878 namespace ovx {
1879     template <> inline bool skipSmallImages<VX_KERNEL_GAUSSIAN_3x3>(int w, int h) { return w*h < 320 * 240; }
1880 }
1881 static bool openvx_gaussianBlur(InputArray _src, OutputArray _dst, Size ksize,
1882                                 double sigma1, double sigma2, int borderType)
1883 {
1884     if (sigma2 <= 0)
1885         sigma2 = sigma1;
1886     // automatic detection of kernel size from sigma
1887     if (ksize.width <= 0 && sigma1 > 0)
1888         ksize.width = cvRound(sigma1*6 + 1) | 1;
1889     if (ksize.height <= 0 && sigma2 > 0)
1890         ksize.height = cvRound(sigma2*6 + 1) | 1;
1891
1892     if (_src.type() != CV_8UC1 ||
1893         _src.cols() < 3 || _src.rows() < 3 ||
1894         ksize.width != 3 || ksize.height != 3)
1895         return false;
1896
1897     sigma1 = std::max(sigma1, 0.);
1898     sigma2 = std::max(sigma2, 0.);
1899
1900     if (!(sigma1 == 0.0 || (sigma1 - 0.8) < DBL_EPSILON) || !(sigma2 == 0.0 || (sigma2 - 0.8) < DBL_EPSILON) ||
1901         ovx::skipSmallImages<VX_KERNEL_GAUSSIAN_3x3>(_src.cols(), _src.rows()))
1902         return false;
1903
1904     Mat src = _src.getMat();
1905     Mat dst = _dst.getMat();
1906
1907     if ((borderType & BORDER_ISOLATED) == 0 && src.isSubmatrix())
1908         return false; //Process isolated borders only
1909     vx_enum border;
1910     switch (borderType & ~BORDER_ISOLATED)
1911     {
1912     case BORDER_CONSTANT:
1913         border = VX_BORDER_CONSTANT;
1914         break;
1915     case BORDER_REPLICATE:
1916         border = VX_BORDER_REPLICATE;
1917         break;
1918     default:
1919         return false;
1920     }
1921
1922     try
1923     {
1924         ivx::Context ctx = ovx::getOpenVXContext();
1925
1926         Mat a;
1927         if (dst.data != src.data)
1928             a = src;
1929         else
1930             src.copyTo(a);
1931
1932         ivx::Image
1933             ia = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
1934                 ivx::Image::createAddressing(a.cols, a.rows, 1, (vx_int32)(a.step)), a.data),
1935             ib = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
1936                 ivx::Image::createAddressing(dst.cols, dst.rows, 1, (vx_int32)(dst.step)), dst.data);
1937
1938         //ATTENTION: VX_CONTEXT_IMMEDIATE_BORDER attribute change could lead to strange issues in multi-threaded environments
1939         //since OpenVX standart says nothing about thread-safety for now
1940         ivx::border_t prevBorder = ctx.immediateBorder();
1941         ctx.setImmediateBorder(border, (vx_uint8)(0));
1942         ivx::IVX_CHECK_STATUS(vxuGaussian3x3(ctx, ia, ib));
1943         ctx.setImmediateBorder(prevBorder);
1944     }
1945     catch (ivx::RuntimeError & e)
1946     {
1947         VX_DbgThrow(e.what());
1948     }
1949     catch (ivx::WrapperError & e)
1950     {
1951         VX_DbgThrow(e.what());
1952     }
1953     return true;
1954 }
1955
1956 #endif
1957
1958 #ifdef HAVE_IPP
1959 #if IPP_VERSION_X100 == 201702  // IW 2017u2 has bug which doesn't allow use of partial inMem with tiling
1960 #define IPP_GAUSSIANBLUR_PARALLEL 0
1961 #else
1962 #define IPP_GAUSSIANBLUR_PARALLEL 1
1963 #endif
1964
1965 #ifdef HAVE_IPP_IW
1966
1967 class ipp_gaussianBlurParallel: public ParallelLoopBody
1968 {
1969 public:
1970     ipp_gaussianBlurParallel(::ipp::IwiImage &src, ::ipp::IwiImage &dst, int kernelSize, float sigma, ::ipp::IwiBorderType &border, bool *pOk):
1971         m_src(src), m_dst(dst), m_kernelSize(kernelSize), m_sigma(sigma), m_border(border), m_pOk(pOk) {
1972         *m_pOk = true;
1973     }
1974     ~ipp_gaussianBlurParallel()
1975     {
1976     }
1977
1978     virtual void operator() (const Range& range) const
1979     {
1980         CV_INSTRUMENT_REGION_IPP()
1981
1982         if(!*m_pOk)
1983             return;
1984
1985         try
1986         {
1987             ::ipp::IwiTile tile = ::ipp::IwiRoi(0, range.start, m_dst.m_size.width, range.end - range.start);
1988             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterGaussian, m_src, m_dst, m_kernelSize, m_sigma, ::ipp::IwDefault(), m_border, tile);
1989         }
1990         catch(::ipp::IwException e)
1991         {
1992             *m_pOk = false;
1993             return;
1994         }
1995     }
1996 private:
1997     ::ipp::IwiImage &m_src;
1998     ::ipp::IwiImage &m_dst;
1999
2000     int m_kernelSize;
2001     float m_sigma;
2002     ::ipp::IwiBorderType &m_border;
2003
2004     volatile bool *m_pOk;
2005     const ipp_gaussianBlurParallel& operator= (const ipp_gaussianBlurParallel&);
2006 };
2007
2008 #endif
2009
2010 static bool ipp_GaussianBlur(InputArray _src, OutputArray _dst, Size ksize,
2011                    double sigma1, double sigma2, int borderType )
2012 {
2013 #ifdef HAVE_IPP_IW
2014     CV_INSTRUMENT_REGION_IPP()
2015
2016 #if IPP_VERSION_X100 < 201800 && ((defined _MSC_VER && defined _M_IX86) || (defined __GNUC__ && defined __i386__))
2017     CV_UNUSED(_src); CV_UNUSED(_dst); CV_UNUSED(ksize); CV_UNUSED(sigma1); CV_UNUSED(sigma2); CV_UNUSED(borderType);
2018     return false; // bug on ia32
2019 #else
2020     if(sigma1 != sigma2)
2021         return false;
2022
2023     if(sigma1 < FLT_EPSILON)
2024         return false;
2025
2026     if(ksize.width != ksize.height)
2027         return false;
2028
2029     // Acquire data and begin processing
2030     try
2031     {
2032         Mat src = _src.getMat();
2033         Mat dst = _dst.getMat();
2034         ::ipp::IwiImage       iwSrc      = ippiGetImage(src);
2035         ::ipp::IwiImage       iwDst      = ippiGetImage(dst);
2036         ::ipp::IwiBorderSize  borderSize = ::ipp::iwiSizeToBorderSize(ippiGetSize(ksize));
2037         ::ipp::IwiBorderType  ippBorder(ippiGetBorder(iwSrc, borderType, borderSize));
2038         if(!ippBorder)
2039             return false;
2040
2041         const int threads = ippiSuggestThreadsNum(iwDst, 2);
2042         if(IPP_GAUSSIANBLUR_PARALLEL && threads > 1) {
2043             bool ok;
2044             ipp_gaussianBlurParallel invoker(iwSrc, iwDst, ksize.width, (float) sigma1, ippBorder, &ok);
2045
2046             if(!ok)
2047                 return false;
2048             const Range range(0, (int) iwDst.m_size.height);
2049             parallel_for_(range, invoker, threads*4);
2050
2051             if(!ok)
2052                 return false;
2053         } else {
2054             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterGaussian, iwSrc, iwDst, ksize.width, sigma1, ::ipp::IwDefault(), ippBorder);
2055         }
2056     }
2057     catch (::ipp::IwException ex)
2058     {
2059         return false;
2060     }
2061
2062     return true;
2063 #endif
2064 #else
2065     CV_UNUSED(_src); CV_UNUSED(_dst); CV_UNUSED(ksize); CV_UNUSED(sigma1); CV_UNUSED(sigma2); CV_UNUSED(borderType);
2066     return false;
2067 #endif
2068 }
2069 #endif
2070 }
2071
2072 void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
2073                    double sigma1, double sigma2,
2074                    int borderType )
2075 {
2076     CV_INSTRUMENT_REGION()
2077
2078     int type = _src.type();
2079     Size size = _src.size();
2080     _dst.create( size, type );
2081
2082     if( borderType != BORDER_CONSTANT && (borderType & BORDER_ISOLATED) != 0 )
2083     {
2084         if( size.height == 1 )
2085             ksize.height = 1;
2086         if( size.width == 1 )
2087             ksize.width = 1;
2088     }
2089
2090     if( ksize.width == 1 && ksize.height == 1 )
2091     {
2092         _src.copyTo(_dst);
2093         return;
2094     }
2095
2096     CV_OVX_RUN(true,
2097                openvx_gaussianBlur(_src, _dst, ksize, sigma1, sigma2, borderType))
2098
2099 #ifdef HAVE_TEGRA_OPTIMIZATION
2100     Mat src = _src.getMat();
2101     Mat dst = _dst.getMat();
2102     if(sigma1 == 0 && sigma2 == 0 && tegra::useTegra() && tegra::gaussian(src, dst, ksize, borderType))
2103         return;
2104 #endif
2105     bool useOpenCL = (ocl::useOpenCL() && _dst.isUMat() && _src.dims() <= 2 &&
2106                ((ksize.width == 3 && ksize.height == 3) ||
2107                (ksize.width == 5 && ksize.height == 5)) &&
2108                _src.rows() > ksize.height && _src.cols() > ksize.width);
2109     (void)useOpenCL;
2110
2111     CV_IPP_RUN(!useOpenCL, ipp_GaussianBlur( _src,  _dst,  ksize, sigma1,  sigma2, borderType));
2112
2113     Mat kx, ky;
2114     createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
2115
2116     CV_OCL_RUN(useOpenCL, ocl_GaussianBlur_8UC1(_src, _dst, ksize, CV_MAT_DEPTH(type), kx, ky, borderType));
2117
2118     sepFilter2D(_src, _dst, CV_MAT_DEPTH(type), kx, ky, Point(-1,-1), 0, borderType );
2119 }
2120
2121 /****************************************************************************************\
2122                                       Median Filter
2123 \****************************************************************************************/
2124
2125 namespace cv
2126 {
2127 typedef ushort HT;
2128
2129 /**
2130  * This structure represents a two-tier histogram. The first tier (known as the
2131  * "coarse" level) is 4 bit wide and the second tier (known as the "fine" level)
2132  * is 8 bit wide. Pixels inserted in the fine level also get inserted into the
2133  * coarse bucket designated by the 4 MSBs of the fine bucket value.
2134  *
2135  * The structure is aligned on 16 bits, which is a prerequisite for SIMD
2136  * instructions. Each bucket is 16 bit wide, which means that extra care must be
2137  * taken to prevent overflow.
2138  */
2139 typedef struct
2140 {
2141     HT coarse[16];
2142     HT fine[16][16];
2143 } Histogram;
2144
2145
2146 #if CV_SIMD128
2147
2148 static inline void histogram_add_simd( const HT x[16], HT y[16] )
2149 {
2150     v_store(y, v_load(x) + v_load(y));
2151     v_store(y + 8, v_load(x + 8) + v_load(y + 8));
2152 }
2153
2154 static inline void histogram_sub_simd( const HT x[16], HT y[16] )
2155 {
2156     v_store(y, v_load(y) - v_load(x));
2157     v_store(y + 8, v_load(y + 8) - v_load(x + 8));
2158 }
2159
2160 #endif
2161
2162
2163 static inline void histogram_add( const HT x[16], HT y[16] )
2164 {
2165     int i;
2166     for( i = 0; i < 16; ++i )
2167         y[i] = (HT)(y[i] + x[i]);
2168 }
2169
2170 static inline void histogram_sub( const HT x[16], HT y[16] )
2171 {
2172     int i;
2173     for( i = 0; i < 16; ++i )
2174         y[i] = (HT)(y[i] - x[i]);
2175 }
2176
2177 static inline void histogram_muladd( int a, const HT x[16],
2178         HT y[16] )
2179 {
2180     for( int i = 0; i < 16; ++i )
2181         y[i] = (HT)(y[i] + a * x[i]);
2182 }
2183
2184 static void
2185 medianBlur_8u_O1( const Mat& _src, Mat& _dst, int ksize )
2186 {
2187 /**
2188  * HOP is short for Histogram OPeration. This macro makes an operation \a op on
2189  * histogram \a h for pixel value \a x. It takes care of handling both levels.
2190  */
2191 #define HOP(h,x,op) \
2192     h.coarse[x>>4] op, \
2193     *((HT*)h.fine + x) op
2194
2195 #define COP(c,j,x,op) \
2196     h_coarse[ 16*(n*c+j) + (x>>4) ] op, \
2197     h_fine[ 16 * (n*(16*c+(x>>4)) + j) + (x & 0xF) ] op
2198
2199     int cn = _dst.channels(), m = _dst.rows, r = (ksize-1)/2;
2200     CV_Assert(cn > 0 && cn <= 4);
2201     size_t sstep = _src.step, dstep = _dst.step;
2202     Histogram CV_DECL_ALIGNED(16) H[4];
2203     HT CV_DECL_ALIGNED(16) luc[4][16];
2204
2205     int STRIPE_SIZE = std::min( _dst.cols, 512/cn );
2206
2207     std::vector<HT> _h_coarse(1 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
2208     std::vector<HT> _h_fine(16 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
2209     HT* h_coarse = alignPtr(&_h_coarse[0], 16);
2210     HT* h_fine = alignPtr(&_h_fine[0], 16);
2211 #if CV_SIMD128
2212     volatile bool useSIMD = hasSIMD128();
2213 #endif
2214
2215     for( int x = 0; x < _dst.cols; x += STRIPE_SIZE )
2216     {
2217         int i, j, k, c, n = std::min(_dst.cols - x, STRIPE_SIZE) + r*2;
2218         const uchar* src = _src.ptr() + x*cn;
2219         uchar* dst = _dst.ptr() + (x - r)*cn;
2220
2221         memset( h_coarse, 0, 16*n*cn*sizeof(h_coarse[0]) );
2222         memset( h_fine, 0, 16*16*n*cn*sizeof(h_fine[0]) );
2223
2224         // First row initialization
2225         for( c = 0; c < cn; c++ )
2226         {
2227             for( j = 0; j < n; j++ )
2228                 COP( c, j, src[cn*j+c], += (cv::HT)(r+2) );
2229
2230             for( i = 1; i < r; i++ )
2231             {
2232                 const uchar* p = src + sstep*std::min(i, m-1);
2233                 for ( j = 0; j < n; j++ )
2234                     COP( c, j, p[cn*j+c], ++ );
2235             }
2236         }
2237
2238         for( i = 0; i < m; i++ )
2239         {
2240             const uchar* p0 = src + sstep * std::max( 0, i-r-1 );
2241             const uchar* p1 = src + sstep * std::min( m-1, i+r );
2242
2243             memset( H, 0, cn*sizeof(H[0]) );
2244             memset( luc, 0, cn*sizeof(luc[0]) );
2245             for( c = 0; c < cn; c++ )
2246             {
2247                 // Update column histograms for the entire row.
2248                 for( j = 0; j < n; j++ )
2249                 {
2250                     COP( c, j, p0[j*cn + c], -- );
2251                     COP( c, j, p1[j*cn + c], ++ );
2252                 }
2253
2254                 // First column initialization
2255                 for( k = 0; k < 16; ++k )
2256                     histogram_muladd( 2*r+1, &h_fine[16*n*(16*c+k)], &H[c].fine[k][0] );
2257
2258 #if CV_SIMD128
2259                 if( useSIMD )
2260                 {
2261                     for( j = 0; j < 2*r; ++j )
2262                         histogram_add_simd( &h_coarse[16*(n*c+j)], H[c].coarse );
2263
2264                     for( j = r; j < n-r; j++ )
2265                     {
2266                         int t = 2*r*r + 2*r, b, sum = 0;
2267                         HT* segment;
2268
2269                         histogram_add_simd( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
2270
2271                         // Find median at coarse level
2272                         for ( k = 0; k < 16 ; ++k )
2273                         {
2274                             sum += H[c].coarse[k];
2275                             if ( sum > t )
2276                             {
2277                                 sum -= H[c].coarse[k];
2278                                 break;
2279                             }
2280                         }
2281                         CV_Assert( k < 16 );
2282
2283                         /* Update corresponding histogram segment */
2284                         if ( luc[c][k] <= j-r )
2285                         {
2286                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
2287                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
2288                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
2289
2290                             if ( luc[c][k] < j+r+1 )
2291                             {
2292                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
2293                                 luc[c][k] = (HT)(j+r+1);
2294                             }
2295                         }
2296                         else
2297                         {
2298                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
2299                             {
2300                                 histogram_sub_simd( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
2301                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
2302                             }
2303                         }
2304
2305                         histogram_sub_simd( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
2306
2307                         /* Find median in segment */
2308                         segment = H[c].fine[k];
2309                         for ( b = 0; b < 16 ; b++ )
2310                         {
2311                             sum += segment[b];
2312                             if ( sum > t )
2313                             {
2314                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
2315                                 break;
2316                             }
2317                         }
2318                         CV_Assert( b < 16 );
2319                     }
2320                 }
2321                 else
2322 #endif
2323                 {
2324                     for( j = 0; j < 2*r; ++j )
2325                         histogram_add( &h_coarse[16*(n*c+j)], H[c].coarse );
2326
2327                     for( j = r; j < n-r; j++ )
2328                     {
2329                         int t = 2*r*r + 2*r, b, sum = 0;
2330                         HT* segment;
2331
2332                         histogram_add( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
2333
2334                         // Find median at coarse level
2335                         for ( k = 0; k < 16 ; ++k )
2336                         {
2337                             sum += H[c].coarse[k];
2338                             if ( sum > t )
2339                             {
2340                                 sum -= H[c].coarse[k];
2341                                 break;
2342                             }
2343                         }
2344                         CV_Assert( k < 16 );
2345
2346                         /* Update corresponding histogram segment */
2347                         if ( luc[c][k] <= j-r )
2348                         {
2349                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
2350                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
2351                                 histogram_add( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
2352
2353                             if ( luc[c][k] < j+r+1 )
2354                             {
2355                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
2356                                 luc[c][k] = (HT)(j+r+1);
2357                             }
2358                         }
2359                         else
2360                         {
2361                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
2362                             {
2363                                 histogram_sub( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
2364                                 histogram_add( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
2365                             }
2366                         }
2367
2368                         histogram_sub( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
2369
2370                         /* Find median in segment */
2371                         segment = H[c].fine[k];
2372                         for ( b = 0; b < 16 ; b++ )
2373                         {
2374                             sum += segment[b];
2375                             if ( sum > t )
2376                             {
2377                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
2378                                 break;
2379                             }
2380                         }
2381                         CV_Assert( b < 16 );
2382                     }
2383                 }
2384             }
2385         }
2386     }
2387
2388 #undef HOP
2389 #undef COP
2390 }
2391
2392 static void
2393 medianBlur_8u_Om( const Mat& _src, Mat& _dst, int m )
2394 {
2395     #define N  16
2396     int     zone0[4][N];
2397     int     zone1[4][N*N];
2398     int     x, y;
2399     int     n2 = m*m/2;
2400     Size    size = _dst.size();
2401     const uchar* src = _src.ptr();
2402     uchar*  dst = _dst.ptr();
2403     int     src_step = (int)_src.step, dst_step = (int)_dst.step;
2404     int     cn = _src.channels();
2405     const uchar*  src_max = src + size.height*src_step;
2406     CV_Assert(cn > 0 && cn <= 4);
2407
2408     #define UPDATE_ACC01( pix, cn, op ) \
2409     {                                   \
2410         int p = (pix);                  \
2411         zone1[cn][p] op;                \
2412         zone0[cn][p >> 4] op;           \
2413     }
2414
2415     //CV_Assert( size.height >= nx && size.width >= nx );
2416     for( x = 0; x < size.width; x++, src += cn, dst += cn )
2417     {
2418         uchar* dst_cur = dst;
2419         const uchar* src_top = src;
2420         const uchar* src_bottom = src;
2421         int k, c;
2422         int src_step1 = src_step, dst_step1 = dst_step;
2423
2424         if( x % 2 != 0 )
2425         {
2426             src_bottom = src_top += src_step*(size.height-1);
2427             dst_cur += dst_step*(size.height-1);
2428             src_step1 = -src_step1;
2429             dst_step1 = -dst_step1;
2430         }
2431
2432         // init accumulator
2433         memset( zone0, 0, sizeof(zone0[0])*cn );
2434         memset( zone1, 0, sizeof(zone1[0])*cn );
2435
2436         for( y = 0; y <= m/2; y++ )
2437         {
2438             for( c = 0; c < cn; c++ )
2439             {
2440                 if( y > 0 )
2441                 {
2442                     for( k = 0; k < m*cn; k += cn )
2443                         UPDATE_ACC01( src_bottom[k+c], c, ++ );
2444                 }
2445                 else
2446                 {
2447                     for( k = 0; k < m*cn; k += cn )
2448                         UPDATE_ACC01( src_bottom[k+c], c, += m/2+1 );
2449                 }
2450             }
2451
2452             if( (src_step1 > 0 && y < size.height-1) ||
2453                 (src_step1 < 0 && size.height-y-1 > 0) )
2454                 src_bottom += src_step1;
2455         }
2456
2457         for( y = 0; y < size.height; y++, dst_cur += dst_step1 )
2458         {
2459             // find median
2460             for( c = 0; c < cn; c++ )
2461             {
2462                 int s = 0;
2463                 for( k = 0; ; k++ )
2464                 {
2465                     int t = s + zone0[c][k];
2466                     if( t > n2 ) break;
2467                     s = t;
2468                 }
2469
2470                 for( k *= N; ;k++ )
2471                 {
2472                     s += zone1[c][k];
2473                     if( s > n2 ) break;
2474                 }
2475
2476                 dst_cur[c] = (uchar)k;
2477             }
2478
2479             if( y+1 == size.height )
2480                 break;
2481
2482             if( cn == 1 )
2483             {
2484                 for( k = 0; k < m; k++ )
2485                 {
2486                     int p = src_top[k];
2487                     int q = src_bottom[k];
2488                     zone1[0][p]--;
2489                     zone0[0][p>>4]--;
2490                     zone1[0][q]++;
2491                     zone0[0][q>>4]++;
2492                 }
2493             }
2494             else if( cn == 3 )
2495             {
2496                 for( k = 0; k < m*3; k += 3 )
2497                 {
2498                     UPDATE_ACC01( src_top[k], 0, -- );
2499                     UPDATE_ACC01( src_top[k+1], 1, -- );
2500                     UPDATE_ACC01( src_top[k+2], 2, -- );
2501
2502                     UPDATE_ACC01( src_bottom[k], 0, ++ );
2503                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
2504                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
2505                 }
2506             }
2507             else
2508             {
2509                 assert( cn == 4 );
2510                 for( k = 0; k < m*4; k += 4 )
2511                 {
2512                     UPDATE_ACC01( src_top[k], 0, -- );
2513                     UPDATE_ACC01( src_top[k+1], 1, -- );
2514                     UPDATE_ACC01( src_top[k+2], 2, -- );
2515                     UPDATE_ACC01( src_top[k+3], 3, -- );
2516
2517                     UPDATE_ACC01( src_bottom[k], 0, ++ );
2518                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
2519                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
2520                     UPDATE_ACC01( src_bottom[k+3], 3, ++ );
2521                 }
2522             }
2523
2524             if( (src_step1 > 0 && src_bottom + src_step1 < src_max) ||
2525                 (src_step1 < 0 && src_bottom + src_step1 >= src) )
2526                 src_bottom += src_step1;
2527
2528             if( y >= m/2 )
2529                 src_top += src_step1;
2530         }
2531     }
2532 #undef N
2533 #undef UPDATE_ACC
2534 }
2535
2536
2537 struct MinMax8u
2538 {
2539     typedef uchar value_type;
2540     typedef int arg_type;
2541     enum { SIZE = 1 };
2542     arg_type load(const uchar* ptr) { return *ptr; }
2543     void store(uchar* ptr, arg_type val) { *ptr = (uchar)val; }
2544     void operator()(arg_type& a, arg_type& b) const
2545     {
2546         int t = CV_FAST_CAST_8U(a - b);
2547         b += t; a -= t;
2548     }
2549 };
2550
2551 struct MinMax16u
2552 {
2553     typedef ushort value_type;
2554     typedef int arg_type;
2555     enum { SIZE = 1 };
2556     arg_type load(const ushort* ptr) { return *ptr; }
2557     void store(ushort* ptr, arg_type val) { *ptr = (ushort)val; }
2558     void operator()(arg_type& a, arg_type& b) const
2559     {
2560         arg_type t = a;
2561         a = std::min(a, b);
2562         b = std::max(b, t);
2563     }
2564 };
2565
2566 struct MinMax16s
2567 {
2568     typedef short value_type;
2569     typedef int arg_type;
2570     enum { SIZE = 1 };
2571     arg_type load(const short* ptr) { return *ptr; }
2572     void store(short* ptr, arg_type val) { *ptr = (short)val; }
2573     void operator()(arg_type& a, arg_type& b) const
2574     {
2575         arg_type t = a;
2576         a = std::min(a, b);
2577         b = std::max(b, t);
2578     }
2579 };
2580
2581 struct MinMax32f
2582 {
2583     typedef float value_type;
2584     typedef float arg_type;
2585     enum { SIZE = 1 };
2586     arg_type load(const float* ptr) { return *ptr; }
2587     void store(float* ptr, arg_type val) { *ptr = val; }
2588     void operator()(arg_type& a, arg_type& b) const
2589     {
2590         arg_type t = a;
2591         a = std::min(a, b);
2592         b = std::max(b, t);
2593     }
2594 };
2595
2596 #if CV_SIMD128
2597
2598 struct MinMaxVec8u
2599 {
2600     typedef uchar value_type;
2601     typedef v_uint8x16 arg_type;
2602     enum { SIZE = 16 };
2603     arg_type load(const uchar* ptr) { return v_load(ptr); }
2604     void store(uchar* ptr, const arg_type &val) { v_store(ptr, val); }
2605     void operator()(arg_type& a, arg_type& b) const
2606     {
2607         arg_type t = a;
2608         a = v_min(a, b);
2609         b = v_max(b, t);
2610     }
2611 };
2612
2613
2614 struct MinMaxVec16u
2615 {
2616     typedef ushort value_type;
2617     typedef v_uint16x8 arg_type;
2618     enum { SIZE = 8 };
2619     arg_type load(const ushort* ptr) { return v_load(ptr); }
2620     void store(ushort* ptr, const arg_type &val) { v_store(ptr, val); }
2621     void operator()(arg_type& a, arg_type& b) const
2622     {
2623         arg_type t = a;
2624         a = v_min(a, b);
2625         b = v_max(b, t);
2626     }
2627 };
2628
2629
2630 struct MinMaxVec16s
2631 {
2632     typedef short value_type;
2633     typedef v_int16x8 arg_type;
2634     enum { SIZE = 8 };
2635     arg_type load(const short* ptr) { return v_load(ptr); }
2636     void store(short* ptr, const arg_type &val) { v_store(ptr, val); }
2637     void operator()(arg_type& a, arg_type& b) const
2638     {
2639         arg_type t = a;
2640         a = v_min(a, b);
2641         b = v_max(b, t);
2642     }
2643 };
2644
2645
2646 struct MinMaxVec32f
2647 {
2648     typedef float value_type;
2649     typedef v_float32x4 arg_type;
2650     enum { SIZE = 4 };
2651     arg_type load(const float* ptr) { return v_load(ptr); }
2652     void store(float* ptr, const arg_type &val) { v_store(ptr, val); }
2653     void operator()(arg_type& a, arg_type& b) const
2654     {
2655         arg_type t = a;
2656         a = v_min(a, b);
2657         b = v_max(b, t);
2658     }
2659 };
2660
2661 #else
2662
2663 typedef MinMax8u MinMaxVec8u;
2664 typedef MinMax16u MinMaxVec16u;
2665 typedef MinMax16s MinMaxVec16s;
2666 typedef MinMax32f MinMaxVec32f;
2667
2668 #endif
2669
2670 template<class Op, class VecOp>
2671 static void
2672 medianBlur_SortNet( const Mat& _src, Mat& _dst, int m )
2673 {
2674     typedef typename Op::value_type T;
2675     typedef typename Op::arg_type WT;
2676     typedef typename VecOp::arg_type VT;
2677
2678     const T* src = _src.ptr<T>();
2679     T* dst = _dst.ptr<T>();
2680     int sstep = (int)(_src.step/sizeof(T));
2681     int dstep = (int)(_dst.step/sizeof(T));
2682     Size size = _dst.size();
2683     int i, j, k, cn = _src.channels();
2684     Op op;
2685     VecOp vop;
2686     volatile bool useSIMD = hasSIMD128();
2687
2688     if( m == 3 )
2689     {
2690         if( size.width == 1 || size.height == 1 )
2691         {
2692             int len = size.width + size.height - 1;
2693             int sdelta = size.height == 1 ? cn : sstep;
2694             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
2695             int ddelta = size.height == 1 ? cn : dstep;
2696
2697             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
2698                 for( j = 0; j < cn; j++, src++ )
2699                 {
2700                     WT p0 = src[i > 0 ? -sdelta : 0];
2701                     WT p1 = src[0];
2702                     WT p2 = src[i < len - 1 ? sdelta : 0];
2703
2704                     op(p0, p1); op(p1, p2); op(p0, p1);
2705                     dst[j] = (T)p1;
2706                 }
2707             return;
2708         }
2709
2710         size.width *= cn;
2711         for( i = 0; i < size.height; i++, dst += dstep )
2712         {
2713             const T* row0 = src + std::max(i - 1, 0)*sstep;
2714             const T* row1 = src + i*sstep;
2715             const T* row2 = src + std::min(i + 1, size.height-1)*sstep;
2716             int limit = useSIMD ? cn : size.width;
2717
2718             for(j = 0;; )
2719             {
2720                 for( ; j < limit; j++ )
2721                 {
2722                     int j0 = j >= cn ? j - cn : j;
2723                     int j2 = j < size.width - cn ? j + cn : j;
2724                     WT p0 = row0[j0], p1 = row0[j], p2 = row0[j2];
2725                     WT p3 = row1[j0], p4 = row1[j], p5 = row1[j2];
2726                     WT p6 = row2[j0], p7 = row2[j], p8 = row2[j2];
2727
2728                     op(p1, p2); op(p4, p5); op(p7, p8); op(p0, p1);
2729                     op(p3, p4); op(p6, p7); op(p1, p2); op(p4, p5);
2730                     op(p7, p8); op(p0, p3); op(p5, p8); op(p4, p7);
2731                     op(p3, p6); op(p1, p4); op(p2, p5); op(p4, p7);
2732                     op(p4, p2); op(p6, p4); op(p4, p2);
2733                     dst[j] = (T)p4;
2734                 }
2735
2736                 if( limit == size.width )
2737                     break;
2738
2739                 for( ; j <= size.width - VecOp::SIZE - cn; j += VecOp::SIZE )
2740                 {
2741                     VT p0 = vop.load(row0+j-cn), p1 = vop.load(row0+j), p2 = vop.load(row0+j+cn);
2742                     VT p3 = vop.load(row1+j-cn), p4 = vop.load(row1+j), p5 = vop.load(row1+j+cn);
2743                     VT p6 = vop.load(row2+j-cn), p7 = vop.load(row2+j), p8 = vop.load(row2+j+cn);
2744
2745                     vop(p1, p2); vop(p4, p5); vop(p7, p8); vop(p0, p1);
2746                     vop(p3, p4); vop(p6, p7); vop(p1, p2); vop(p4, p5);
2747                     vop(p7, p8); vop(p0, p3); vop(p5, p8); vop(p4, p7);
2748                     vop(p3, p6); vop(p1, p4); vop(p2, p5); vop(p4, p7);
2749                     vop(p4, p2); vop(p6, p4); vop(p4, p2);
2750                     vop.store(dst+j, p4);
2751                 }
2752
2753                 limit = size.width;
2754             }
2755         }
2756     }
2757     else if( m == 5 )
2758     {
2759         if( size.width == 1 || size.height == 1 )
2760         {
2761             int len = size.width + size.height - 1;
2762             int sdelta = size.height == 1 ? cn : sstep;
2763             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
2764             int ddelta = size.height == 1 ? cn : dstep;
2765
2766             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
2767                 for( j = 0; j < cn; j++, src++ )
2768                 {
2769                     int i1 = i > 0 ? -sdelta : 0;
2770                     int i0 = i > 1 ? -sdelta*2 : i1;
2771                     int i3 = i < len-1 ? sdelta : 0;
2772                     int i4 = i < len-2 ? sdelta*2 : i3;
2773                     WT p0 = src[i0], p1 = src[i1], p2 = src[0], p3 = src[i3], p4 = src[i4];
2774
2775                     op(p0, p1); op(p3, p4); op(p2, p3); op(p3, p4); op(p0, p2);
2776                     op(p2, p4); op(p1, p3); op(p1, p2);
2777                     dst[j] = (T)p2;
2778                 }
2779             return;
2780         }
2781
2782         size.width *= cn;
2783         for( i = 0; i < size.height; i++, dst += dstep )
2784         {
2785             const T* row[5];
2786             row[0] = src + std::max(i - 2, 0)*sstep;
2787             row[1] = src + std::max(i - 1, 0)*sstep;
2788             row[2] = src + i*sstep;
2789             row[3] = src + std::min(i + 1, size.height-1)*sstep;
2790             row[4] = src + std::min(i + 2, size.height-1)*sstep;
2791             int limit = useSIMD ? cn*2 : size.width;
2792
2793             for(j = 0;; )
2794             {
2795                 for( ; j < limit; j++ )
2796                 {
2797                     WT p[25];
2798                     int j1 = j >= cn ? j - cn : j;
2799                     int j0 = j >= cn*2 ? j - cn*2 : j1;
2800                     int j3 = j < size.width - cn ? j + cn : j;
2801                     int j4 = j < size.width - cn*2 ? j + cn*2 : j3;
2802                     for( k = 0; k < 5; k++ )
2803                     {
2804                         const T* rowk = row[k];
2805                         p[k*5] = rowk[j0]; p[k*5+1] = rowk[j1];
2806                         p[k*5+2] = rowk[j]; p[k*5+3] = rowk[j3];
2807                         p[k*5+4] = rowk[j4];
2808                     }
2809
2810                     op(p[1], p[2]); op(p[0], p[1]); op(p[1], p[2]); op(p[4], p[5]); op(p[3], p[4]);
2811                     op(p[4], p[5]); op(p[0], p[3]); op(p[2], p[5]); op(p[2], p[3]); op(p[1], p[4]);
2812                     op(p[1], p[2]); op(p[3], p[4]); op(p[7], p[8]); op(p[6], p[7]); op(p[7], p[8]);
2813                     op(p[10], p[11]); op(p[9], p[10]); op(p[10], p[11]); op(p[6], p[9]); op(p[8], p[11]);
2814                     op(p[8], p[9]); op(p[7], p[10]); op(p[7], p[8]); op(p[9], p[10]); op(p[0], p[6]);
2815                     op(p[4], p[10]); op(p[4], p[6]); op(p[2], p[8]); op(p[2], p[4]); op(p[6], p[8]);
2816                     op(p[1], p[7]); op(p[5], p[11]); op(p[5], p[7]); op(p[3], p[9]); op(p[3], p[5]);
2817                     op(p[7], p[9]); op(p[1], p[2]); op(p[3], p[4]); op(p[5], p[6]); op(p[7], p[8]);
2818                     op(p[9], p[10]); op(p[13], p[14]); op(p[12], p[13]); op(p[13], p[14]); op(p[16], p[17]);
2819                     op(p[15], p[16]); op(p[16], p[17]); op(p[12], p[15]); op(p[14], p[17]); op(p[14], p[15]);
2820                     op(p[13], p[16]); op(p[13], p[14]); op(p[15], p[16]); op(p[19], p[20]); op(p[18], p[19]);
2821                     op(p[19], p[20]); op(p[21], p[22]); op(p[23], p[24]); op(p[21], p[23]); op(p[22], p[24]);
2822                     op(p[22], p[23]); op(p[18], p[21]); op(p[20], p[23]); op(p[20], p[21]); op(p[19], p[22]);
2823                     op(p[22], p[24]); op(p[19], p[20]); op(p[21], p[22]); op(p[23], p[24]); op(p[12], p[18]);
2824                     op(p[16], p[22]); op(p[16], p[18]); op(p[14], p[20]); op(p[20], p[24]); op(p[14], p[16]);
2825                     op(p[18], p[20]); op(p[22], p[24]); op(p[13], p[19]); op(p[17], p[23]); op(p[17], p[19]);
2826                     op(p[15], p[21]); op(p[15], p[17]); op(p[19], p[21]); op(p[13], p[14]); op(p[15], p[16]);
2827                     op(p[17], p[18]); op(p[19], p[20]); op(p[21], p[22]); op(p[23], p[24]); op(p[0], p[12]);
2828                     op(p[8], p[20]); op(p[8], p[12]); op(p[4], p[16]); op(p[16], p[24]); op(p[12], p[16]);
2829                     op(p[2], p[14]); op(p[10], p[22]); op(p[10], p[14]); op(p[6], p[18]); op(p[6], p[10]);
2830                     op(p[10], p[12]); op(p[1], p[13]); op(p[9], p[21]); op(p[9], p[13]); op(p[5], p[17]);
2831                     op(p[13], p[17]); op(p[3], p[15]); op(p[11], p[23]); op(p[11], p[15]); op(p[7], p[19]);
2832                     op(p[7], p[11]); op(p[11], p[13]); op(p[11], p[12]);
2833                     dst[j] = (T)p[12];
2834                 }
2835
2836                 if( limit == size.width )
2837                     break;
2838
2839                 for( ; j <= size.width - VecOp::SIZE - cn*2; j += VecOp::SIZE )
2840                 {
2841                     VT p[25];
2842                     for( k = 0; k < 5; k++ )
2843                     {
2844                         const T* rowk = row[k];
2845                         p[k*5] = vop.load(rowk+j-cn*2); p[k*5+1] = vop.load(rowk+j-cn);
2846                         p[k*5+2] = vop.load(rowk+j); p[k*5+3] = vop.load(rowk+j+cn);
2847                         p[k*5+4] = vop.load(rowk+j+cn*2);
2848                     }
2849
2850                     vop(p[1], p[2]); vop(p[0], p[1]); vop(p[1], p[2]); vop(p[4], p[5]); vop(p[3], p[4]);
2851                     vop(p[4], p[5]); vop(p[0], p[3]); vop(p[2], p[5]); vop(p[2], p[3]); vop(p[1], p[4]);
2852                     vop(p[1], p[2]); vop(p[3], p[4]); vop(p[7], p[8]); vop(p[6], p[7]); vop(p[7], p[8]);
2853                     vop(p[10], p[11]); vop(p[9], p[10]); vop(p[10], p[11]); vop(p[6], p[9]); vop(p[8], p[11]);
2854                     vop(p[8], p[9]); vop(p[7], p[10]); vop(p[7], p[8]); vop(p[9], p[10]); vop(p[0], p[6]);
2855                     vop(p[4], p[10]); vop(p[4], p[6]); vop(p[2], p[8]); vop(p[2], p[4]); vop(p[6], p[8]);
2856                     vop(p[1], p[7]); vop(p[5], p[11]); vop(p[5], p[7]); vop(p[3], p[9]); vop(p[3], p[5]);
2857                     vop(p[7], p[9]); vop(p[1], p[2]); vop(p[3], p[4]); vop(p[5], p[6]); vop(p[7], p[8]);
2858                     vop(p[9], p[10]); vop(p[13], p[14]); vop(p[12], p[13]); vop(p[13], p[14]); vop(p[16], p[17]);
2859                     vop(p[15], p[16]); vop(p[16], p[17]); vop(p[12], p[15]); vop(p[14], p[17]); vop(p[14], p[15]);
2860                     vop(p[13], p[16]); vop(p[13], p[14]); vop(p[15], p[16]); vop(p[19], p[20]); vop(p[18], p[19]);
2861                     vop(p[19], p[20]); vop(p[21], p[22]); vop(p[23], p[24]); vop(p[21], p[23]); vop(p[22], p[24]);
2862                     vop(p[22], p[23]); vop(p[18], p[21]); vop(p[20], p[23]); vop(p[20], p[21]); vop(p[19], p[22]);
2863                     vop(p[22], p[24]); vop(p[19], p[20]); vop(p[21], p[22]); vop(p[23], p[24]); vop(p[12], p[18]);
2864                     vop(p[16], p[22]); vop(p[16], p[18]); vop(p[14], p[20]); vop(p[20], p[24]); vop(p[14], p[16]);
2865                     vop(p[18], p[20]); vop(p[22], p[24]); vop(p[13], p[19]); vop(p[17], p[23]); vop(p[17], p[19]);
2866                     vop(p[15], p[21]); vop(p[15], p[17]); vop(p[19], p[21]); vop(p[13], p[14]); vop(p[15], p[16]);
2867                     vop(p[17], p[18]); vop(p[19], p[20]); vop(p[21], p[22]); vop(p[23], p[24]); vop(p[0], p[12]);
2868                     vop(p[8], p[20]); vop(p[8], p[12]); vop(p[4], p[16]); vop(p[16], p[24]); vop(p[12], p[16]);
2869                     vop(p[2], p[14]); vop(p[10], p[22]); vop(p[10], p[14]); vop(p[6], p[18]); vop(p[6], p[10]);
2870                     vop(p[10], p[12]); vop(p[1], p[13]); vop(p[9], p[21]); vop(p[9], p[13]); vop(p[5], p[17]);
2871                     vop(p[13], p[17]); vop(p[3], p[15]); vop(p[11], p[23]); vop(p[11], p[15]); vop(p[7], p[19]);
2872                     vop(p[7], p[11]); vop(p[11], p[13]); vop(p[11], p[12]);
2873                     vop.store(dst+j, p[12]);
2874                 }
2875
2876                 limit = size.width;
2877             }
2878         }
2879     }
2880 }
2881
2882 #ifdef HAVE_OPENCL
2883
2884 static bool ocl_medianFilter(InputArray _src, OutputArray _dst, int m)
2885 {
2886     size_t localsize[2] = { 16, 16 };
2887     size_t globalsize[2];
2888     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2889
2890     if ( !((depth == CV_8U || depth == CV_16U || depth == CV_16S || depth == CV_32F) && cn <= 4 && (m == 3 || m == 5)) )
2891         return false;
2892
2893     Size imgSize = _src.size();
2894     bool useOptimized = (1 == cn) &&
2895                         (size_t)imgSize.width >= localsize[0] * 8  &&
2896                         (size_t)imgSize.height >= localsize[1] * 8 &&
2897                         imgSize.width % 4 == 0 &&
2898                         imgSize.height % 4 == 0 &&
2899                         (ocl::Device::getDefault().isIntel());
2900
2901     cv::String kname = format( useOptimized ? "medianFilter%d_u" : "medianFilter%d", m) ;
2902     cv::String kdefs = useOptimized ?
2903                          format("-D T=%s -D T1=%s -D T4=%s%d -D cn=%d -D USE_4OPT", ocl::typeToStr(type),
2904                          ocl::typeToStr(depth), ocl::typeToStr(depth), cn*4, cn)
2905                          :
2906                          format("-D T=%s -D T1=%s -D cn=%d", ocl::typeToStr(type), ocl::typeToStr(depth), cn) ;
2907
2908     ocl::Kernel k(kname.c_str(), ocl::imgproc::medianFilter_oclsrc, kdefs.c_str() );
2909
2910     if (k.empty())
2911         return false;
2912
2913     UMat src = _src.getUMat();
2914     _dst.create(src.size(), type);
2915     UMat dst = _dst.getUMat();
2916
2917     k.args(ocl::KernelArg::ReadOnlyNoSize(src), ocl::KernelArg::WriteOnly(dst));
2918
2919     if( useOptimized )
2920     {
2921         globalsize[0] = DIVUP(src.cols / 4, localsize[0]) * localsize[0];
2922         globalsize[1] = DIVUP(src.rows / 4, localsize[1]) * localsize[1];
2923     }
2924     else
2925     {
2926         globalsize[0] = (src.cols + localsize[0] + 2) / localsize[0] * localsize[0];
2927         globalsize[1] = (src.rows + localsize[1] - 1) / localsize[1] * localsize[1];
2928     }
2929
2930     return k.run(2, globalsize, localsize, false);
2931 }
2932
2933 #endif
2934
2935 }
2936
2937 #ifdef HAVE_OPENVX
2938 namespace cv
2939 {
2940     namespace ovx {
2941         template <> inline bool skipSmallImages<VX_KERNEL_MEDIAN_3x3>(int w, int h) { return w*h < 1280 * 720; }
2942     }
2943     static bool openvx_medianFilter(InputArray _src, OutputArray _dst, int ksize)
2944     {
2945         if (_src.type() != CV_8UC1 || _dst.type() != CV_8U
2946 #ifndef VX_VERSION_1_1
2947             || ksize != 3
2948 #endif
2949             )
2950             return false;
2951
2952         Mat src = _src.getMat();
2953         Mat dst = _dst.getMat();
2954
2955         if (
2956 #ifdef VX_VERSION_1_1
2957              ksize != 3 ? ovx::skipSmallImages<VX_KERNEL_NON_LINEAR_FILTER>(src.cols, src.rows) :
2958 #endif
2959              ovx::skipSmallImages<VX_KERNEL_MEDIAN_3x3>(src.cols, src.rows)
2960            )
2961             return false;
2962
2963         try
2964         {
2965             ivx::Context ctx = ovx::getOpenVXContext();
2966 #ifdef VX_VERSION_1_1
2967             if ((vx_size)ksize > ctx.nonlinearMaxDimension())
2968                 return false;
2969 #endif
2970
2971             Mat a;
2972             if (dst.data != src.data)
2973                 a = src;
2974             else
2975                 src.copyTo(a);
2976
2977             ivx::Image
2978                 ia = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
2979                     ivx::Image::createAddressing(a.cols, a.rows, 1, (vx_int32)(a.step)), a.data),
2980                 ib = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
2981                     ivx::Image::createAddressing(dst.cols, dst.rows, 1, (vx_int32)(dst.step)), dst.data);
2982
2983             //ATTENTION: VX_CONTEXT_IMMEDIATE_BORDER attribute change could lead to strange issues in multi-threaded environments
2984             //since OpenVX standart says nothing about thread-safety for now
2985             ivx::border_t prevBorder = ctx.immediateBorder();
2986             ctx.setImmediateBorder(VX_BORDER_REPLICATE);
2987 #ifdef VX_VERSION_1_1
2988             if (ksize == 3)
2989 #endif
2990             {
2991                 ivx::IVX_CHECK_STATUS(vxuMedian3x3(ctx, ia, ib));
2992             }
2993 #ifdef VX_VERSION_1_1
2994             else
2995             {
2996                 ivx::Matrix mtx;
2997                 if(ksize == 5)
2998                     mtx = ivx::Matrix::createFromPattern(ctx, VX_PATTERN_BOX, ksize, ksize);
2999                 else
3000                 {
3001                     vx_size supportedSize;
3002                     ivx::IVX_CHECK_STATUS(vxQueryContext(ctx, VX_CONTEXT_NONLINEAR_MAX_DIMENSION, &supportedSize, sizeof(supportedSize)));
3003                     if ((vx_size)ksize > supportedSize)
3004                     {
3005                         ctx.setImmediateBorder(prevBorder);
3006                         return false;
3007                     }
3008                     Mat mask(ksize, ksize, CV_8UC1, Scalar(255));
3009                     mtx = ivx::Matrix::create(ctx, VX_TYPE_UINT8, ksize, ksize);
3010                     mtx.copyFrom(mask);
3011                 }
3012                 ivx::IVX_CHECK_STATUS(vxuNonLinearFilter(ctx, VX_NONLINEAR_FILTER_MEDIAN, ia, mtx, ib));
3013             }
3014 #endif
3015             ctx.setImmediateBorder(prevBorder);
3016         }
3017         catch (ivx::RuntimeError & e)
3018         {
3019             VX_DbgThrow(e.what());
3020         }
3021         catch (ivx::WrapperError & e)
3022         {
3023             VX_DbgThrow(e.what());
3024         }
3025
3026         return true;
3027     }
3028 }
3029 #endif
3030
3031 #ifdef HAVE_IPP
3032 namespace cv
3033 {
3034 static bool ipp_medianFilter(Mat &src0, Mat &dst, int ksize)
3035 {
3036     CV_INSTRUMENT_REGION_IPP()
3037
3038 #if IPP_VERSION_X100 < 201801
3039     // Degradations for big kernel
3040     if(ksize > 7)
3041         return false;
3042 #endif
3043
3044     {
3045         int         bufSize;
3046         IppiSize    dstRoiSize = ippiSize(dst.cols, dst.rows), maskSize = ippiSize(ksize, ksize);
3047         IppDataType ippType = ippiGetDataType(src0.type());
3048         int         channels = src0.channels();
3049         IppAutoBuffer<Ipp8u> buffer;
3050
3051         if(src0.isSubmatrix())
3052             return false;
3053
3054         Mat src;
3055         if(dst.data != src0.data)
3056             src = src0;
3057         else
3058             src0.copyTo(src);
3059
3060         if(ippiFilterMedianBorderGetBufferSize(dstRoiSize, maskSize, ippType, channels, &bufSize) < 0)
3061             return false;
3062
3063         buffer.allocate(bufSize);
3064
3065         switch(ippType)
3066         {
3067         case ipp8u:
3068             if(channels == 1)
3069                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_8u_C1R, src.ptr<Ipp8u>(), (int)src.step, dst.ptr<Ipp8u>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3070             else if(channels == 3)
3071                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_8u_C3R, src.ptr<Ipp8u>(), (int)src.step, dst.ptr<Ipp8u>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3072             else if(channels == 4)
3073                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_8u_C4R, src.ptr<Ipp8u>(), (int)src.step, dst.ptr<Ipp8u>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3074             else
3075                 return false;
3076         case ipp16u:
3077             if(channels == 1)
3078                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_16u_C1R, src.ptr<Ipp16u>(), (int)src.step, dst.ptr<Ipp16u>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3079             else if(channels == 3)
3080                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_16u_C3R, src.ptr<Ipp16u>(), (int)src.step, dst.ptr<Ipp16u>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3081             else if(channels == 4)
3082                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_16u_C4R, src.ptr<Ipp16u>(), (int)src.step, dst.ptr<Ipp16u>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3083             else
3084                 return false;
3085         case ipp16s:
3086             if(channels == 1)
3087                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_16s_C1R, src.ptr<Ipp16s>(), (int)src.step, dst.ptr<Ipp16s>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3088             else if(channels == 3)
3089                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_16s_C3R, src.ptr<Ipp16s>(), (int)src.step, dst.ptr<Ipp16s>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3090             else if(channels == 4)
3091                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_16s_C4R, src.ptr<Ipp16s>(), (int)src.step, dst.ptr<Ipp16s>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3092             else
3093                 return false;
3094         case ipp32f:
3095             if(channels == 1)
3096                 return CV_INSTRUMENT_FUN_IPP(ippiFilterMedianBorder_32f_C1R, src.ptr<Ipp32f>(), (int)src.step, dst.ptr<Ipp32f>(), (int)dst.step, dstRoiSize, maskSize, ippBorderRepl, 0, buffer) >= 0;
3097             else
3098                 return false;
3099         default:
3100             return false;
3101         }
3102     }
3103 }
3104 }
3105 #endif
3106
3107 void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize )
3108 {
3109     CV_INSTRUMENT_REGION()
3110
3111     CV_Assert( (ksize % 2 == 1) && (_src0.dims() <= 2 ));
3112
3113     if( ksize <= 1 || _src0.empty() )
3114     {
3115         _src0.copyTo(_dst);
3116         return;
3117     }
3118
3119     CV_OCL_RUN(_dst.isUMat(),
3120                ocl_medianFilter(_src0,_dst, ksize))
3121
3122     Mat src0 = _src0.getMat();
3123     _dst.create( src0.size(), src0.type() );
3124     Mat dst = _dst.getMat();
3125
3126     CV_OVX_RUN(true,
3127                openvx_medianFilter(_src0, _dst, ksize))
3128
3129     CV_IPP_RUN_FAST(ipp_medianFilter(src0, dst, ksize));
3130
3131 #ifdef HAVE_TEGRA_OPTIMIZATION
3132     if (tegra::useTegra() && tegra::medianBlur(src0, dst, ksize))
3133         return;
3134 #endif
3135
3136     bool useSortNet = ksize == 3 || (ksize == 5
3137 #if !(CV_SIMD128)
3138             && ( src0.depth() > CV_8U || src0.channels() == 2 || src0.channels() > 4 )
3139 #endif
3140         );
3141
3142     Mat src;
3143     if( useSortNet )
3144     {
3145         if( dst.data != src0.data )
3146             src = src0;
3147         else
3148             src0.copyTo(src);
3149
3150         if( src.depth() == CV_8U )
3151             medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize );
3152         else if( src.depth() == CV_16U )
3153             medianBlur_SortNet<MinMax16u, MinMaxVec16u>( src, dst, ksize );
3154         else if( src.depth() == CV_16S )
3155             medianBlur_SortNet<MinMax16s, MinMaxVec16s>( src, dst, ksize );
3156         else if( src.depth() == CV_32F )
3157             medianBlur_SortNet<MinMax32f, MinMaxVec32f>( src, dst, ksize );
3158         else
3159             CV_Error(CV_StsUnsupportedFormat, "");
3160
3161         return;
3162     }
3163     else
3164     {
3165         cv::copyMakeBorder( src0, src, 0, 0, ksize/2, ksize/2, BORDER_REPLICATE );
3166
3167         int cn = src0.channels();
3168         CV_Assert( src.depth() == CV_8U && (cn == 1 || cn == 3 || cn == 4) );
3169
3170         double img_size_mp = (double)(src0.total())/(1 << 20);
3171         if( ksize <= 3 + (img_size_mp < 1 ? 12 : img_size_mp < 4 ? 6 : 2)*
3172             (CV_SIMD128 && hasSIMD128() ? 1 : 3))
3173             medianBlur_8u_Om( src, dst, ksize );
3174         else
3175             medianBlur_8u_O1( src, dst, ksize );
3176     }
3177 }
3178
3179 /****************************************************************************************\
3180                                    Bilateral Filtering
3181 \****************************************************************************************/
3182
3183 namespace cv
3184 {
3185
3186 class BilateralFilter_8u_Invoker :
3187     public ParallelLoopBody
3188 {
3189 public:
3190     BilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, int _radius, int _maxk,
3191         int* _space_ofs, float *_space_weight, float *_color_weight) :
3192         temp(&_temp), dest(&_dest), radius(_radius),
3193         maxk(_maxk), space_ofs(_space_ofs), space_weight(_space_weight), color_weight(_color_weight)
3194     {
3195     }
3196
3197     virtual void operator() (const Range& range) const
3198     {
3199         int i, j, cn = dest->channels(), k;
3200         Size size = dest->size();
3201 #if CV_SIMD128
3202         int CV_DECL_ALIGNED(16) buf[4];
3203         bool haveSIMD128 = hasSIMD128();
3204 #endif
3205
3206         for( i = range.start; i < range.end; i++ )
3207         {
3208             const uchar* sptr = temp->ptr(i+radius) + radius*cn;
3209             uchar* dptr = dest->ptr(i);
3210
3211             if( cn == 1 )
3212             {
3213                 for( j = 0; j < size.width; j++ )
3214                 {
3215                     float sum = 0, wsum = 0;
3216                     int val0 = sptr[j];
3217                     k = 0;
3218 #if CV_SIMD128
3219                     if( haveSIMD128 )
3220                     {
3221                         v_float32x4 _val0 = v_setall_f32(static_cast<float>(val0));
3222                         v_float32x4 vsumw = v_setzero_f32();
3223                         v_float32x4 vsumc = v_setzero_f32();
3224
3225                         for( ; k <= maxk - 4; k += 4 )
3226                         {
3227                             v_float32x4 _valF = v_float32x4(sptr[j + space_ofs[k]],
3228                                 sptr[j + space_ofs[k + 1]],
3229                                 sptr[j + space_ofs[k + 2]],
3230                                 sptr[j + space_ofs[k + 3]]);
3231                             v_float32x4 _val = v_abs(_valF - _val0);
3232                             v_store(buf, v_round(_val));
3233
3234                             v_float32x4 _cw = v_float32x4(color_weight[buf[0]],
3235                                 color_weight[buf[1]],
3236                                 color_weight[buf[2]],
3237                                 color_weight[buf[3]]);
3238                             v_float32x4 _sw = v_load(space_weight+k);
3239                             v_float32x4 _w = _cw * _sw;
3240                             _cw = _w * _valF;
3241
3242                             vsumw += _w;
3243                             vsumc += _cw;
3244                         }
3245                         float *bufFloat = (float*)buf;
3246                         v_float32x4 sum4 = v_reduce_sum4(vsumw, vsumc, vsumw, vsumc);
3247                         v_store(bufFloat, sum4);
3248                         sum += bufFloat[1];
3249                         wsum += bufFloat[0];
3250                     }
3251 #endif
3252                     for( ; k < maxk; k++ )
3253                     {
3254                         int val = sptr[j + space_ofs[k]];
3255                         float w = space_weight[k]*color_weight[std::abs(val - val0)];
3256                         sum += val*w;
3257                         wsum += w;
3258                     }
3259                     // overflow is not possible here => there is no need to use cv::saturate_cast
3260                     dptr[j] = (uchar)cvRound(sum/wsum);
3261                 }
3262             }
3263             else
3264             {
3265                 assert( cn == 3 );
3266                 for( j = 0; j < size.width*3; j += 3 )
3267                 {
3268                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
3269                     int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
3270                     k = 0;
3271 #if CV_SIMD128
3272                     if( haveSIMD128 )
3273                     {
3274                         v_float32x4 vsumw = v_setzero_f32();
3275                         v_float32x4 vsumb = v_setzero_f32();
3276                         v_float32x4 vsumg = v_setzero_f32();
3277                         v_float32x4 vsumr = v_setzero_f32();
3278                         const v_float32x4 _b0 = v_setall_f32(static_cast<float>(b0));
3279                         const v_float32x4 _g0 = v_setall_f32(static_cast<float>(g0));
3280                         const v_float32x4 _r0 = v_setall_f32(static_cast<float>(r0));
3281
3282                         for( ; k <= maxk - 4; k += 4 )
3283                         {
3284                             const uchar* const sptr_k0  = sptr + j + space_ofs[k];
3285                             const uchar* const sptr_k1  = sptr + j + space_ofs[k+1];
3286                             const uchar* const sptr_k2  = sptr + j + space_ofs[k+2];
3287                             const uchar* const sptr_k3  = sptr + j + space_ofs[k+3];
3288
3289                             v_float32x4 __b = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k0)));
3290                             v_float32x4 __g = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k1)));
3291                             v_float32x4 __r = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k2)));
3292                             v_float32x4 __z = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k3)));
3293                             v_float32x4 _b, _g, _r, _z;
3294
3295                             v_transpose4x4(__b, __g, __r, __z, _b, _g, _r, _z);
3296
3297                             v_float32x4 bt = v_abs(_b -_b0);
3298                             v_float32x4 gt = v_abs(_g -_g0);
3299                             v_float32x4 rt = v_abs(_r -_r0);
3300
3301                             bt = rt + bt + gt;
3302                             v_store(buf, v_round(bt));
3303
3304                             v_float32x4 _w  = v_float32x4(color_weight[buf[0]],color_weight[buf[1]],
3305                                                     color_weight[buf[2]],color_weight[buf[3]]);
3306                             v_float32x4 _sw = v_load(space_weight+k);
3307
3308                             _w *= _sw;
3309                             _b *=  _w;
3310                             _g *=  _w;
3311                             _r *=  _w;
3312
3313                             vsumw += _w;
3314                             vsumb += _b;
3315                             vsumg += _g;
3316                             vsumr += _r;
3317                         }
3318                         float *bufFloat = (float*)buf;
3319                         v_float32x4 sum4 = v_reduce_sum4(vsumw, vsumb, vsumg, vsumr);
3320                         v_store(bufFloat, sum4);
3321                         wsum += bufFloat[0];
3322                         sum_b += bufFloat[1];
3323                         sum_g += bufFloat[2];
3324                         sum_r += bufFloat[3];
3325                     }
3326 #endif
3327
3328                     for( ; k < maxk; k++ )
3329                     {
3330                         const uchar* sptr_k = sptr + j + space_ofs[k];
3331                         int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
3332                         float w = space_weight[k]*color_weight[std::abs(b - b0) +
3333                                                                std::abs(g - g0) + std::abs(r - r0)];
3334                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
3335                         wsum += w;
3336                     }
3337                     wsum = 1.f/wsum;
3338                     b0 = cvRound(sum_b*wsum);
3339                     g0 = cvRound(sum_g*wsum);
3340                     r0 = cvRound(sum_r*wsum);
3341                     dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
3342                 }
3343             }
3344         }
3345     }
3346
3347 private:
3348     const Mat *temp;
3349     Mat *dest;
3350     int radius, maxk, *space_ofs;
3351     float *space_weight, *color_weight;
3352 };
3353
3354 #ifdef HAVE_OPENCL
3355
3356 static bool ocl_bilateralFilter_8u(InputArray _src, OutputArray _dst, int d,
3357                                    double sigma_color, double sigma_space,
3358                                    int borderType)
3359 {
3360 #ifdef __ANDROID__
3361     if (ocl::Device::getDefault().isNVidia())
3362         return false;
3363 #endif
3364
3365     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
3366     int i, j, maxk, radius;
3367
3368     if (depth != CV_8U || cn > 4)
3369         return false;
3370
3371     if (sigma_color <= 0)
3372         sigma_color = 1;
3373     if (sigma_space <= 0)
3374         sigma_space = 1;
3375
3376     double gauss_color_coeff = -0.5 / (sigma_color * sigma_color);
3377     double gauss_space_coeff = -0.5 / (sigma_space * sigma_space);
3378
3379     if ( d <= 0 )
3380         radius = cvRound(sigma_space * 1.5);
3381     else
3382         radius = d / 2;
3383     radius = MAX(radius, 1);
3384     d = radius * 2 + 1;
3385
3386     UMat src = _src.getUMat(), dst = _dst.getUMat(), temp;
3387     if (src.u == dst.u)
3388         return false;
3389
3390     copyMakeBorder(src, temp, radius, radius, radius, radius, borderType);
3391     std::vector<float> _space_weight(d * d);
3392     std::vector<int> _space_ofs(d * d);
3393     float * const space_weight = &_space_weight[0];
3394     int * const space_ofs = &_space_ofs[0];
3395
3396     // initialize space-related bilateral filter coefficients
3397     for( i = -radius, maxk = 0; i <= radius; i++ )
3398         for( j = -radius; j <= radius; j++ )
3399         {
3400             double r = std::sqrt((double)i * i + (double)j * j);
3401             if ( r > radius )
3402                 continue;
3403             space_weight[maxk] = (float)std::exp(r * r * gauss_space_coeff);
3404             space_ofs[maxk++] = (int)(i * temp.step + j * cn);
3405         }
3406
3407     char cvt[3][40];
3408     String cnstr = cn > 1 ? format("%d", cn) : "";
3409     String kernelName("bilateral");
3410     size_t sizeDiv = 1;
3411     if ((ocl::Device::getDefault().isIntel()) &&
3412         (ocl::Device::getDefault().type() == ocl::Device::TYPE_GPU))
3413     {
3414             //Intel GPU
3415             if (dst.cols % 4 == 0 && cn == 1) // For single channel x4 sized images.
3416             {
3417                 kernelName = "bilateral_float4";
3418                 sizeDiv = 4;
3419             }
3420      }
3421      ocl::Kernel k(kernelName.c_str(), ocl::imgproc::bilateral_oclsrc,
3422             format("-D radius=%d -D maxk=%d -D cn=%d -D int_t=%s -D uint_t=uint%s -D convert_int_t=%s"
3423             " -D uchar_t=%s -D float_t=%s -D convert_float_t=%s -D convert_uchar_t=%s -D gauss_color_coeff=(float)%f",
3424             radius, maxk, cn, ocl::typeToStr(CV_32SC(cn)), cnstr.c_str(),
3425             ocl::convertTypeStr(CV_8U, CV_32S, cn, cvt[0]),
3426             ocl::typeToStr(type), ocl::typeToStr(CV_32FC(cn)),
3427             ocl::convertTypeStr(CV_32S, CV_32F, cn, cvt[1]),
3428             ocl::convertTypeStr(CV_32F, CV_8U, cn, cvt[2]), gauss_color_coeff));
3429     if (k.empty())
3430         return false;
3431
3432     Mat mspace_weight(1, d * d, CV_32FC1, space_weight);
3433     Mat mspace_ofs(1, d * d, CV_32SC1, space_ofs);
3434     UMat ucolor_weight, uspace_weight, uspace_ofs;
3435
3436     mspace_weight.copyTo(uspace_weight);
3437     mspace_ofs.copyTo(uspace_ofs);
3438
3439     k.args(ocl::KernelArg::ReadOnlyNoSize(temp), ocl::KernelArg::WriteOnly(dst),
3440            ocl::KernelArg::PtrReadOnly(uspace_weight),
3441            ocl::KernelArg::PtrReadOnly(uspace_ofs));
3442
3443     size_t globalsize[2] = { (size_t)dst.cols / sizeDiv, (size_t)dst.rows };
3444     return k.run(2, globalsize, NULL, false);
3445 }
3446
3447 #endif
3448 static void
3449 bilateralFilter_8u( const Mat& src, Mat& dst, int d,
3450     double sigma_color, double sigma_space,
3451     int borderType )
3452 {
3453     int cn = src.channels();
3454     int i, j, maxk, radius;
3455     Size size = src.size();
3456
3457     CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) && src.data != dst.data );
3458
3459     if( sigma_color <= 0 )
3460         sigma_color = 1;
3461     if( sigma_space <= 0 )
3462         sigma_space = 1;
3463
3464     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
3465     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
3466
3467     if( d <= 0 )
3468         radius = cvRound(sigma_space*1.5);
3469     else
3470         radius = d/2;
3471     radius = MAX(radius, 1);
3472     d = radius*2 + 1;
3473
3474     Mat temp;
3475     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
3476
3477     std::vector<float> _color_weight(cn*256);
3478     std::vector<float> _space_weight(d*d);
3479     std::vector<int> _space_ofs(d*d);
3480     float* color_weight = &_color_weight[0];
3481     float* space_weight = &_space_weight[0];
3482     int* space_ofs = &_space_ofs[0];
3483
3484     // initialize color-related bilateral filter coefficients
3485
3486     for( i = 0; i < 256*cn; i++ )
3487         color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);
3488
3489     // initialize space-related bilateral filter coefficients
3490     for( i = -radius, maxk = 0; i <= radius; i++ )
3491     {
3492         j = -radius;
3493
3494         for( ; j <= radius; j++ )
3495         {
3496             double r = std::sqrt((double)i*i + (double)j*j);
3497             if( r > radius )
3498                 continue;
3499             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
3500             space_ofs[maxk++] = (int)(i*temp.step + j*cn);
3501         }
3502     }
3503
3504     BilateralFilter_8u_Invoker body(dst, temp, radius, maxk, space_ofs, space_weight, color_weight);
3505     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
3506 }
3507
3508
3509 class BilateralFilter_32f_Invoker :
3510     public ParallelLoopBody
3511 {
3512 public:
3513
3514     BilateralFilter_32f_Invoker(int _cn, int _radius, int _maxk, int *_space_ofs,
3515         const Mat& _temp, Mat& _dest, float _scale_index, float *_space_weight, float *_expLUT) :
3516         cn(_cn), radius(_radius), maxk(_maxk), space_ofs(_space_ofs),
3517         temp(&_temp), dest(&_dest), scale_index(_scale_index), space_weight(_space_weight), expLUT(_expLUT)
3518     {
3519     }
3520
3521     virtual void operator() (const Range& range) const
3522     {
3523         int i, j, k;
3524         Size size = dest->size();
3525 #if CV_SIMD128
3526         int CV_DECL_ALIGNED(16) idxBuf[4];
3527         bool haveSIMD128 = hasSIMD128();
3528 #endif
3529
3530         for( i = range.start; i < range.end; i++ )
3531         {
3532             const float* sptr = temp->ptr<float>(i+radius) + radius*cn;
3533             float* dptr = dest->ptr<float>(i);
3534
3535             if( cn == 1 )
3536             {
3537                 for( j = 0; j < size.width; j++ )
3538                 {
3539                     float sum = 0, wsum = 0;
3540                     float val0 = sptr[j];
3541                     k = 0;
3542 #if CV_SIMD128
3543                     if( haveSIMD128 )
3544                     {
3545                         v_float32x4 vecwsum = v_setzero_f32();
3546                         v_float32x4 vecvsum = v_setzero_f32();
3547                         const v_float32x4 _val0 = v_setall_f32(sptr[j]);
3548                         const v_float32x4 _scale_index = v_setall_f32(scale_index);
3549
3550                         for (; k <= maxk - 4; k += 4)
3551                         {
3552                             v_float32x4 _sw = v_load(space_weight + k);
3553                             v_float32x4 _val = v_float32x4(sptr[j + space_ofs[k]],
3554                                 sptr[j + space_ofs[k + 1]],
3555                                 sptr[j + space_ofs[k + 2]],
3556                                 sptr[j + space_ofs[k + 3]]);
3557                             v_float32x4 _alpha = v_abs(_val - _val0) * _scale_index;
3558
3559                             v_int32x4 _idx = v_round(_alpha);
3560                             v_store(idxBuf, _idx);
3561                             _alpha -= v_cvt_f32(_idx);
3562
3563                             v_float32x4 _explut = v_float32x4(expLUT[idxBuf[0]],
3564                                 expLUT[idxBuf[1]],
3565                                 expLUT[idxBuf[2]],
3566                                 expLUT[idxBuf[3]]);
3567                             v_float32x4 _explut1 = v_float32x4(expLUT[idxBuf[0] + 1],
3568                                 expLUT[idxBuf[1] + 1],
3569                                 expLUT[idxBuf[2] + 1],
3570                                 expLUT[idxBuf[3] + 1]);
3571
3572                             v_float32x4 _w = _sw * (_explut + (_alpha * (_explut1 - _explut)));
3573                             _val *= _w;
3574
3575                             vecwsum += _w;
3576                             vecvsum += _val;
3577                         }
3578                         float *bufFloat = (float*)idxBuf;
3579                         v_float32x4 sum4 = v_reduce_sum4(vecwsum, vecvsum, vecwsum, vecvsum);
3580                         v_store(bufFloat, sum4);
3581                         sum += bufFloat[1];
3582                         wsum += bufFloat[0];
3583                     }
3584 #endif
3585
3586                     for( ; k < maxk; k++ )
3587                     {
3588                         float val = sptr[j + space_ofs[k]];
3589                         float alpha = (float)(std::abs(val - val0)*scale_index);
3590                         int idx = cvFloor(alpha);
3591                         alpha -= idx;
3592                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
3593                         sum += val*w;
3594                         wsum += w;
3595                     }
3596                     dptr[j] = (float)(sum/wsum);
3597                 }
3598             }
3599             else
3600             {
3601                 CV_Assert( cn == 3 );
3602                 for( j = 0; j < size.width*3; j += 3 )
3603                 {
3604                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
3605                     float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
3606                     k = 0;
3607 #if CV_SIMD128
3608                     if( haveSIMD128 )
3609                     {
3610                         v_float32x4 sumw = v_setzero_f32();
3611                         v_float32x4 sumb = v_setzero_f32();
3612                         v_float32x4 sumg = v_setzero_f32();
3613                         v_float32x4 sumr = v_setzero_f32();
3614                         const v_float32x4 _b0 = v_setall_f32(b0);
3615                         const v_float32x4 _g0 = v_setall_f32(g0);
3616                         const v_float32x4 _r0 = v_setall_f32(r0);
3617                         const v_float32x4 _scale_index = v_setall_f32(scale_index);
3618
3619                         for( ; k <= maxk-4; k += 4 )
3620                         {
3621                             v_float32x4 _sw = v_load(space_weight + k);
3622
3623                             const float* const sptr_k0 = sptr + j + space_ofs[k];
3624                             const float* const sptr_k1 = sptr + j + space_ofs[k+1];
3625                             const float* const sptr_k2 = sptr + j + space_ofs[k+2];
3626                             const float* const sptr_k3 = sptr + j + space_ofs[k+3];
3627
3628                             v_float32x4 _v0 = v_load(sptr_k0);
3629                             v_float32x4 _v1 = v_load(sptr_k1);
3630                             v_float32x4 _v2 = v_load(sptr_k2);
3631                             v_float32x4 _v3 = v_load(sptr_k3);
3632                             v_float32x4 _b, _g, _r, _dummy;
3633
3634                             v_transpose4x4(_v0, _v1, _v2, _v3, _b, _g, _r, _dummy);
3635
3636                             v_float32x4 _bt = v_abs(_b - _b0);
3637                             v_float32x4 _gt = v_abs(_g - _g0);
3638                             v_float32x4 _rt = v_abs(_r - _r0);
3639                             v_float32x4 _alpha = _scale_index * (_bt + _gt + _rt);
3640
3641                             v_int32x4 _idx = v_round(_alpha);
3642                             v_store((int*)idxBuf, _idx);
3643                             _alpha -= v_cvt_f32(_idx);
3644
3645                             v_float32x4 _explut = v_float32x4(expLUT[idxBuf[0]],
3646                                 expLUT[idxBuf[1]],
3647                                 expLUT[idxBuf[2]],
3648                                 expLUT[idxBuf[3]]);
3649                             v_float32x4 _explut1 = v_float32x4(expLUT[idxBuf[0] + 1],
3650                                 expLUT[idxBuf[1] + 1],
3651                                 expLUT[idxBuf[2] + 1],
3652                                 expLUT[idxBuf[3] + 1]);
3653
3654                             v_float32x4 _w = _sw * (_explut + (_alpha * (_explut1 - _explut)));
3655
3656                             _b *=  _w;
3657                             _g *=  _w;
3658                             _r *=  _w;
3659                             sumw += _w;
3660                             sumb += _b;
3661                             sumg += _g;
3662                             sumr += _r;
3663                         }
3664                         v_float32x4 sum4 = v_reduce_sum4(sumw, sumb, sumg, sumr);
3665                         float *bufFloat = (float*)idxBuf;
3666                         v_store(bufFloat, sum4);
3667                         wsum += bufFloat[0];
3668                         sum_b += bufFloat[1];
3669                         sum_g += bufFloat[2];
3670                         sum_r += bufFloat[3];
3671                     }
3672 #endif
3673
3674                     for(; k < maxk; k++ )
3675                     {
3676                         const float* sptr_k = sptr + j + space_ofs[k];
3677                         float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
3678                         float alpha = (float)((std::abs(b - b0) +
3679                             std::abs(g - g0) + std::abs(r - r0))*scale_index);
3680                         int idx = cvFloor(alpha);
3681                         alpha -= idx;
3682                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
3683                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
3684                         wsum += w;
3685                     }
3686                     wsum = 1.f/wsum;
3687                     b0 = sum_b*wsum;
3688                     g0 = sum_g*wsum;
3689                     r0 = sum_r*wsum;
3690                     dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
3691                 }
3692             }
3693         }
3694     }
3695
3696 private:
3697     int cn, radius, maxk, *space_ofs;
3698     const Mat* temp;
3699     Mat *dest;
3700     float scale_index, *space_weight, *expLUT;
3701 };
3702
3703
3704 static void
3705 bilateralFilter_32f( const Mat& src, Mat& dst, int d,
3706                      double sigma_color, double sigma_space,
3707                      int borderType )
3708 {
3709     int cn = src.channels();
3710     int i, j, maxk, radius;
3711     double minValSrc=-1, maxValSrc=1;
3712     const int kExpNumBinsPerChannel = 1 << 12;
3713     int kExpNumBins = 0;
3714     float lastExpVal = 1.f;
3715     float len, scale_index;
3716     Size size = src.size();
3717
3718     CV_Assert( (src.type() == CV_32FC1 || src.type() == CV_32FC3) && src.data != dst.data );
3719
3720     if( sigma_color <= 0 )
3721         sigma_color = 1;
3722     if( sigma_space <= 0 )
3723         sigma_space = 1;
3724
3725     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
3726     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
3727
3728     if( d <= 0 )
3729         radius = cvRound(sigma_space*1.5);
3730     else
3731         radius = d/2;
3732     radius = MAX(radius, 1);
3733     d = radius*2 + 1;
3734     // compute the min/max range for the input image (even if multichannel)
3735
3736     minMaxLoc( src.reshape(1), &minValSrc, &maxValSrc );
3737     if(std::abs(minValSrc - maxValSrc) < FLT_EPSILON)
3738     {
3739         src.copyTo(dst);
3740         return;
3741     }
3742
3743     // temporary copy of the image with borders for easy processing
3744     Mat temp;
3745     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
3746     const double insteadNaNValue = -5. * sigma_color;
3747     patchNaNs( temp, insteadNaNValue ); // this replacement of NaNs makes the assumption that depth values are nonnegative
3748                                         // TODO: make insteadNaNValue avalible in the outside function interface to control the cases breaking the assumption
3749     // allocate lookup tables
3750     std::vector<float> _space_weight(d*d);
3751     std::vector<int> _space_ofs(d*d);
3752     float* space_weight = &_space_weight[0];
3753     int* space_ofs = &_space_ofs[0];
3754
3755     // assign a length which is slightly more than needed
3756     len = (float)(maxValSrc - minValSrc) * cn;
3757     kExpNumBins = kExpNumBinsPerChannel * cn;
3758     std::vector<float> _expLUT(kExpNumBins+2);
3759     float* expLUT = &_expLUT[0];
3760
3761     scale_index = kExpNumBins/len;
3762
3763     // initialize the exp LUT
3764     for( i = 0; i < kExpNumBins+2; i++ )
3765     {
3766         if( lastExpVal > 0.f )
3767         {
3768             double val =  i / scale_index;
3769             expLUT[i] = (float)std::exp(val * val * gauss_color_coeff);
3770             lastExpVal = expLUT[i];
3771         }
3772         else
3773             expLUT[i] = 0.f;
3774     }
3775
3776     // initialize space-related bilateral filter coefficients
3777     for( i = -radius, maxk = 0; i <= radius; i++ )
3778         for( j = -radius; j <= radius; j++ )
3779         {
3780             double r = std::sqrt((double)i*i + (double)j*j);
3781             if( r > radius )
3782                 continue;
3783             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
3784             space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
3785         }
3786
3787     // parallel_for usage
3788
3789     BilateralFilter_32f_Invoker body(cn, radius, maxk, space_ofs, temp, dst, scale_index, space_weight, expLUT);
3790     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
3791 }
3792
3793 #ifdef HAVE_IPP
3794 #define IPP_BILATERAL_PARALLEL 1
3795
3796 #ifdef HAVE_IPP_IW
3797 class ipp_bilateralFilterParallel: public ParallelLoopBody
3798 {
3799 public:
3800     ipp_bilateralFilterParallel(::ipp::IwiImage &_src, ::ipp::IwiImage &_dst, int _radius, Ipp32f _valSquareSigma, Ipp32f _posSquareSigma, ::ipp::IwiBorderType _borderType, bool *_ok):
3801         src(_src), dst(_dst)
3802     {
3803         pOk = _ok;
3804
3805         radius          = _radius;
3806         valSquareSigma  = _valSquareSigma;
3807         posSquareSigma  = _posSquareSigma;
3808         borderType      = _borderType;
3809
3810         *pOk = true;
3811     }
3812     ~ipp_bilateralFilterParallel() {}
3813
3814     virtual void operator() (const Range& range) const
3815     {
3816         if(*pOk == false)
3817             return;
3818
3819         try
3820         {
3821             ::ipp::IwiTile tile = ::ipp::IwiRoi(0, range.start, dst.m_size.width, range.end - range.start);
3822             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBilateral, src, dst, radius, valSquareSigma, posSquareSigma, ::ipp::IwDefault(), borderType, tile);
3823         }
3824         catch(::ipp::IwException)
3825         {
3826             *pOk = false;
3827             return;
3828         }
3829     }
3830 private:
3831     ::ipp::IwiImage &src;
3832     ::ipp::IwiImage &dst;
3833
3834     int                  radius;
3835     Ipp32f               valSquareSigma;
3836     Ipp32f               posSquareSigma;
3837     ::ipp::IwiBorderType borderType;
3838
3839     bool  *pOk;
3840     const ipp_bilateralFilterParallel& operator= (const ipp_bilateralFilterParallel&);
3841 };
3842 #endif
3843
3844 static bool ipp_bilateralFilter(Mat &src, Mat &dst, int d, double sigmaColor, double sigmaSpace, int borderType)
3845 {
3846 #ifdef HAVE_IPP_IW
3847     CV_INSTRUMENT_REGION_IPP()
3848
3849     int         radius         = IPP_MAX(((d <= 0)?cvRound(sigmaSpace*1.5):d/2), 1);
3850     Ipp32f      valSquareSigma = (Ipp32f)((sigmaColor <= 0)?1:sigmaColor*sigmaColor);
3851     Ipp32f      posSquareSigma = (Ipp32f)((sigmaSpace <= 0)?1:sigmaSpace*sigmaSpace);
3852
3853     // Acquire data and begin processing
3854     try
3855     {
3856         ::ipp::IwiImage      iwSrc = ippiGetImage(src);
3857         ::ipp::IwiImage      iwDst = ippiGetImage(dst);
3858         ::ipp::IwiBorderSize borderSize(radius);
3859         ::ipp::IwiBorderType ippBorder(ippiGetBorder(iwSrc, borderType, borderSize));
3860         if(!ippBorder)
3861             return false;
3862
3863         const int threads = ippiSuggestThreadsNum(iwDst, 2);
3864         if(IPP_BILATERAL_PARALLEL && threads > 1) {
3865             bool  ok      = true;
3866             Range range(0, (int)iwDst.m_size.height);
3867             ipp_bilateralFilterParallel invoker(iwSrc, iwDst, radius, valSquareSigma, posSquareSigma, ippBorder, &ok);
3868             if(!ok)
3869                 return false;
3870
3871             parallel_for_(range, invoker, threads*4);
3872
3873             if(!ok)
3874                 return false;
3875         } else {
3876             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBilateral, iwSrc, iwDst, radius, valSquareSigma, posSquareSigma, ::ipp::IwDefault(), ippBorder);
3877         }
3878     }
3879     catch (::ipp::IwException)
3880     {
3881         return false;
3882     }
3883     return true;
3884 #else
3885     CV_UNUSED(src); CV_UNUSED(dst); CV_UNUSED(d); CV_UNUSED(sigmaColor); CV_UNUSED(sigmaSpace); CV_UNUSED(borderType);
3886     return false;
3887 #endif
3888 }
3889 #endif
3890
3891 }
3892
3893 void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,
3894                       double sigmaColor, double sigmaSpace,
3895                       int borderType )
3896 {
3897     CV_INSTRUMENT_REGION()
3898
3899     _dst.create( _src.size(), _src.type() );
3900
3901     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
3902                ocl_bilateralFilter_8u(_src, _dst, d, sigmaColor, sigmaSpace, borderType))
3903
3904     Mat src = _src.getMat(), dst = _dst.getMat();
3905
3906     CV_IPP_RUN_FAST(ipp_bilateralFilter(src, dst, d, sigmaColor, sigmaSpace, borderType));
3907
3908     if( src.depth() == CV_8U )
3909         bilateralFilter_8u( src, dst, d, sigmaColor, sigmaSpace, borderType );
3910     else if( src.depth() == CV_32F )
3911         bilateralFilter_32f( src, dst, d, sigmaColor, sigmaSpace, borderType );
3912     else
3913         CV_Error( CV_StsUnsupportedFormat,
3914         "Bilateral filtering is only implemented for 8u and 32f images" );
3915 }
3916
3917 //////////////////////////////////////////////////////////////////////////////////////////
3918
3919 CV_IMPL void
3920 cvSmooth( const void* srcarr, void* dstarr, int smooth_type,
3921           int param1, int param2, double param3, double param4 )
3922 {
3923     cv::Mat src = cv::cvarrToMat(srcarr), dst0 = cv::cvarrToMat(dstarr), dst = dst0;
3924
3925     CV_Assert( dst.size() == src.size() &&
3926         (smooth_type == CV_BLUR_NO_SCALE || dst.type() == src.type()) );
3927
3928     if( param2 <= 0 )
3929         param2 = param1;
3930
3931     if( smooth_type == CV_BLUR || smooth_type == CV_BLUR_NO_SCALE )
3932         cv::boxFilter( src, dst, dst.depth(), cv::Size(param1, param2), cv::Point(-1,-1),
3933             smooth_type == CV_BLUR, cv::BORDER_REPLICATE );
3934     else if( smooth_type == CV_GAUSSIAN )
3935         cv::GaussianBlur( src, dst, cv::Size(param1, param2), param3, param4, cv::BORDER_REPLICATE );
3936     else if( smooth_type == CV_MEDIAN )
3937         cv::medianBlur( src, dst, param1 );
3938     else
3939         cv::bilateralFilter( src, dst, param1, param3, param4, cv::BORDER_REPLICATE );
3940
3941     if( dst.data != dst0.data )
3942         CV_Error( CV_StsUnmatchedFormats, "The destination image does not have the proper type" );
3943 }
3944
3945 /* End of file. */