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
  46 #include <vector>
  47
  48 #include "opencv2/core/hal/intrin.hpp"
  49 #include "opencl_kernels_imgproc.hpp"
  50
  51 #include "opencv2/core/openvx/ovx_defs.hpp"
  52
  53 #include "filter.hpp"
  54
  55 #include "fixedpoint.inl.hpp"
  56 /*
  57  * This file includes the code, contributed by Simon Perreault
  58  * (the function icvMedianBlur_8u_O1)
  59  *
  60  * Constant-time median filtering -- http://nomis80.org/ctmf.html
  61  * Copyright (C) 2006 Simon Perreault
  62  *
  63  * Contact:
  64  *  Laboratoire de vision et systemes numeriques
  65  *  Pavillon Adrien-Pouliot
  66  *  Universite Laval
  67  *  Sainte-Foy, Quebec, Canada
  68  *  G1K 7P4
  69  *
  70  *  perreaul@gel.ulaval.ca
  71  */
  72
  73 namespace cv
  74 {
  75
  76 /****************************************************************************************\
  77                                          Box Filter
  78 \****************************************************************************************/
  79
  80 template<typename T, typename ST>
  81 struct RowSum :
  82         public BaseRowFilter
  83 {
  84     RowSum( int _ksize, int _anchor ) :
  85         BaseRowFilter()
  86     {
  87         ksize = _ksize;
  88         anchor = _anchor;
  89     }
  90
  91     virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
  92     {
  93         const T* S = (const T*)src;
  94         ST* D = (ST*)dst;
  95         int i = 0, k, ksz_cn = ksize*cn;
  96
  97         width = (width - 1)*cn;
  98         if( ksize == 3 )
  99         {
 100             for( i = 0; i < width + cn; i++ )
 101             {
 102                 D[i] = (ST)S[i] + (ST)S[i+cn] + (ST)S[i+cn*2];
 103             }
 104         }
 105         else if( ksize == 5 )
 106         {
 107             for( i = 0; i < width + cn; i++ )
 108             {
 109                 D[i] = (ST)S[i] + (ST)S[i+cn] + (ST)S[i+cn*2] + (ST)S[i + cn*3] + (ST)S[i + cn*4];
 110             }
 111         }
 112         else if( cn == 1 )
 113         {
 114             ST s = 0;
 115             for( i = 0; i < ksz_cn; i++ )
 116                 s += (ST)S[i];
 117             D[0] = s;
 118             for( i = 0; i < width; i++ )
 119             {
 120                 s += (ST)S[i + ksz_cn] - (ST)S[i];
 121                 D[i+1] = s;
 122             }
 123         }
 124         else if( cn == 3 )
 125         {
 126             ST s0 = 0, s1 = 0, s2 = 0;
 127             for( i = 0; i < ksz_cn; i += 3 )
 128             {
 129                 s0 += (ST)S[i];
 130                 s1 += (ST)S[i+1];
 131                 s2 += (ST)S[i+2];
 132             }
 133             D[0] = s0;
 134             D[1] = s1;
 135             D[2] = s2;
 136             for( i = 0; i < width; i += 3 )
 137             {
 138                 s0 += (ST)S[i + ksz_cn] - (ST)S[i];
 139                 s1 += (ST)S[i + ksz_cn + 1] - (ST)S[i + 1];
 140                 s2 += (ST)S[i + ksz_cn + 2] - (ST)S[i + 2];
 141                 D[i+3] = s0;
 142                 D[i+4] = s1;
 143                 D[i+5] = s2;
 144             }
 145         }
 146         else if( cn == 4 )
 147         {
 148             ST s0 = 0, s1 = 0, s2 = 0, s3 = 0;
 149             for( i = 0; i < ksz_cn; i += 4 )
 150             {
 151                 s0 += (ST)S[i];
 152                 s1 += (ST)S[i+1];
 153                 s2 += (ST)S[i+2];
 154                 s3 += (ST)S[i+3];
 155             }
 156             D[0] = s0;
 157             D[1] = s1;
 158             D[2] = s2;
 159             D[3] = s3;
 160             for( i = 0; i < width; i += 4 )
 161             {
 162                 s0 += (ST)S[i + ksz_cn] - (ST)S[i];
 163                 s1 += (ST)S[i + ksz_cn + 1] - (ST)S[i + 1];
 164                 s2 += (ST)S[i + ksz_cn + 2] - (ST)S[i + 2];
 165                 s3 += (ST)S[i + ksz_cn + 3] - (ST)S[i + 3];
 166                 D[i+4] = s0;
 167                 D[i+5] = s1;
 168                 D[i+6] = s2;
 169                 D[i+7] = s3;
 170             }
 171         }
 172         else
 173             for( k = 0; k < cn; k++, S++, D++ )
 174             {
 175                 ST s = 0;
 176                 for( i = 0; i < ksz_cn; i += cn )
 177                     s += (ST)S[i];
 178                 D[0] = s;
 179                 for( i = 0; i < width; i += cn )
 180                 {
 181                     s += (ST)S[i + ksz_cn] - (ST)S[i];
 182                     D[i+cn] = s;
 183                 }
 184             }
 185     }
 186 };
 187
 188
 189 template<typename ST, typename T>
 190 struct ColumnSum :
 191         public BaseColumnFilter
 192 {
 193     ColumnSum( int _ksize, int _anchor, double _scale ) :
 194         BaseColumnFilter()
 195     {
 196         ksize = _ksize;
 197         anchor = _anchor;
 198         scale = _scale;
 199         sumCount = 0;
 200     }
 201
 202     virtual void reset() { sumCount = 0; }
 203
 204     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 205     {
 206         int i;
 207         ST* SUM;
 208         bool haveScale = scale != 1;
 209         double _scale = scale;
 210
 211         if( width != (int)sum.size() )
 212         {
 213             sum.resize(width);
 214             sumCount = 0;
 215         }
 216
 217         SUM = &sum[0];
 218         if( sumCount == 0 )
 219         {
 220             memset((void*)SUM, 0, width*sizeof(ST));
 221
 222             for( ; sumCount < ksize - 1; sumCount++, src++ )
 223             {
 224                 const ST* Sp = (const ST*)src[0];
 225
 226                 for( i = 0; i < width; i++ )
 227                     SUM[i] += Sp[i];
 228             }
 229         }
 230         else
 231         {
 232             CV_Assert( sumCount == ksize-1 );
 233             src += ksize-1;
 234         }
 235
 236         for( ; count--; src++ )
 237         {
 238             const ST* Sp = (const ST*)src[0];
 239             const ST* Sm = (const ST*)src[1-ksize];
 240             T* D = (T*)dst;
 241             if( haveScale )
 242             {
 243                 for( i = 0; i <= width - 2; i += 2 )
 244                 {
 245                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 246                     D[i] = saturate_cast<T>(s0*_scale);
 247                     D[i+1] = saturate_cast<T>(s1*_scale);
 248                     s0 -= Sm[i]; s1 -= Sm[i+1];
 249                     SUM[i] = s0; SUM[i+1] = s1;
 250                 }
 251
 252                 for( ; i < width; i++ )
 253                 {
 254                     ST s0 = SUM[i] + Sp[i];
 255                     D[i] = saturate_cast<T>(s0*_scale);
 256                     SUM[i] = s0 - Sm[i];
 257                 }
 258             }
 259             else
 260             {
 261                 for( i = 0; i <= width - 2; i += 2 )
 262                 {
 263                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 264                     D[i] = saturate_cast<T>(s0);
 265                     D[i+1] = saturate_cast<T>(s1);
 266                     s0 -= Sm[i]; s1 -= Sm[i+1];
 267                     SUM[i] = s0; SUM[i+1] = s1;
 268                 }
 269
 270                 for( ; i < width; i++ )
 271                 {
 272                     ST s0 = SUM[i] + Sp[i];
 273                     D[i] = saturate_cast<T>(s0);
 274                     SUM[i] = s0 - Sm[i];
 275                 }
 276             }
 277             dst += dststep;
 278         }
 279     }
 280
 281     double scale;
 282     int sumCount;
 283     std::vector<ST> sum;
 284 };
 285
 286
 287 template<>
 288 struct ColumnSum<int, uchar> :
 289         public BaseColumnFilter
 290 {
 291     ColumnSum( int _ksize, int _anchor, double _scale ) :
 292         BaseColumnFilter()
 293     {
 294         ksize = _ksize;
 295         anchor = _anchor;
 296         scale = _scale;
 297         sumCount = 0;
 298     }
 299
 300     virtual void reset() { sumCount = 0; }
 301
 302     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 303     {
 304         int* SUM;
 305         bool haveScale = scale != 1;
 306         double _scale = scale;
 307
 308 #if CV_SIMD128
 309         bool haveSIMD128 = hasSIMD128();
 310 #endif
 311
 312         if( width != (int)sum.size() )
 313         {
 314             sum.resize(width);
 315             sumCount = 0;
 316         }
 317
 318         SUM = &sum[0];
 319         if( sumCount == 0 )
 320         {
 321             memset((void*)SUM, 0, width*sizeof(int));
 322             for( ; sumCount < ksize - 1; sumCount++, src++ )
 323             {
 324                 const int* Sp = (const int*)src[0];
 325                 int i = 0;
 326 #if CV_SIMD128
 327                 if( haveSIMD128 )
 328                 {
 329                     for (; i <= width - 4; i += 4)
 330                     {
 331                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 332                     }
 333                 }
 334 #endif
 335                 for( ; i < width; i++ )
 336                     SUM[i] += Sp[i];
 337             }
 338         }
 339         else
 340         {
 341             CV_Assert( sumCount == ksize-1 );
 342             src += ksize-1;
 343         }
 344
 345         for( ; count--; src++ )
 346         {
 347             const int* Sp = (const int*)src[0];
 348             const int* Sm = (const int*)src[1-ksize];
 349             uchar* D = (uchar*)dst;
 350             if( haveScale )
 351             {
 352                 int i = 0;
 353 #if CV_SIMD128
 354                 if( haveSIMD128 )
 355                 {
 356
 357                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 358                     for( ; i <= width-8; i+=8 )
 359                     {
 360                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 361                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 362
 363                         v_uint32x4 v_s0d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s0) * v_scale));
 364                         v_uint32x4 v_s01d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s01) * v_scale));
 365
 366                         v_uint16x8 v_dst = v_pack(v_s0d, v_s01d);
 367                         v_pack_store(D + i, v_dst);
 368
 369                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 370                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 371                     }
 372                 }
 373 #endif
 374                 for( ; i < width; i++ )
 375                 {
 376                     int s0 = SUM[i] + Sp[i];
 377                     D[i] = saturate_cast<uchar>(s0*_scale);
 378                     SUM[i] = s0 - Sm[i];
 379                 }
 380             }
 381             else
 382             {
 383                 int i = 0;
 384 #if CV_SIMD128
 385                 if( haveSIMD128 )
 386                 {
 387                     for( ; i <= width-8; i+=8 )
 388                     {
 389                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 390                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 391
 392                         v_uint16x8 v_dst = v_pack(v_reinterpret_as_u32(v_s0), v_reinterpret_as_u32(v_s01));
 393                         v_pack_store(D + i, v_dst);
 394
 395                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 396                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 397                     }
 398                 }
 399 #endif
 400
 401                 for( ; i < width; i++ )
 402                 {
 403                     int s0 = SUM[i] + Sp[i];
 404                     D[i] = saturate_cast<uchar>(s0);
 405                     SUM[i] = s0 - Sm[i];
 406                 }
 407             }
 408             dst += dststep;
 409         }
 410     }
 411
 412     double scale;
 413     int sumCount;
 414     std::vector<int> sum;
 415 };
 416
 417
 418 template<>
 419 struct ColumnSum<ushort, uchar> :
 420 public BaseColumnFilter
 421 {
 422     enum { SHIFT = 23 };
 423
 424     ColumnSum( int _ksize, int _anchor, double _scale ) :
 425     BaseColumnFilter()
 426     {
 427         ksize = _ksize;
 428         anchor = _anchor;
 429         scale = _scale;
 430         sumCount = 0;
 431         divDelta = 0;
 432         divScale = 1;
 433         if( scale != 1 )
 434         {
 435             int d = cvRound(1./scale);
 436             double scalef = ((double)(1 << SHIFT))/d;
 437             divScale = cvFloor(scalef);
 438             scalef -= divScale;
 439             divDelta = d/2;
 440             if( scalef < 0.5 )
 441                 divDelta++;
 442             else
 443                 divScale++;
 444         }
 445     }
 446
 447     virtual void reset() { sumCount = 0; }
 448
 449     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 450     {
 451         const int ds = divScale;
 452         const int dd = divDelta;
 453         ushort* SUM;
 454         const bool haveScale = scale != 1;
 455
 456 #if CV_SIMD128
 457         bool haveSIMD128 = hasSIMD128();
 458 #endif
 459
 460         if( width != (int)sum.size() )
 461         {
 462             sum.resize(width);
 463             sumCount = 0;
 464         }
 465
 466         SUM = &sum[0];
 467         if( sumCount == 0 )
 468         {
 469             memset((void*)SUM, 0, width*sizeof(SUM[0]));
 470             for( ; sumCount < ksize - 1; sumCount++, src++ )
 471             {
 472                 const ushort* Sp = (const ushort*)src[0];
 473                 int i = 0;
 474 #if CV_SIMD128
 475                 if( haveSIMD128 )
 476                 {
 477                     for( ; i <= width - 8; i += 8 )
 478                     {
 479                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 480                     }
 481                 }
 482 #endif
 483                 for( ; i < width; i++ )
 484                     SUM[i] += Sp[i];
 485             }
 486         }
 487         else
 488         {
 489             CV_Assert( sumCount == ksize-1 );
 490             src += ksize-1;
 491         }
 492
 493         for( ; count--; src++ )
 494         {
 495             const ushort* Sp = (const ushort*)src[0];
 496             const ushort* Sm = (const ushort*)src[1-ksize];
 497             uchar* D = (uchar*)dst;
 498             if( haveScale )
 499             {
 500                 int i = 0;
 501             #if CV_SIMD128
 502                 v_uint32x4 ds4 = v_setall_u32((unsigned)ds);
 503                 v_uint16x8 dd8 = v_setall_u16((ushort)dd);
 504
 505                 for( ; i <= width-16; i+=16 )
 506                 {
 507                     v_uint16x8 _sm0 = v_load(Sm + i);
 508                     v_uint16x8 _sm1 = v_load(Sm + i + 8);
 509
 510                     v_uint16x8 _s0 = v_add_wrap(v_load(SUM + i), v_load(Sp + i));
 511                     v_uint16x8 _s1 = v_add_wrap(v_load(SUM + i + 8), v_load(Sp + i + 8));
 512
 513                     v_uint32x4 _s00, _s01, _s10, _s11;
 514
 515                     v_expand(_s0 + dd8, _s00, _s01);
 516                     v_expand(_s1 + dd8, _s10, _s11);
 517
 518                     _s00 = v_shr<SHIFT>(_s00*ds4);
 519                     _s01 = v_shr<SHIFT>(_s01*ds4);
 520                     _s10 = v_shr<SHIFT>(_s10*ds4);
 521                     _s11 = v_shr<SHIFT>(_s11*ds4);
 522
 523                     v_int16x8 r0 = v_pack(v_reinterpret_as_s32(_s00), v_reinterpret_as_s32(_s01));
 524                     v_int16x8 r1 = v_pack(v_reinterpret_as_s32(_s10), v_reinterpret_as_s32(_s11));
 525
 526                     _s0 = v_sub_wrap(_s0, _sm0);
 527                     _s1 = v_sub_wrap(_s1, _sm1);
 528
 529                     v_store(D + i, v_pack_u(r0, r1));
 530                     v_store(SUM + i, _s0);
 531                     v_store(SUM + i + 8, _s1);
 532                 }
 533             #endif
 534                 for( ; i < width; i++ )
 535                 {
 536                     int s0 = SUM[i] + Sp[i];
 537                     D[i] = (uchar)((s0 + dd)*ds >> SHIFT);
 538                     SUM[i] = (ushort)(s0 - Sm[i]);
 539                 }
 540             }
 541             else
 542             {
 543                 int i = 0;
 544                 for( ; i < width; i++ )
 545                 {
 546                     int s0 = SUM[i] + Sp[i];
 547                     D[i] = saturate_cast<uchar>(s0);
 548                     SUM[i] = (ushort)(s0 - Sm[i]);
 549                 }
 550             }
 551             dst += dststep;
 552         }
 553     }
 554
 555     double scale;
 556     int sumCount;
 557     int divDelta;
 558     int divScale;
 559     std::vector<ushort> sum;
 560 };
 561
 562
 563 template<>
 564 struct ColumnSum<int, short> :
 565         public BaseColumnFilter
 566 {
 567     ColumnSum( int _ksize, int _anchor, double _scale ) :
 568         BaseColumnFilter()
 569     {
 570         ksize = _ksize;
 571         anchor = _anchor;
 572         scale = _scale;
 573         sumCount = 0;
 574     }
 575
 576     virtual void reset() { sumCount = 0; }
 577
 578     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 579     {
 580         int i;
 581         int* SUM;
 582         bool haveScale = scale != 1;
 583         double _scale = scale;
 584
 585 #if CV_SIMD128
 586         bool haveSIMD128 = hasSIMD128();
 587 #endif
 588
 589         if( width != (int)sum.size() )
 590         {
 591             sum.resize(width);
 592             sumCount = 0;
 593         }
 594
 595         SUM = &sum[0];
 596         if( sumCount == 0 )
 597         {
 598             memset((void*)SUM, 0, width*sizeof(int));
 599             for( ; sumCount < ksize - 1; sumCount++, src++ )
 600             {
 601                 const int* Sp = (const int*)src[0];
 602                 i = 0;
 603 #if CV_SIMD128
 604                 if( haveSIMD128 )
 605                 {
 606                     for( ; i <= width - 4; i+=4 )
 607                     {
 608                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 609                     }
 610                 }
 611                 #endif
 612                 for( ; i < width; i++ )
 613                     SUM[i] += Sp[i];
 614             }
 615         }
 616         else
 617         {
 618             CV_Assert( sumCount == ksize-1 );
 619             src += ksize-1;
 620         }
 621
 622         for( ; count--; src++ )
 623         {
 624             const int* Sp = (const int*)src[0];
 625             const int* Sm = (const int*)src[1-ksize];
 626             short* D = (short*)dst;
 627             if( haveScale )
 628             {
 629                 i = 0;
 630 #if CV_SIMD128
 631                 if( haveSIMD128 )
 632                 {
 633                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 634                     for( ; i <= width-8; i+=8 )
 635                     {
 636                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 637                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 638
 639                         v_int32x4 v_s0d =  v_round(v_cvt_f32(v_s0) * v_scale);
 640                         v_int32x4 v_s01d = v_round(v_cvt_f32(v_s01) * v_scale);
 641                         v_store(D + i, v_pack(v_s0d, v_s01d));
 642
 643                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 644                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 645                     }
 646                 }
 647 #endif
 648                 for( ; i < width; i++ )
 649                 {
 650                     int s0 = SUM[i] + Sp[i];
 651                     D[i] = saturate_cast<short>(s0*_scale);
 652                     SUM[i] = s0 - Sm[i];
 653                 }
 654             }
 655             else
 656             {
 657                 i = 0;
 658 #if CV_SIMD128
 659                 if( haveSIMD128 )
 660                 {
 661                     for( ; i <= width-8; i+=8 )
 662                     {
 663                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 664                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 665
 666                         v_store(D + i, v_pack(v_s0, v_s01));
 667
 668                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 669                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 670                     }
 671                 }
 672 #endif
 673
 674                 for( ; i < width; i++ )
 675                 {
 676                     int s0 = SUM[i] + Sp[i];
 677                     D[i] = saturate_cast<short>(s0);
 678                     SUM[i] = s0 - Sm[i];
 679                 }
 680             }
 681             dst += dststep;
 682         }
 683     }
 684
 685     double scale;
 686     int sumCount;
 687     std::vector<int> sum;
 688 };
 689
 690
 691 template<>
 692 struct ColumnSum<int, ushort> :
 693         public BaseColumnFilter
 694 {
 695     ColumnSum( int _ksize, int _anchor, double _scale ) :
 696         BaseColumnFilter()
 697     {
 698         ksize = _ksize;
 699         anchor = _anchor;
 700         scale = _scale;
 701         sumCount = 0;
 702     }
 703
 704     virtual void reset() { sumCount = 0; }
 705
 706     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 707     {
 708         int* SUM;
 709         bool haveScale = scale != 1;
 710         double _scale = scale;
 711
 712 #if CV_SIMD128
 713         bool haveSIMD128 = hasSIMD128();
 714 #endif
 715
 716         if( width != (int)sum.size() )
 717         {
 718             sum.resize(width);
 719             sumCount = 0;
 720         }
 721
 722         SUM = &sum[0];
 723         if( sumCount == 0 )
 724         {
 725             memset((void*)SUM, 0, width*sizeof(int));
 726             for( ; sumCount < ksize - 1; sumCount++, src++ )
 727             {
 728                 const int* Sp = (const int*)src[0];
 729                 int i = 0;
 730 #if CV_SIMD128
 731                 if( haveSIMD128 )
 732                 {
 733                     for (; i <= width - 4; i += 4)
 734                     {
 735                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 736                     }
 737                 }
 738 #endif
 739                 for( ; i < width; i++ )
 740                     SUM[i] += Sp[i];
 741             }
 742         }
 743         else
 744         {
 745             CV_Assert( sumCount == ksize-1 );
 746             src += ksize-1;
 747         }
 748
 749         for( ; count--; src++ )
 750         {
 751             const int* Sp = (const int*)src[0];
 752             const int* Sm = (const int*)src[1-ksize];
 753             ushort* D = (ushort*)dst;
 754             if( haveScale )
 755             {
 756                 int i = 0;
 757 #if CV_SIMD128
 758                 if( haveSIMD128 )
 759                 {
 760                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 761                     for( ; i <= width-8; i+=8 )
 762                     {
 763                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 764                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 765
 766                         v_uint32x4 v_s0d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s0) * v_scale));
 767                         v_uint32x4 v_s01d = v_reinterpret_as_u32(v_round(v_cvt_f32(v_s01) * v_scale));
 768                         v_store(D + i, v_pack(v_s0d, v_s01d));
 769
 770                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 771                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 772                     }
 773                 }
 774 #endif
 775                 for( ; i < width; i++ )
 776                 {
 777                     int s0 = SUM[i] + Sp[i];
 778                     D[i] = saturate_cast<ushort>(s0*_scale);
 779                     SUM[i] = s0 - Sm[i];
 780                 }
 781             }
 782             else
 783             {
 784                 int i = 0;
 785 #if CV_SIMD128
 786                 if( haveSIMD128 )
 787                 {
 788                     for( ; i <= width-8; i+=8 )
 789                     {
 790                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 791                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
 792
 793                         v_store(D + i, v_pack(v_reinterpret_as_u32(v_s0), v_reinterpret_as_u32(v_s01)));
 794
 795                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 796                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
 797                     }
 798                 }
 799 #endif
 800                 for( ; i < width; i++ )
 801                 {
 802                     int s0 = SUM[i] + Sp[i];
 803                     D[i] = saturate_cast<ushort>(s0);
 804                     SUM[i] = s0 - Sm[i];
 805                 }
 806             }
 807             dst += dststep;
 808         }
 809     }
 810
 811     double scale;
 812     int sumCount;
 813     std::vector<int> sum;
 814 };
 815
 816 template<>
 817 struct ColumnSum<int, int> :
 818         public BaseColumnFilter
 819 {
 820     ColumnSum( int _ksize, int _anchor, double _scale ) :
 821         BaseColumnFilter()
 822     {
 823         ksize = _ksize;
 824         anchor = _anchor;
 825         scale = _scale;
 826         sumCount = 0;
 827     }
 828
 829     virtual void reset() { sumCount = 0; }
 830
 831     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 832     {
 833         int* SUM;
 834         bool haveScale = scale != 1;
 835         double _scale = scale;
 836
 837 #if CV_SIMD128
 838         bool haveSIMD128 = hasSIMD128();
 839 #endif
 840
 841         if( width != (int)sum.size() )
 842         {
 843             sum.resize(width);
 844             sumCount = 0;
 845         }
 846
 847         SUM = &sum[0];
 848         if( sumCount == 0 )
 849         {
 850             memset((void*)SUM, 0, width*sizeof(int));
 851             for( ; sumCount < ksize - 1; sumCount++, src++ )
 852             {
 853                 const int* Sp = (const int*)src[0];
 854                 int i = 0;
 855 #if CV_SIMD128
 856                 if( haveSIMD128 )
 857                 {
 858                     for( ; i <= width - 4; i+=4 )
 859                     {
 860                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 861                     }
 862                 }
 863 #endif
 864                 for( ; i < width; i++ )
 865                     SUM[i] += Sp[i];
 866             }
 867         }
 868         else
 869         {
 870             CV_Assert( sumCount == ksize-1 );
 871             src += ksize-1;
 872         }
 873
 874         for( ; count--; src++ )
 875         {
 876             const int* Sp = (const int*)src[0];
 877             const int* Sm = (const int*)src[1-ksize];
 878             int* D = (int*)dst;
 879             if( haveScale )
 880             {
 881                 int i = 0;
 882 #if CV_SIMD128
 883                 if( haveSIMD128 )
 884                 {
 885                     v_float32x4 v_scale = v_setall_f32((float)_scale);
 886                     for( ; i <= width-4; i+=4 )
 887                     {
 888                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 889                         v_int32x4 v_s0d = v_round(v_cvt_f32(v_s0) * v_scale);
 890
 891                         v_store(D + i, v_s0d);
 892                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 893                     }
 894                 }
 895 #endif
 896                 for( ; i < width; i++ )
 897                 {
 898                     int s0 = SUM[i] + Sp[i];
 899                     D[i] = saturate_cast<int>(s0*_scale);
 900                     SUM[i] = s0 - Sm[i];
 901                 }
 902             }
 903             else
 904             {
 905                 int i = 0;
 906 #if CV_SIMD128
 907                 if( haveSIMD128 )
 908                 {
 909                     for( ; i <= width-4; i+=4 )
 910                     {
 911                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
 912
 913                         v_store(D + i, v_s0);
 914                         v_store(SUM + i, v_s0 - v_load(Sm + i));
 915                     }
 916                 }
 917 #endif
 918                 for( ; i < width; i++ )
 919                 {
 920                     int s0 = SUM[i] + Sp[i];
 921                     D[i] = s0;
 922                     SUM[i] = s0 - Sm[i];
 923                 }
 924             }
 925             dst += dststep;
 926         }
 927     }
 928
 929     double scale;
 930     int sumCount;
 931     std::vector<int> sum;
 932 };
 933
 934
 935 template<>
 936 struct ColumnSum<int, float> :
 937         public BaseColumnFilter
 938 {
 939     ColumnSum( int _ksize, int _anchor, double _scale ) :
 940         BaseColumnFilter()
 941     {
 942         ksize = _ksize;
 943         anchor = _anchor;
 944         scale = _scale;
 945         sumCount = 0;
 946     }
 947
 948     virtual void reset() { sumCount = 0; }
 949
 950     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 951     {
 952         int* SUM;
 953         bool haveScale = scale != 1;
 954         double _scale = scale;
 955
 956 #if CV_SIMD128
 957         bool haveSIMD128 = hasSIMD128();
 958 #endif
 959
 960         if( width != (int)sum.size() )
 961         {
 962             sum.resize(width);
 963             sumCount = 0;
 964         }
 965
 966         SUM = &sum[0];
 967         if( sumCount == 0 )
 968         {
 969             memset((void*)SUM, 0, width*sizeof(int));
 970             for( ; sumCount < ksize - 1; sumCount++, src++ )
 971             {
 972                 const int* Sp = (const int*)src[0];
 973                 int i = 0;
 974 #if CV_SIMD128
 975                 if( haveSIMD128 )
 976                 {
 977                     for( ; i <= width - 4; i+=4 )
 978                     {
 979                         v_store(SUM + i, v_load(SUM + i) + v_load(Sp + i));
 980                     }
 981                 }
 982 #endif
 983
 984                 for( ; i < width; i++ )
 985                     SUM[i] += Sp[i];
 986             }
 987         }
 988         else
 989         {
 990             CV_Assert( sumCount == ksize-1 );
 991             src += ksize-1;
 992         }
 993
 994         for( ; count--; src++ )
 995         {
 996             const int * Sp = (const int*)src[0];
 997             const int * Sm = (const int*)src[1-ksize];
 998             float* D = (float*)dst;
 999             if( haveScale )
1000             {
1001                 int i = 0;
1002
1003 #if CV_SIMD128
1004                 if( haveSIMD128 )
1005                 {
1006                     v_float32x4 v_scale = v_setall_f32((float)_scale);
1007                     for (; i <= width - 8; i += 8)
1008                     {
1009                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
1010                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
1011
1012                         v_store(D + i, v_cvt_f32(v_s0) * v_scale);
1013                         v_store(D + i + 4, v_cvt_f32(v_s01) * v_scale);
1014
1015                         v_store(SUM + i, v_s0 - v_load(Sm + i));
1016                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
1017                     }
1018                 }
1019 #endif
1020                 for( ; i < width; i++ )
1021                 {
1022                     int s0 = SUM[i] + Sp[i];
1023                     D[i] = (float)(s0*_scale);
1024                     SUM[i] = s0 - Sm[i];
1025                 }
1026             }
1027             else
1028             {
1029                 int i = 0;
1030
1031 #if CV_SIMD128
1032                 if( haveSIMD128 )
1033                 {
1034                     for( ; i <= width-8; i+=8 )
1035                     {
1036                         v_int32x4 v_s0 = v_load(SUM + i) + v_load(Sp + i);
1037                         v_int32x4 v_s01 = v_load(SUM + i + 4) + v_load(Sp + i + 4);
1038
1039                         v_store(D + i, v_cvt_f32(v_s0));
1040                         v_store(D + i + 4, v_cvt_f32(v_s01));
1041
1042                         v_store(SUM + i, v_s0 - v_load(Sm + i));
1043                         v_store(SUM + i + 4, v_s01 - v_load(Sm + i + 4));
1044                     }
1045                 }
1046 #endif
1047                 for( ; i < width; i++ )
1048                 {
1049                     int s0 = SUM[i] + Sp[i];
1050                     D[i] = (float)(s0);
1051                     SUM[i] = s0 - Sm[i];
1052                 }
1053             }
1054             dst += dststep;
1055         }
1056     }
1057
1058     double scale;
1059     int sumCount;
1060     std::vector<int> sum;
1061 };
1062
1063 #ifdef HAVE_OPENCL
1064
1065 static bool ocl_boxFilter3x3_8UC1( InputArray _src, OutputArray _dst, int ddepth,
1066                                    Size ksize, Point anchor, int borderType, bool normalize )
1067 {
1068     const ocl::Device & dev = ocl::Device::getDefault();
1069     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
1070
1071     if (ddepth < 0)
1072         ddepth = sdepth;
1073
1074     if (anchor.x < 0)
1075         anchor.x = ksize.width / 2;
1076     if (anchor.y < 0)
1077         anchor.y = ksize.height / 2;
1078
1079     if ( !(dev.isIntel() && (type == CV_8UC1) &&
1080          (_src.offset() == 0) && (_src.step() % 4 == 0) &&
1081          (_src.cols() % 16 == 0) && (_src.rows() % 2 == 0) &&
1082          (anchor.x == 1) && (anchor.y == 1) &&
1083          (ksize.width == 3) && (ksize.height == 3)) )
1084         return false;
1085
1086     float alpha = 1.0f / (ksize.height * ksize.width);
1087     Size size = _src.size();
1088     size_t globalsize[2] = { 0, 0 };
1089     size_t localsize[2] = { 0, 0 };
1090     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
1091
1092     globalsize[0] = size.width / 16;
1093     globalsize[1] = size.height / 2;
1094
1095     char build_opts[1024];
1096     sprintf(build_opts, "-D %s %s", borderMap[borderType], normalize ? "-D NORMALIZE" : "");
1097
1098     ocl::Kernel kernel("boxFilter3x3_8UC1_cols16_rows2", cv::ocl::imgproc::boxFilter3x3_oclsrc, build_opts);
1099     if (kernel.empty())
1100         return false;
1101
1102     UMat src = _src.getUMat();
1103     _dst.create(size, CV_MAKETYPE(ddepth, cn));
1104     if (!(_dst.offset() == 0 && _dst.step() % 4 == 0))
1105         return false;
1106     UMat dst = _dst.getUMat();
1107
1108     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
1109     idxArg = kernel.set(idxArg, (int)src.step);
1110     idxArg = kernel.set(idxArg, ocl::KernelArg::PtrWriteOnly(dst));
1111     idxArg = kernel.set(idxArg, (int)dst.step);
1112     idxArg = kernel.set(idxArg, (int)dst.rows);
1113     idxArg = kernel.set(idxArg, (int)dst.cols);
1114     if (normalize)
1115         idxArg = kernel.set(idxArg, (float)alpha);
1116
1117     return kernel.run(2, globalsize, (localsize[0] == 0) ? NULL : localsize, false);
1118 }
1119
1120 #define DIVUP(total, grain) ((total + grain - 1) / (grain))
1121 #define ROUNDUP(sz, n)      ((sz) + (n) - 1 - (((sz) + (n) - 1) % (n)))
1122
1123 static bool ocl_boxFilter( InputArray _src, OutputArray _dst, int ddepth,
1124                            Size ksize, Point anchor, int borderType, bool normalize, bool sqr = false )
1125 {
1126     const ocl::Device & dev = ocl::Device::getDefault();
1127     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type), esz = CV_ELEM_SIZE(type);
1128     bool doubleSupport = dev.doubleFPConfig() > 0;
1129
1130     if (ddepth < 0)
1131         ddepth = sdepth;
1132
1133     if (cn > 4 || (!doubleSupport && (sdepth == CV_64F || ddepth == CV_64F)) ||
1134         _src.offset() % esz != 0 || _src.step() % esz != 0)
1135         return false;
1136
1137     if (anchor.x < 0)
1138         anchor.x = ksize.width / 2;
1139     if (anchor.y < 0)
1140         anchor.y = ksize.height / 2;
1141
1142     int computeUnits = ocl::Device::getDefault().maxComputeUnits();
1143     float alpha = 1.0f / (ksize.height * ksize.width);
1144     Size size = _src.size(), wholeSize;
1145     bool isolated = (borderType & BORDER_ISOLATED) != 0;
1146     borderType &= ~BORDER_ISOLATED;
1147     int wdepth = std::max(CV_32F, std::max(ddepth, sdepth)),
1148         wtype = CV_MAKE_TYPE(wdepth, cn), dtype = CV_MAKE_TYPE(ddepth, cn);
1149
1150     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
1151     size_t globalsize[2] = { (size_t)size.width, (size_t)size.height };
1152     size_t localsize_general[2] = { 0, 1 }, * localsize = NULL;
1153
1154     UMat src = _src.getUMat();
1155     if (!isolated)
1156     {
1157         Point ofs;
1158         src.locateROI(wholeSize, ofs);
1159     }
1160
1161     int h = isolated ? size.height : wholeSize.height;
1162     int w = isolated ? size.width : wholeSize.width;
1163
1164     size_t maxWorkItemSizes[32];
1165     ocl::Device::getDefault().maxWorkItemSizes(maxWorkItemSizes);
1166     int tryWorkItems = (int)maxWorkItemSizes[0];
1167
1168     ocl::Kernel kernel;
1169
1170     if (dev.isIntel() && !(dev.type() & ocl::Device::TYPE_CPU) &&
1171         ((ksize.width < 5 && ksize.height < 5 && esz <= 4) ||
1172          (ksize.width == 5 && ksize.height == 5 && cn == 1)))
1173     {
1174         if (w < ksize.width || h < ksize.height)
1175             return false;
1176
1177         // Figure out what vector size to use for loading the pixels.
1178         int pxLoadNumPixels = cn != 1 || size.width % 4 ? 1 : 4;
1179         int pxLoadVecSize = cn * pxLoadNumPixels;
1180
1181         // Figure out how many pixels per work item to compute in X and Y
1182         // directions.  Too many and we run out of registers.
1183         int pxPerWorkItemX = 1, pxPerWorkItemY = 1;
1184         if (cn <= 2 && ksize.width <= 4 && ksize.height <= 4)
1185         {
1186             pxPerWorkItemX = size.width % 8 ? size.width % 4 ? size.width % 2 ? 1 : 2 : 4 : 8;
1187             pxPerWorkItemY = size.height % 2 ? 1 : 2;
1188         }
1189         else if (cn < 4 || (ksize.width <= 4 && ksize.height <= 4))
1190         {
1191             pxPerWorkItemX = size.width % 2 ? 1 : 2;
1192             pxPerWorkItemY = size.height % 2 ? 1 : 2;
1193         }
1194         globalsize[0] = size.width / pxPerWorkItemX;
1195         globalsize[1] = size.height / pxPerWorkItemY;
1196
1197         // Need some padding in the private array for pixels
1198         int privDataWidth = ROUNDUP(pxPerWorkItemX + ksize.width - 1, pxLoadNumPixels);
1199
1200         // Make the global size a nice round number so the runtime can pick
1201         // from reasonable choices for the workgroup size
1202         const int wgRound = 256;
1203         globalsize[0] = ROUNDUP(globalsize[0], wgRound);
1204
1205         char build_options[1024], cvt[2][40];
1206         sprintf(build_options, "-D cn=%d "
1207                 "-D ANCHOR_X=%d -D ANCHOR_Y=%d -D KERNEL_SIZE_X=%d -D KERNEL_SIZE_Y=%d "
1208                 "-D PX_LOAD_VEC_SIZE=%d -D PX_LOAD_NUM_PX=%d "
1209                 "-D PX_PER_WI_X=%d -D PX_PER_WI_Y=%d -D PRIV_DATA_WIDTH=%d -D %s -D %s "
1210                 "-D PX_LOAD_X_ITERATIONS=%d -D PX_LOAD_Y_ITERATIONS=%d "
1211                 "-D srcT=%s -D srcT1=%s -D dstT=%s -D dstT1=%s -D WT=%s -D WT1=%s "
1212                 "-D convertToWT=%s -D convertToDstT=%s%s%s -D PX_LOAD_FLOAT_VEC_CONV=convert_%s -D OP_BOX_FILTER",
1213                 cn, anchor.x, anchor.y, ksize.width, ksize.height,
1214                 pxLoadVecSize, pxLoadNumPixels,
1215                 pxPerWorkItemX, pxPerWorkItemY, privDataWidth, borderMap[borderType],
1216                 isolated ? "BORDER_ISOLATED" : "NO_BORDER_ISOLATED",
1217                 privDataWidth / pxLoadNumPixels, pxPerWorkItemY + ksize.height - 1,
1218                 ocl::typeToStr(type), ocl::typeToStr(sdepth), ocl::typeToStr(dtype),
1219                 ocl::typeToStr(ddepth), ocl::typeToStr(wtype), ocl::typeToStr(wdepth),
1220                 ocl::convertTypeStr(sdepth, wdepth, cn, cvt[0]),
1221                 ocl::convertTypeStr(wdepth, ddepth, cn, cvt[1]),
1222                 normalize ? " -D NORMALIZE" : "", sqr ? " -D SQR" : "",
1223                 ocl::typeToStr(CV_MAKE_TYPE(wdepth, pxLoadVecSize)) //PX_LOAD_FLOAT_VEC_CONV
1224                 );
1225
1226
1227         if (!kernel.create("filterSmall", cv::ocl::imgproc::filterSmall_oclsrc, build_options))
1228             return false;
1229     }
1230     else
1231     {
1232         localsize = localsize_general;
1233         for ( ; ; )
1234         {
1235             int BLOCK_SIZE_X = tryWorkItems, BLOCK_SIZE_Y = std::min(ksize.height * 10, size.height);
1236
1237             while (BLOCK_SIZE_X > 32 && BLOCK_SIZE_X >= ksize.width * 2 && BLOCK_SIZE_X > size.width * 2)
1238                 BLOCK_SIZE_X /= 2;
1239             while (BLOCK_SIZE_Y < BLOCK_SIZE_X / 8 && BLOCK_SIZE_Y * computeUnits * 32 < size.height)
1240                 BLOCK_SIZE_Y *= 2;
1241
1242             if (ksize.width > BLOCK_SIZE_X || w < ksize.width || h < ksize.height)
1243                 return false;
1244
1245             char cvt[2][50];
1246             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"
1247                                  " -D ANCHOR_X=%d -D ANCHOR_Y=%d -D KERNEL_SIZE_X=%d -D KERNEL_SIZE_Y=%d -D %s%s%s%s%s"
1248                                  " -D ST1=%s -D DT1=%s -D cn=%d",
1249                                  BLOCK_SIZE_X, BLOCK_SIZE_Y, ocl::typeToStr(type), ocl::typeToStr(CV_MAKE_TYPE(ddepth, cn)),
1250                                  ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)),
1251                                  ocl::convertTypeStr(wdepth, ddepth, cn, cvt[0]),
1252                                  ocl::convertTypeStr(sdepth, wdepth, cn, cvt[1]),
1253                                  anchor.x, anchor.y, ksize.width, ksize.height, borderMap[borderType],
1254                                  isolated ? " -D BORDER_ISOLATED" : "", doubleSupport ? " -D DOUBLE_SUPPORT" : "",
1255                                  normalize ? " -D NORMALIZE" : "", sqr ? " -D SQR" : "",
1256                                  ocl::typeToStr(sdepth), ocl::typeToStr(ddepth), cn);
1257
1258             localsize[0] = BLOCK_SIZE_X;
1259             globalsize[0] = DIVUP(size.width, BLOCK_SIZE_X - (ksize.width - 1)) * BLOCK_SIZE_X;
1260             globalsize[1] = DIVUP(size.height, BLOCK_SIZE_Y);
1261
1262             kernel.create("boxFilter", cv::ocl::imgproc::boxFilter_oclsrc, opts);
1263             if (kernel.empty())
1264                 return false;
1265
1266             size_t kernelWorkGroupSize = kernel.workGroupSize();
1267             if (localsize[0] <= kernelWorkGroupSize)
1268                 break;
1269             if (BLOCK_SIZE_X < (int)kernelWorkGroupSize)
1270                 return false;
1271
1272             tryWorkItems = (int)kernelWorkGroupSize;
1273         }
1274     }
1275
1276     _dst.create(size, CV_MAKETYPE(ddepth, cn));
1277     UMat dst = _dst.getUMat();
1278
1279     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
1280     idxArg = kernel.set(idxArg, (int)src.step);
1281     int srcOffsetX = (int)((src.offset % src.step) / src.elemSize());
1282     int srcOffsetY = (int)(src.offset / src.step);
1283     int srcEndX = isolated ? srcOffsetX + size.width : wholeSize.width;
1284     int srcEndY = isolated ? srcOffsetY + size.height : wholeSize.height;
1285     idxArg = kernel.set(idxArg, srcOffsetX);
1286     idxArg = kernel.set(idxArg, srcOffsetY);
1287     idxArg = kernel.set(idxArg, srcEndX);
1288     idxArg = kernel.set(idxArg, srcEndY);
1289     idxArg = kernel.set(idxArg, ocl::KernelArg::WriteOnly(dst));
1290     if (normalize)
1291         idxArg = kernel.set(idxArg, (float)alpha);
1292
1293     return kernel.run(2, globalsize, localsize, false);
1294 }
1295
1296 #undef ROUNDUP
1297
1298 #endif
1299
1300 }
1301
1302
1303 cv::Ptr<cv::BaseRowFilter> cv::getRowSumFilter(int srcType, int sumType, int ksize, int anchor)
1304 {
1305     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
1306     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
1307
1308     if( anchor < 0 )
1309         anchor = ksize/2;
1310
1311     if( sdepth == CV_8U && ddepth == CV_32S )
1312         return makePtr<RowSum<uchar, int> >(ksize, anchor);
1313     if( sdepth == CV_8U && ddepth == CV_16U )
1314         return makePtr<RowSum<uchar, ushort> >(ksize, anchor);
1315     if( sdepth == CV_8U && ddepth == CV_64F )
1316         return makePtr<RowSum<uchar, double> >(ksize, anchor);
1317     if( sdepth == CV_16U && ddepth == CV_32S )
1318         return makePtr<RowSum<ushort, int> >(ksize, anchor);
1319     if( sdepth == CV_16U && ddepth == CV_64F )
1320         return makePtr<RowSum<ushort, double> >(ksize, anchor);
1321     if( sdepth == CV_16S && ddepth == CV_32S )
1322         return makePtr<RowSum<short, int> >(ksize, anchor);
1323     if( sdepth == CV_32S && ddepth == CV_32S )
1324         return makePtr<RowSum<int, int> >(ksize, anchor);
1325     if( sdepth == CV_16S && ddepth == CV_64F )
1326         return makePtr<RowSum<short, double> >(ksize, anchor);
1327     if( sdepth == CV_32F && ddepth == CV_64F )
1328         return makePtr<RowSum<float, double> >(ksize, anchor);
1329     if( sdepth == CV_64F && ddepth == CV_64F )
1330         return makePtr<RowSum<double, double> >(ksize, anchor);
1331
1332     CV_Error_( CV_StsNotImplemented,
1333         ("Unsupported combination of source format (=%d), and buffer format (=%d)",
1334         srcType, sumType));
1335
1336     return Ptr<BaseRowFilter>();
1337 }
1338
1339
1340 cv::Ptr<cv::BaseColumnFilter> cv::getColumnSumFilter(int sumType, int dstType, int ksize,
1341                                                      int anchor, double scale)
1342 {
1343     int sdepth = CV_MAT_DEPTH(sumType), ddepth = CV_MAT_DEPTH(dstType);
1344     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(dstType) );
1345
1346     if( anchor < 0 )
1347         anchor = ksize/2;
1348
1349     if( ddepth == CV_8U && sdepth == CV_32S )
1350         return makePtr<ColumnSum<int, uchar> >(ksize, anchor, scale);
1351     if( ddepth == CV_8U && sdepth == CV_16U )
1352         return makePtr<ColumnSum<ushort, uchar> >(ksize, anchor, scale);
1353     if( ddepth == CV_8U && sdepth == CV_64F )
1354         return makePtr<ColumnSum<double, uchar> >(ksize, anchor, scale);
1355     if( ddepth == CV_16U && sdepth == CV_32S )
1356         return makePtr<ColumnSum<int, ushort> >(ksize, anchor, scale);
1357     if( ddepth == CV_16U && sdepth == CV_64F )
1358         return makePtr<ColumnSum<double, ushort> >(ksize, anchor, scale);
1359     if( ddepth == CV_16S && sdepth == CV_32S )
1360         return makePtr<ColumnSum<int, short> >(ksize, anchor, scale);
1361     if( ddepth == CV_16S && sdepth == CV_64F )
1362         return makePtr<ColumnSum<double, short> >(ksize, anchor, scale);
1363     if( ddepth == CV_32S && sdepth == CV_32S )
1364         return makePtr<ColumnSum<int, int> >(ksize, anchor, scale);
1365     if( ddepth == CV_32F && sdepth == CV_32S )
1366         return makePtr<ColumnSum<int, float> >(ksize, anchor, scale);
1367     if( ddepth == CV_32F && sdepth == CV_64F )
1368         return makePtr<ColumnSum<double, float> >(ksize, anchor, scale);
1369     if( ddepth == CV_64F && sdepth == CV_32S )
1370         return makePtr<ColumnSum<int, double> >(ksize, anchor, scale);
1371     if( ddepth == CV_64F && sdepth == CV_64F )
1372         return makePtr<ColumnSum<double, double> >(ksize, anchor, scale);
1373
1374     CV_Error_( CV_StsNotImplemented,
1375         ("Unsupported combination of sum format (=%d), and destination format (=%d)",
1376         sumType, dstType));
1377
1378     return Ptr<BaseColumnFilter>();
1379 }
1380
1381
1382 cv::Ptr<cv::FilterEngine> cv::createBoxFilter( int srcType, int dstType, Size ksize,
1383                     Point anchor, bool normalize, int borderType )
1384 {
1385     int sdepth = CV_MAT_DEPTH(srcType);
1386     int cn = CV_MAT_CN(srcType), sumType = CV_64F;
1387     if( sdepth == CV_8U && CV_MAT_DEPTH(dstType) == CV_8U &&
1388         ksize.width*ksize.height <= 256 )
1389         sumType = CV_16U;
1390     else if( sdepth <= CV_32S && (!normalize ||
1391         ksize.width*ksize.height <= (sdepth == CV_8U ? (1<<23) :
1392             sdepth == CV_16U ? (1 << 15) : (1 << 16))) )
1393         sumType = CV_32S;
1394     sumType = CV_MAKETYPE( sumType, cn );
1395
1396     Ptr<BaseRowFilter> rowFilter = getRowSumFilter(srcType, sumType, ksize.width, anchor.x );
1397     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
1398         dstType, ksize.height, anchor.y, normalize ? 1./(ksize.width*ksize.height) : 1);
1399
1400     return makePtr<FilterEngine>(Ptr<BaseFilter>(), rowFilter, columnFilter,
1401            srcType, dstType, sumType, borderType );
1402 }
1403
1404 #ifdef HAVE_OPENVX
1405 namespace cv
1406 {
1407     namespace ovx {
1408         template <> inline bool skipSmallImages<VX_KERNEL_BOX_3x3>(int w, int h) { return w*h < 640 * 480; }
1409     }
1410     static bool openvx_boxfilter(InputArray _src, OutputArray _dst, int ddepth,
1411                                  Size ksize, Point anchor,
1412                                  bool normalize, int borderType)
1413     {
1414         if (ddepth < 0)
1415             ddepth = CV_8UC1;
1416         if (_src.type() != CV_8UC1 || ddepth != CV_8U || !normalize ||
1417             _src.cols() < 3 || _src.rows() < 3 ||
1418             ksize.width != 3 || ksize.height != 3 ||
1419             (anchor.x >= 0 && anchor.x != 1) ||
1420             (anchor.y >= 0 && anchor.y != 1) ||
1421             ovx::skipSmallImages<VX_KERNEL_BOX_3x3>(_src.cols(), _src.rows()))
1422             return false;
1423
1424         Mat src = _src.getMat();
1425
1426         if ((borderType & BORDER_ISOLATED) == 0 && src.isSubmatrix())
1427             return false; //Process isolated borders only
1428         vx_enum border;
1429         switch (borderType & ~BORDER_ISOLATED)
1430         {
1431         case BORDER_CONSTANT:
1432             border = VX_BORDER_CONSTANT;
1433             break;
1434         case BORDER_REPLICATE:
1435             border = VX_BORDER_REPLICATE;
1436             break;
1437         default:
1438             return false;
1439         }
1440
1441         _dst.create(src.size(), CV_8UC1);
1442         Mat dst = _dst.getMat();
1443
1444         try
1445         {
1446             ivx::Context ctx = ovx::getOpenVXContext();
1447
1448             Mat a;
1449             if (dst.data != src.data)
1450                 a = src;
1451             else
1452                 src.copyTo(a);
1453
1454             ivx::Image
1455                 ia = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
1456                                                   ivx::Image::createAddressing(a.cols, a.rows, 1, (vx_int32)(a.step)), a.data),
1457                 ib = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
1458                                                   ivx::Image::createAddressing(dst.cols, dst.rows, 1, (vx_int32)(dst.step)), dst.data);
1459
1460             //ATTENTION: VX_CONTEXT_IMMEDIATE_BORDER attribute change could lead to strange issues in multi-threaded environments
1461             //since OpenVX standard says nothing about thread-safety for now
1462             ivx::border_t prevBorder = ctx.immediateBorder();
1463             ctx.setImmediateBorder(border, (vx_uint8)(0));
1464             ivx::IVX_CHECK_STATUS(vxuBox3x3(ctx, ia, ib));
1465             ctx.setImmediateBorder(prevBorder);
1466         }
1467         catch (ivx::RuntimeError & e)
1468         {
1469             VX_DbgThrow(e.what());
1470         }
1471         catch (ivx::WrapperError & e)
1472         {
1473             VX_DbgThrow(e.what());
1474         }
1475
1476         return true;
1477     }
1478 }
1479 #endif
1480
1481 #if defined(HAVE_IPP)
1482 namespace cv
1483 {
1484 static bool ipp_boxfilter(Mat &src, Mat &dst, Size ksize, Point anchor, bool normalize, int borderType)
1485 {
1486 #ifdef HAVE_IPP_IW
1487     CV_INSTRUMENT_REGION_IPP()
1488
1489 #if IPP_VERSION_X100 < 201801
1490     // Problem with SSE42 optimization for 16s and some 8u modes
1491     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))))
1492         return false;
1493
1494     // Other optimizations has some degradations too
1495     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))))
1496         return false;
1497 #endif
1498
1499     if(!normalize)
1500         return false;
1501
1502     if(!ippiCheckAnchor(anchor, ksize))
1503         return false;
1504
1505     try
1506     {
1507         ::ipp::IwiImage       iwSrc      = ippiGetImage(src);
1508         ::ipp::IwiImage       iwDst      = ippiGetImage(dst);
1509         ::ipp::IwiSize        iwKSize    = ippiGetSize(ksize);
1510         ::ipp::IwiBorderSize  borderSize(iwKSize);
1511         ::ipp::IwiBorderType  ippBorder(ippiGetBorder(iwSrc, borderType, borderSize));
1512         if(!ippBorder)
1513             return false;
1514
1515         CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBox, iwSrc, iwDst, iwKSize, ::ipp::IwDefault(), ippBorder);
1516     }
1517     catch (::ipp::IwException)
1518     {
1519         return false;
1520     }
1521
1522     return true;
1523 #else
1524     CV_UNUSED(src); CV_UNUSED(dst); CV_UNUSED(ksize); CV_UNUSED(anchor); CV_UNUSED(normalize); CV_UNUSED(borderType);
1525     return false;
1526 #endif
1527 }
1528 }
1529 #endif
1530
1531
1532 void cv::boxFilter( InputArray _src, OutputArray _dst, int ddepth,
1533                 Size ksize, Point anchor,
1534                 bool normalize, int borderType )
1535 {
1536     CV_INSTRUMENT_REGION()
1537
1538     CV_OCL_RUN(_dst.isUMat() &&
1539                (borderType == BORDER_REPLICATE || borderType == BORDER_CONSTANT ||
1540                 borderType == BORDER_REFLECT || borderType == BORDER_REFLECT_101),
1541                ocl_boxFilter3x3_8UC1(_src, _dst, ddepth, ksize, anchor, borderType, normalize))
1542
1543     CV_OCL_RUN(_dst.isUMat(), ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize))
1544
1545     Mat src = _src.getMat();
1546     int stype = src.type(), sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype);
1547     if( ddepth < 0 )
1548         ddepth = sdepth;
1549     _dst.create( src.size(), CV_MAKETYPE(ddepth, cn) );
1550     Mat dst = _dst.getMat();
1551     if( borderType != BORDER_CONSTANT && normalize && (borderType & BORDER_ISOLATED) != 0 )
1552     {
1553         if( src.rows == 1 )
1554             ksize.height = 1;
1555         if( src.cols == 1 )
1556             ksize.width = 1;
1557     }
1558
1559     Point ofs;
1560     Size wsz(src.cols, src.rows);
1561     if(!(borderType&BORDER_ISOLATED))
1562         src.locateROI( wsz, ofs );
1563
1564     CALL_HAL(boxFilter, cv_hal_boxFilter, src.ptr(), src.step, dst.ptr(), dst.step, src.cols, src.rows, sdepth, ddepth, cn,
1565              ofs.x, ofs.y, wsz.width - src.cols - ofs.x, wsz.height - src.rows - ofs.y, ksize.width, ksize.height,
1566              anchor.x, anchor.y, normalize, borderType&~BORDER_ISOLATED);
1567
1568     CV_OVX_RUN(true,
1569                openvx_boxfilter(src, dst, ddepth, ksize, anchor, normalize, borderType))
1570
1571     CV_IPP_RUN_FAST(ipp_boxfilter(src, dst, ksize, anchor, normalize, borderType));
1572
1573     borderType = (borderType&~BORDER_ISOLATED);
1574
1575     Ptr<FilterEngine> f = createBoxFilter( src.type(), dst.type(),
1576                         ksize, anchor, normalize, borderType );
1577
1578     f->apply( src, dst, wsz, ofs );
1579 }
1580
1581
1582 void cv::blur( InputArray src, OutputArray dst,
1583            Size ksize, Point anchor, int borderType )
1584 {
1585     CV_INSTRUMENT_REGION()
1586
1587     boxFilter( src, dst, -1, ksize, anchor, true, borderType );
1588 }
1589
1590
1591 /****************************************************************************************\
1592                                     Squared Box Filter
1593 \****************************************************************************************/
1594
1595 namespace cv
1596 {
1597
1598 template<typename T, typename ST>
1599 struct SqrRowSum :
1600         public BaseRowFilter
1601 {
1602     SqrRowSum( int _ksize, int _anchor ) :
1603         BaseRowFilter()
1604     {
1605         ksize = _ksize;
1606         anchor = _anchor;
1607     }
1608
1609     virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
1610     {
1611         const T* S = (const T*)src;
1612         ST* D = (ST*)dst;
1613         int i = 0, k, ksz_cn = ksize*cn;
1614
1615         width = (width - 1)*cn;
1616         for( k = 0; k < cn; k++, S++, D++ )
1617         {
1618             ST s = 0;
1619             for( i = 0; i < ksz_cn; i += cn )
1620             {
1621                 ST val = (ST)S[i];
1622                 s += val*val;
1623             }
1624             D[0] = s;
1625             for( i = 0; i < width; i += cn )
1626             {
1627                 ST val0 = (ST)S[i], val1 = (ST)S[i + ksz_cn];
1628                 s += val1*val1 - val0*val0;
1629                 D[i+cn] = s;
1630             }
1631         }
1632     }
1633 };
1634
1635 static Ptr<BaseRowFilter> getSqrRowSumFilter(int srcType, int sumType, int ksize, int anchor)
1636 {
1637     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
1638     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
1639
1640     if( anchor < 0 )
1641         anchor = ksize/2;
1642
1643     if( sdepth == CV_8U && ddepth == CV_32S )
1644         return makePtr<SqrRowSum<uchar, int> >(ksize, anchor);
1645     if( sdepth == CV_8U && ddepth == CV_64F )
1646         return makePtr<SqrRowSum<uchar, double> >(ksize, anchor);
1647     if( sdepth == CV_16U && ddepth == CV_64F )
1648         return makePtr<SqrRowSum<ushort, double> >(ksize, anchor);
1649     if( sdepth == CV_16S && ddepth == CV_64F )
1650         return makePtr<SqrRowSum<short, double> >(ksize, anchor);
1651     if( sdepth == CV_32F && ddepth == CV_64F )
1652         return makePtr<SqrRowSum<float, double> >(ksize, anchor);
1653     if( sdepth == CV_64F && ddepth == CV_64F )
1654         return makePtr<SqrRowSum<double, double> >(ksize, anchor);
1655
1656     CV_Error_( CV_StsNotImplemented,
1657               ("Unsupported combination of source format (=%d), and buffer format (=%d)",
1658                srcType, sumType));
1659
1660     return Ptr<BaseRowFilter>();
1661 }
1662
1663 }
1664
1665 void cv::sqrBoxFilter( InputArray _src, OutputArray _dst, int ddepth,
1666                        Size ksize, Point anchor,
1667                        bool normalize, int borderType )
1668 {
1669     CV_INSTRUMENT_REGION()
1670
1671     int srcType = _src.type(), sdepth = CV_MAT_DEPTH(srcType), cn = CV_MAT_CN(srcType);
1672     Size size = _src.size();
1673
1674     if( ddepth < 0 )
1675         ddepth = sdepth < CV_32F ? CV_32F : CV_64F;
1676
1677     if( borderType != BORDER_CONSTANT && normalize )
1678     {
1679         if( size.height == 1 )
1680             ksize.height = 1;
1681         if( size.width == 1 )
1682             ksize.width = 1;
1683     }
1684
1685     CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2,
1686                ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize, true))
1687
1688     int sumDepth = CV_64F;
1689     if( sdepth == CV_8U )
1690         sumDepth = CV_32S;
1691     int sumType = CV_MAKETYPE( sumDepth, cn ), dstType = CV_MAKETYPE(ddepth, cn);
1692
1693     Mat src = _src.getMat();
1694     _dst.create( size, dstType );
1695     Mat dst = _dst.getMat();
1696
1697     Ptr<BaseRowFilter> rowFilter = getSqrRowSumFilter(srcType, sumType, ksize.width, anchor.x );
1698     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
1699                                                             dstType, ksize.height, anchor.y,
1700                                                             normalize ? 1./(ksize.width*ksize.height) : 1);
1701
1702     Ptr<FilterEngine> f = makePtr<FilterEngine>(Ptr<BaseFilter>(), rowFilter, columnFilter,
1703                                                 srcType, dstType, sumType, borderType );
1704     Point ofs;
1705     Size wsz(src.cols, src.rows);
1706     src.locateROI( wsz, ofs );
1707
1708     f->apply( src, dst, wsz, ofs );
1709 }
1710
1711
1712 /****************************************************************************************\
1713                                      Gaussian Blur
1714 \****************************************************************************************/
1715
1716 cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )
1717 {
1718     const int SMALL_GAUSSIAN_SIZE = 7;
1719     static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
1720     {
1721         {1.f},
1722         {0.25f, 0.5f, 0.25f},
1723         {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
1724         {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
1725     };
1726
1727     const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
1728         small_gaussian_tab[n>>1] : 0;
1729
1730     CV_Assert( ktype == CV_32F || ktype == CV_64F );
1731     Mat kernel(n, 1, ktype);
1732     float* cf = kernel.ptr<float>();
1733     double* cd = kernel.ptr<double>();
1734
1735     double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
1736     double scale2X = -0.5/(sigmaX*sigmaX);
1737     double sum = 0;
1738
1739     int i;
1740     for( i = 0; i < n; i++ )
1741     {
1742         double x = i - (n-1)*0.5;
1743         double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);
1744         if( ktype == CV_32F )
1745         {
1746             cf[i] = (float)t;
1747             sum += cf[i];
1748         }
1749         else
1750         {
1751             cd[i] = t;
1752             sum += cd[i];
1753         }
1754     }
1755
1756     sum = 1./sum;
1757     for( i = 0; i < n; i++ )
1758     {
1759         if( ktype == CV_32F )
1760             cf[i] = (float)(cf[i]*sum);
1761         else
1762             cd[i] *= sum;
1763     }
1764
1765     return kernel;
1766 }
1767
1768 namespace cv {
1769
1770 template <typename T>
1771 static std::vector<T> getFixedpointGaussianKernel( int n, double sigma )
1772 {
1773     if (sigma <= 0)
1774     {
1775         if(n == 1)
1776             return std::vector<T>(1, softdouble(1.0));
1777         else if(n == 3)
1778         {
1779             T v3[] = { softdouble(0.25), softdouble(0.5), softdouble(0.25) };
1780             return std::vector<T>(v3, v3 + 3);
1781         }
1782         else if(n == 5)
1783         {
1784             T v5[] = { softdouble(0.0625), softdouble(0.25), softdouble(0.375), softdouble(0.25), softdouble(0.0625) };
1785             return std::vector<T>(v5, v5 + 5);
1786         }
1787         else if(n == 7)
1788         {
1789             T v7[] = { softdouble(0.03125), softdouble(0.109375), softdouble(0.21875), softdouble(0.28125), softdouble(0.21875), softdouble(0.109375), softdouble(0.03125) };
1790             return std::vector<T>(v7, v7 + 7);
1791         }
1792     }
1793
1794
1795     softdouble sigmaX = sigma > 0 ? softdouble(sigma) : mulAdd(softdouble(n),softdouble(0.15),softdouble(0.35));// softdouble(((n-1)*0.5 - 1)*0.3 + 0.8)
1796     softdouble scale2X = softdouble(-0.5*0.25)/(sigmaX*sigmaX);
1797     std::vector<softdouble> values(n);
1798     softdouble sum(0.);
1799     for(int i = 0, x = 1 - n; i < n; i++, x+=2 )
1800     {
1801         // x = i - (n - 1)*0.5
1802         // t = std::exp(scale2X*x*x)
1803         values[i] = exp(softdouble(x*x)*scale2X);
1804         sum += values[i];
1805     }
1806     sum = softdouble::one()/sum;
1807
1808     std::vector<T> kernel(n);
1809     for(int i = 0; i < n; i++ )
1810     {
1811         kernel[i] = values[i] * sum;
1812     }
1813
1814     return kernel;
1815 };
1816
1817 template <typename ET, typename FT>
1818 void hlineSmooth1N(const ET* src, int cn, const FT* m, int, FT* dst, int len, int)
1819 {
1820     for (int i = 0; i < len*cn; i++, src++, dst++)
1821         *dst = (*m) * (*src);
1822 }
1823 template <>
1824 void hlineSmooth1N<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16* m, int, ufixedpoint16* dst, int len, int)
1825 {
1826     int lencn = len*cn;
1827     v_uint16x8 v_mul = v_setall_u16(*((uint16_t*)m));
1828     int i = 0;
1829     for (; i < lencn - 15; i += 16)
1830     {
1831         v_uint8x16 v_src = v_load(src + i);
1832         v_uint16x8 v_tmp0, v_tmp1;
1833         v_expand(v_src, v_tmp0, v_tmp1);
1834         v_store((uint16_t*)dst + i, v_mul*v_tmp0);
1835         v_store((uint16_t*)dst + i + 8, v_mul*v_tmp1);
1836     }
1837     if (i < lencn - 7)
1838     {
1839         v_uint16x8 v_src = v_load_expand(src + i);
1840         v_store((uint16_t*)dst + i, v_mul*v_src);
1841         i += 8;
1842     }
1843     for (; i < lencn; i++)
1844         dst[i] = m[0] * src[i];
1845 }
1846 template <typename ET, typename FT>
1847 void hlineSmooth1N1(const ET* src, int cn, const FT*, int, FT* dst, int len, int)
1848 {
1849     for (int i = 0; i < len*cn; i++, src++, dst++)
1850         *dst = *src;
1851 }
1852 template <>
1853 void hlineSmooth1N1<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16*, int, ufixedpoint16* dst, int len, int)
1854 {
1855     int lencn = len*cn;
1856     int i = 0;
1857     for (; i < lencn - 15; i += 16)
1858     {
1859         v_uint8x16 v_src = v_load(src + i);
1860         v_uint16x8 v_tmp0, v_tmp1;
1861         v_expand(v_src, v_tmp0, v_tmp1);
1862         v_store((uint16_t*)dst + i, v_shl<8>(v_tmp0));
1863         v_store((uint16_t*)dst + i + 8, v_shl<8>(v_tmp1));
1864     }
1865     if (i < lencn - 7)
1866     {
1867         v_uint16x8 v_src = v_load_expand(src + i);
1868         v_store((uint16_t*)dst + i, v_shl<8>(v_src));
1869         i += 8;
1870     }
1871     for (; i < lencn; i++)
1872         dst[i] = src[i];
1873 }
1874 template <typename ET, typename FT>
1875 void hlineSmooth3N(const ET* src, int cn, const FT* m, int, FT* dst, int len, int borderType)
1876 {
1877     if (len == 1)
1878     {
1879         FT msum = borderType != BORDER_CONSTANT ? m[0] + m[1] + m[2] : m[1];
1880         for (int k = 0; k < cn; k++)
1881             dst[k] = msum * src[k];
1882     }
1883     else
1884     {
1885         // Point that fall left from border
1886         for (int k = 0; k < cn; k++)
1887             dst[k] = m[1] * src[k] + m[2] * src[cn + k];
1888         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
1889         {
1890             int src_idx = borderInterpolate(-1, len, borderType);
1891             for (int k = 0; k < cn; k++)
1892                 dst[k] = dst[k] + m[0] * src[src_idx*cn + k];
1893         }
1894
1895         src += cn; dst += cn;
1896         for (int i = cn; i < (len - 1)*cn; i++, src++, dst++)
1897             *dst = m[0] * src[-cn] + m[1] * src[0] + m[2] * src[cn];
1898
1899         // Point that fall right from border
1900         for (int k = 0; k < cn; k++)
1901             dst[k] = m[0] * src[k - cn] + m[1] * src[k];
1902         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
1903         {
1904             int src_idx = (borderInterpolate(len, len, borderType) - (len - 1))*cn;
1905             for (int k = 0; k < cn; k++)
1906                 dst[k] = dst[k] + m[2] * src[src_idx + k];
1907         }
1908     }
1909 }
1910 template <>
1911 void hlineSmooth3N<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16* m, int, ufixedpoint16* dst, int len, int borderType)
1912 {
1913     if (len == 1)
1914     {
1915         ufixedpoint16 msum = borderType != BORDER_CONSTANT ? m[0] + m[1] + m[2] : m[1];
1916         for (int k = 0; k < cn; k++)
1917             dst[k] = msum * src[k];
1918     }
1919     else
1920     {
1921         // Point that fall left from border
1922         for (int k = 0; k < cn; k++)
1923             dst[k] = m[1] * src[k] + m[2] * src[cn + k];
1924         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
1925         {
1926             int src_idx = borderInterpolate(-1, len, borderType);
1927             for (int k = 0; k < cn; k++)
1928                 dst[k] = dst[k] + m[0] * src[src_idx*cn + k];
1929         }
1930
1931         src += cn; dst += cn;
1932         int i = cn, lencn = (len - 1)*cn;
1933         v_int16x8 v_mul01 = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)m)));
1934         v_int16x8 v_mul2 = v_reinterpret_as_s16(v_setall_u16(*((uint16_t*)(m + 2))));
1935         for (; i < lencn - 15; i += 16, src += 16, dst += 16)
1936         {
1937             v_uint16x8 v_src00, v_src01, v_src10, v_src11;
1938             v_int16x8 v_tmp0, v_tmp1;
1939
1940             v_expand(v_load(src - cn), v_src00, v_src01);
1941             v_expand(v_load(src), v_src10, v_src11);
1942             v_zip(v_reinterpret_as_s16(v_src00), v_reinterpret_as_s16(v_src10), v_tmp0, v_tmp1);
1943             v_int32x4 v_res0 = v_dotprod(v_tmp0, v_mul01);
1944             v_int32x4 v_res1 = v_dotprod(v_tmp1, v_mul01);
1945             v_zip(v_reinterpret_as_s16(v_src01), v_reinterpret_as_s16(v_src11), v_tmp0, v_tmp1);
1946             v_int32x4 v_res2 = v_dotprod(v_tmp0, v_mul01);
1947             v_int32x4 v_res3 = v_dotprod(v_tmp1, v_mul01);
1948
1949             v_int32x4 v_resj0, v_resj1, v_resj2, v_resj3;
1950             v_expand(v_load(src + cn), v_src00, v_src01);
1951             v_mul_expand(v_reinterpret_as_s16(v_src00), v_mul2, v_resj0, v_resj1);
1952             v_mul_expand(v_reinterpret_as_s16(v_src01), v_mul2, v_resj2, v_resj3);
1953             v_res0 += v_resj0;
1954             v_res1 += v_resj1;
1955             v_res2 += v_resj2;
1956             v_res3 += v_resj3;
1957
1958             v_store((uint16_t*)dst, v_pack(v_reinterpret_as_u32(v_res0), v_reinterpret_as_u32(v_res1)));
1959             v_store((uint16_t*)dst + 8, v_pack(v_reinterpret_as_u32(v_res2), v_reinterpret_as_u32(v_res3)));
1960         }
1961         for (; i < lencn; i++, src++, dst++)
1962             *dst = m[0] * src[-cn] + m[1] * src[0] + m[2] * src[cn];
1963
1964         // Point that fall right from border
1965         for (int k = 0; k < cn; k++)
1966             dst[k] = m[0] * src[k - cn] + m[1] * src[k];
1967         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
1968         {
1969             int src_idx = (borderInterpolate(len, len, borderType) - (len - 1))*cn;
1970             for (int k = 0; k < cn; k++)
1971                 dst[k] = dst[k] + m[2] * src[src_idx + k];
1972         }
1973     }
1974 }
1975 template <typename ET, typename FT>
1976 void hlineSmooth3N121(const ET* src, int cn, const FT*, int, FT* dst, int len, int borderType)
1977 {
1978     if (len == 1)
1979     {
1980         if(borderType != BORDER_CONSTANT)
1981             for (int k = 0; k < cn; k++)
1982                 dst[k] = FT(src[k]);
1983         else
1984             for (int k = 0; k < cn; k++)
1985                 dst[k] = FT(src[k])>>1;
1986     }
1987     else
1988     {
1989         // Point that fall left from border
1990         for (int k = 0; k < cn; k++)
1991             dst[k] = (FT(src[k])>>1) + (FT(src[cn + k])>>2);
1992         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
1993         {
1994             int src_idx = borderInterpolate(-1, len, borderType);
1995             for (int k = 0; k < cn; k++)
1996                 dst[k] = dst[k] + (FT(src[src_idx*cn + k])>>2);
1997         }
1998
1999         src += cn; dst += cn;
2000         for (int i = cn; i < (len - 1)*cn; i++, src++, dst++)
2001             *dst = ((FT(src[-cn]) + FT(src[cn]))>>2) + (FT(src[0])>>1);
2002
2003         // Point that fall right from border
2004         for (int k = 0; k < cn; k++)
2005             dst[k] = (FT(src[k - cn])>>2) + (FT(src[k])>>1);
2006         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2007         {
2008             int src_idx = (borderInterpolate(len, len, borderType) - (len - 1))*cn;
2009             for (int k = 0; k < cn; k++)
2010                 dst[k] = dst[k] + (FT(src[src_idx + k])>>2);
2011         }
2012     }
2013 }
2014 template <>
2015 void hlineSmooth3N121<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16*, int, ufixedpoint16* dst, int len, int borderType)
2016 {
2017     if (len == 1)
2018     {
2019         if (borderType != BORDER_CONSTANT)
2020             for (int k = 0; k < cn; k++)
2021                 dst[k] = ufixedpoint16(src[k]);
2022         else
2023             for (int k = 0; k < cn; k++)
2024                 dst[k] = ufixedpoint16(src[k]) >> 1;
2025     }
2026     else
2027     {
2028         // Point that fall left from border
2029         for (int k = 0; k < cn; k++)
2030             dst[k] = (ufixedpoint16(src[k])>>1) + (ufixedpoint16(src[cn + k])>>2);
2031         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2032         {
2033             int src_idx = borderInterpolate(-1, len, borderType);
2034             for (int k = 0; k < cn; k++)
2035                 dst[k] = dst[k] + (ufixedpoint16(src[src_idx*cn + k])>>2);
2036         }
2037
2038         src += cn; dst += cn;
2039         int i = cn, lencn = (len - 1)*cn;
2040         for (; i < lencn - 15; i += 16, src += 16, dst += 16)
2041         {
2042             v_uint16x8 v_src00, v_src01, v_src10, v_src11, v_src20, v_src21;
2043             v_expand(v_load(src - cn), v_src00, v_src01);
2044             v_expand(v_load(src), v_src10, v_src11);
2045             v_expand(v_load(src + cn), v_src20, v_src21);
2046             v_store((uint16_t*)dst, (v_src00 + v_src20 + (v_src10 << 1)) << 6);
2047             v_store((uint16_t*)dst + 8, (v_src01 + v_src21 + (v_src11 << 1)) << 6);
2048         }
2049         for (; i < lencn; i++, src++, dst++)
2050             *((uint16_t*)dst) = (uint16_t(src[-cn]) + uint16_t(src[cn]) + (uint16_t(src[0]) << 1)) << 6;
2051
2052         // Point that fall right from border
2053         for (int k = 0; k < cn; k++)
2054             dst[k] = (ufixedpoint16(src[k - cn])>>2) + (ufixedpoint16(src[k])>>1);
2055         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2056         {
2057             int src_idx = (borderInterpolate(len, len, borderType) - (len - 1))*cn;
2058             for (int k = 0; k < cn; k++)
2059                 dst[k] = dst[k] + (ufixedpoint16(src[src_idx + k])>>2);
2060         }
2061     }
2062 }
2063 template <typename ET, typename FT>
2064 void hlineSmooth5N(const ET* src, int cn, const FT* m, int, FT* dst, int len, int borderType)
2065 {
2066     if (len == 1)
2067     {
2068         ufixedpoint16 msum = borderType != BORDER_CONSTANT ? m[0] + m[1] + m[2] + m[3] + m[4] : m[2];
2069         for (int k = 0; k < cn; k++)
2070             dst[k] = msum * src[k];
2071     }
2072     else if (len == 2)
2073     {
2074         if (borderType == BORDER_CONSTANT)
2075             for (int k = 0; k < cn; k++)
2076             {
2077                 dst[k   ] = m[2] * src[k] + m[3] * src[k+cn];
2078                 dst[k+cn] = m[1] * src[k] + m[2] * src[k+cn];
2079             }
2080         else
2081         {
2082             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2083             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2084             int idxp1 = borderInterpolate(2, len, borderType)*cn;
2085             int idxp2 = borderInterpolate(3, len, borderType)*cn;
2086             for (int k = 0; k < cn; k++)
2087             {
2088                 dst[k     ] = m[1] * src[k + idxm1] + m[2] * src[k] + m[3] * src[k + cn] + m[4] * src[k + idxp1] + m[0] * src[k + idxm2];
2089                 dst[k + cn] = m[0] * src[k + idxm1] + m[1] * src[k] + m[2] * src[k + cn] + m[3] * src[k + idxp1] + m[4] * src[k + idxp2];
2090             }
2091         }
2092     }
2093     else if (len == 3)
2094     {
2095         if (borderType == BORDER_CONSTANT)
2096             for (int k = 0; k < cn; k++)
2097             {
2098                 dst[k       ] = m[2] * src[k] + m[3] * src[k + cn] + m[4] * src[k + 2*cn];
2099                 dst[k +   cn] = m[1] * src[k] + m[2] * src[k + cn] + m[3] * src[k + 2*cn];
2100                 dst[k + 2*cn] = m[0] * src[k] + m[1] * src[k + cn] + m[2] * src[k + 2*cn];
2101             }
2102         else
2103         {
2104             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2105             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2106             int idxp1 = borderInterpolate(3, len, borderType)*cn;
2107             int idxp2 = borderInterpolate(4, len, borderType)*cn;
2108             for (int k = 0; k < cn; k++)
2109             {
2110                 dst[k       ] = m[2] * src[k] + m[3] * src[k + cn] + m[4] * src[k + 2*cn] + m[0] * src[k + idxm2] + m[1] * src[k + idxm1];
2111                 dst[k +   cn] = m[1] * src[k] + m[2] * src[k + cn] + m[3] * src[k + 2*cn] + m[0] * src[k + idxm1] + m[4] * src[k + idxp1];
2112                 dst[k + 2*cn] = m[0] * src[k] + m[1] * src[k + cn] + m[2] * src[k + 2*cn] + m[3] * src[k + idxp1] + m[4] * src[k + idxp2];
2113             }
2114         }
2115     }
2116     else
2117     {
2118         // Points that fall left from border
2119         for (int k = 0; k < cn; k++)
2120         {
2121             dst[k] = m[2] * src[k] + m[3] * src[cn + k] + m[4] * src[2*cn + k];
2122             dst[k + cn] = m[1] * src[k] + m[2] * src[cn + k] + m[3] * src[2*cn + k] + m[4] * src[3*cn + k];
2123         }
2124         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2125         {
2126             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2127             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2128             for (int k = 0; k < cn; k++)
2129             {
2130                 dst[k] = dst[k] + m[0] * src[idxm2 + k] + m[1] * src[idxm1 + k];
2131                 dst[k + cn] = dst[k + cn] + m[0] * src[idxm1 + k];
2132             }
2133         }
2134
2135         src += 2*cn; dst += 2*cn;
2136         for (int i = 2*cn; i < (len - 2)*cn; i++, src++, dst++)
2137             *dst = m[0] * src[-2*cn] + m[1] * src[-cn] + m[2] * src[0] + m[3] * src[cn] + m[4] * src[2*cn];
2138
2139         // Points that fall right from border
2140         for (int k = 0; k < cn; k++)
2141         {
2142             dst[k] = m[0] * src[k - 2*cn] + m[1] * src[k - cn] + m[2] * src[k] + m[3] * src[k + cn];
2143             dst[k + cn] = m[0] * src[k - cn] + m[1] * src[k] + m[2] * src[k + cn];
2144         }
2145         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2146         {
2147             int idxp1 = (borderInterpolate(len, len, borderType) - (len - 2))*cn;
2148             int idxp2 = (borderInterpolate(len+1, len, borderType) - (len - 2))*cn;
2149             for (int k = 0; k < cn; k++)
2150             {
2151                 dst[k] = dst[k] + m[4] * src[idxp1 + k];
2152                 dst[k + cn] = dst[k + cn] + m[3] * src[idxp1 + k] + m[4] * src[idxp2 + k];
2153             }
2154         }
2155     }
2156 }
2157 template <>
2158 void hlineSmooth5N<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16* m, int, ufixedpoint16* dst, int len, int borderType)
2159 {
2160     if (len == 1)
2161     {
2162         ufixedpoint16 msum = borderType != BORDER_CONSTANT ? m[0] + m[1] + m[2] + m[3] + m[4] : m[2];
2163         for (int k = 0; k < cn; k++)
2164             dst[k] = msum * src[k];
2165     }
2166     else if (len == 2)
2167     {
2168         if (borderType == BORDER_CONSTANT)
2169             for (int k = 0; k < cn; k++)
2170             {
2171                 dst[k] = m[2] * src[k] + m[3] * src[k + cn];
2172                 dst[k + cn] = m[1] * src[k] + m[2] * src[k + cn];
2173             }
2174         else
2175         {
2176             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2177             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2178             int idxp1 = borderInterpolate(2, len, borderType)*cn;
2179             int idxp2 = borderInterpolate(3, len, borderType)*cn;
2180             for (int k = 0; k < cn; k++)
2181             {
2182                 dst[k] = m[1] * src[k + idxm1] + m[2] * src[k] + m[3] * src[k + cn] + m[4] * src[k + idxp1] + m[0] * src[k + idxm2];
2183                 dst[k + cn] = m[0] * src[k + idxm1] + m[1] * src[k] + m[2] * src[k + cn] + m[3] * src[k + idxp1] + m[4] * src[k + idxp2];
2184             }
2185         }
2186     }
2187     else if (len == 3)
2188     {
2189         if (borderType == BORDER_CONSTANT)
2190             for (int k = 0; k < cn; k++)
2191             {
2192                 dst[k] = m[2] * src[k] + m[3] * src[k + cn] + m[4] * src[k + 2 * cn];
2193                 dst[k + cn] = m[1] * src[k] + m[2] * src[k + cn] + m[3] * src[k + 2 * cn];
2194                 dst[k + 2 * cn] = m[0] * src[k] + m[1] * src[k + cn] + m[2] * src[k + 2 * cn];
2195             }
2196         else
2197         {
2198             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2199             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2200             int idxp1 = borderInterpolate(3, len, borderType)*cn;
2201             int idxp2 = borderInterpolate(4, len, borderType)*cn;
2202             for (int k = 0; k < cn; k++)
2203             {
2204                 dst[k] = m[2] * src[k] + m[3] * src[k + cn] + m[4] * src[k + 2 * cn] + m[0] * src[k + idxm2] + m[1] * src[k + idxm1];
2205                 dst[k + cn] = m[1] * src[k] + m[2] * src[k + cn] + m[3] * src[k + 2 * cn] + m[0] * src[k + idxm1] + m[4] * src[k + idxp1];
2206                 dst[k + 2 * cn] = m[0] * src[k] + m[1] * src[k + cn] + m[2] * src[k + 2 * cn] + m[3] * src[k + idxp1] + m[4] * src[k + idxp2];
2207             }
2208         }
2209     }
2210     else
2211     {
2212         // Points that fall left from border
2213         for (int k = 0; k < cn; k++)
2214         {
2215             dst[k] = m[2] * src[k] + m[3] * src[cn + k] + m[4] * src[2 * cn + k];
2216             dst[k + cn] = m[1] * src[k] + m[2] * src[cn + k] + m[3] * src[2 * cn + k] + m[4] * src[3 * cn + k];
2217         }
2218         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2219         {
2220             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2221             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2222             for (int k = 0; k < cn; k++)
2223             {
2224                 dst[k] = dst[k] + m[0] * src[idxm2 + k] + m[1] * src[idxm1 + k];
2225                 dst[k + cn] = dst[k + cn] + m[0] * src[idxm1 + k];
2226             }
2227         }
2228
2229         src += 2 * cn; dst += 2 * cn;
2230         int i = 2*cn, lencn = (len - 2)*cn;
2231         v_int16x8 v_mul01 = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)m)));
2232         v_int16x8 v_mul23 = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)(m + 2))));
2233         v_int16x8 v_mul4 = v_reinterpret_as_s16(v_setall_u16(*((uint16_t*)(m + 4))));
2234         for (; i < lencn - 15; i += 16, src += 16, dst += 16)
2235         {
2236             v_uint16x8 v_src00, v_src01, v_src10, v_src11;
2237             v_int16x8 v_tmp0, v_tmp1;
2238
2239             v_expand(v_load(src - 2*cn), v_src00, v_src01);
2240             v_expand(v_load(src - cn), v_src10, v_src11);
2241             v_zip(v_reinterpret_as_s16(v_src00), v_reinterpret_as_s16(v_src10), v_tmp0, v_tmp1);
2242             v_int32x4 v_res0 = v_dotprod(v_tmp0, v_mul01);
2243             v_int32x4 v_res1 = v_dotprod(v_tmp1, v_mul01);
2244             v_zip(v_reinterpret_as_s16(v_src01), v_reinterpret_as_s16(v_src11), v_tmp0, v_tmp1);
2245             v_int32x4 v_res2 = v_dotprod(v_tmp0, v_mul01);
2246             v_int32x4 v_res3 = v_dotprod(v_tmp1, v_mul01);
2247
2248
2249             v_expand(v_load(src), v_src00, v_src01);
2250             v_expand(v_load(src + cn), v_src10, v_src11);
2251             v_zip(v_reinterpret_as_s16(v_src00), v_reinterpret_as_s16(v_src10), v_tmp0, v_tmp1);
2252             v_res0 += v_dotprod(v_tmp0, v_mul23);
2253             v_res1 += v_dotprod(v_tmp1, v_mul23);
2254             v_zip(v_reinterpret_as_s16(v_src01), v_reinterpret_as_s16(v_src11), v_tmp0, v_tmp1);
2255             v_res2 += v_dotprod(v_tmp0, v_mul23);
2256             v_res3 += v_dotprod(v_tmp1, v_mul23);
2257
2258             v_int32x4 v_resj0, v_resj1, v_resj2, v_resj3;
2259             v_expand(v_load(src + 2*cn), v_src00, v_src01);
2260             v_mul_expand(v_reinterpret_as_s16(v_src00), v_mul4, v_resj0, v_resj1);
2261             v_mul_expand(v_reinterpret_as_s16(v_src01), v_mul4, v_resj2, v_resj3);
2262             v_res0 += v_resj0;
2263             v_res1 += v_resj1;
2264             v_res2 += v_resj2;
2265             v_res3 += v_resj3;
2266
2267             v_store((uint16_t*)dst, v_pack(v_reinterpret_as_u32(v_res0), v_reinterpret_as_u32(v_res1)));
2268             v_store((uint16_t*)dst + 8, v_pack(v_reinterpret_as_u32(v_res2), v_reinterpret_as_u32(v_res3)));
2269         }
2270         for (; i < lencn; i++, src++, dst++)
2271             *dst = m[0] * src[-2*cn] + m[1] * src[-cn] + m[2] * src[0] + m[3] * src[cn] + m[4] * src[2*cn];
2272
2273         // Points that fall right from border
2274         for (int k = 0; k < cn; k++)
2275         {
2276             dst[k] = m[0] * src[k - 2 * cn] + m[1] * src[k - cn] + m[2] * src[k] + m[3] * src[k + cn];
2277             dst[k + cn] = m[0] * src[k - cn] + m[1] * src[k] + m[2] * src[k + cn];
2278         }
2279         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2280         {
2281             int idxp1 = (borderInterpolate(len, len, borderType) - (len - 2))*cn;
2282             int idxp2 = (borderInterpolate(len + 1, len, borderType) - (len - 2))*cn;
2283             for (int k = 0; k < cn; k++)
2284             {
2285                 dst[k] = dst[k] + m[4] * src[idxp1 + k];
2286                 dst[k + cn] = dst[k + cn] + m[3] * src[idxp1 + k] + m[4] * src[idxp2 + k];
2287             }
2288         }
2289     }
2290 }
2291 template <typename ET, typename FT>
2292 void hlineSmooth5N14641(const ET* src, int cn, const FT*, int, FT* dst, int len, int borderType)
2293 {
2294     if (len == 1)
2295     {
2296         if (borderType == BORDER_CONSTANT)
2297             for (int k = 0; k < cn; k++)
2298                 dst[k] = (FT(src[k])>>3)*3;
2299         else
2300             for (int k = 0; k < cn; k++)
2301                 dst[k] = src[k];
2302     }
2303     else if (len == 2)
2304     {
2305         if (borderType == BORDER_CONSTANT)
2306             for (int k = 0; k < cn; k++)
2307             {
2308                 dst[k] = (FT(src[k])>>4)*6 + (FT(src[k + cn])>>2);
2309                 dst[k + cn] = (FT(src[k]) >> 2) + (FT(src[k + cn])>>4)*6;
2310             }
2311         else
2312         {
2313             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2314             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2315             int idxp1 = borderInterpolate(2, len, borderType)*cn;
2316             int idxp2 = borderInterpolate(3, len, borderType)*cn;
2317             for (int k = 0; k < cn; k++)
2318             {
2319                 dst[k] = (FT(src[k])>>4)*6 + (FT(src[k + idxm1])>>2) + (FT(src[k + cn])>>2) + (FT(src[k + idxp1])>>4) + (FT(src[k + idxm2])>>4);
2320                 dst[k + cn] = (FT(src[k + cn])>>4)*6 + (FT(src[k])>>2) + (FT(src[k + idxp1])>>2) + (FT(src[k + idxm1])>>4) + (FT(src[k + idxp2])>>4);
2321             }
2322         }
2323     }
2324     else if (len == 3)
2325     {
2326         if (borderType == BORDER_CONSTANT)
2327             for (int k = 0; k < cn; k++)
2328             {
2329                 dst[k] = (FT(src[k])>>4)*6 + (FT(src[k + cn])>>2) + (FT(src[k + 2 * cn])>>4);
2330                 dst[k + cn] = (FT(src[k + cn])>>4)*6 + (FT(src[k])>>2) + (FT(src[k + 2 * cn])>>2);
2331                 dst[k + 2 * cn] = (FT(src[k + 2 * cn])>>4)*6 + (FT(src[k + cn])>>2) + (FT(src[k])>>4);
2332             }
2333         else
2334         {
2335             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2336             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2337             int idxp1 = borderInterpolate(3, len, borderType)*cn;
2338             int idxp2 = borderInterpolate(4, len, borderType)*cn;
2339             for (int k = 0; k < cn; k++)
2340             {
2341                 dst[k] = (FT(src[k])>>4)*6 + (FT(src[k + cn])>>2) + (FT(src[k + idxm1])>>2) + (FT(src[k + 2 * cn])>>4) + (FT(src[k + idxm2])>>4);
2342                 dst[k + cn] = (FT(src[k + cn])>>4)*6 + (FT(src[k])>>2) + (FT(src[k + 2 * cn])>>2) + (FT(src[k + idxm1])>>4) + (FT(src[k + idxp1])>>4);
2343                 dst[k + 2 * cn] = (FT(src[k + 2 * cn])>>4)*6 + (FT(src[k + cn])>>2) + (FT(src[k + idxp1])>>2) + (FT(src[k])>>4) + (FT(src[k + idxp2])>>4);
2344             }
2345         }
2346     }
2347     else
2348     {
2349         // Points that fall left from border
2350         for (int k = 0; k < cn; k++)
2351         {
2352             dst[k] = (FT(src[k])>>4)*6 + (FT(src[cn + k])>>2) + (FT(src[2 * cn + k])>>4);
2353             dst[k + cn] = (FT(src[cn + k])>>4)*6 + (FT(src[k])>>2) + (FT(src[2 * cn + k])>>2) + (FT(src[3 * cn + k])>>4);
2354         }
2355         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2356         {
2357             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2358             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2359             for (int k = 0; k < cn; k++)
2360             {
2361                 dst[k] = dst[k] + (FT(src[idxm2 + k])>>4) + (FT(src[idxm1 + k])>>2);
2362                 dst[k + cn] = dst[k + cn] + (FT(src[idxm1 + k])>>4);
2363             }
2364         }
2365
2366         src += 2 * cn; dst += 2 * cn;
2367         for (int i = 2 * cn; i < (len - 2)*cn; i++, src++, dst++)
2368             *dst = (FT(src[0])>>4)*6 + (FT(src[-cn])>>2) + (FT(src[cn])>>2) + (FT(src[-2 * cn])>>4) + (FT(src[2 * cn])>>4);
2369
2370         // Points that fall right from border
2371         for (int k = 0; k < cn; k++)
2372         {
2373             dst[k] = (FT(src[k])>>4)*6 + (FT(src[k - cn])>>2) + (FT(src[k + cn])>>2) + (FT(src[k - 2 * cn])>>4);
2374             dst[k + cn] = (FT(src[k + cn])>>4)*6 + (FT(src[k])>>2) + (FT(src[k - cn])>>4);
2375         }
2376         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2377         {
2378             int idxp1 = (borderInterpolate(len, len, borderType) - (len - 2))*cn;
2379             int idxp2 = (borderInterpolate(len + 1, len, borderType) - (len - 2))*cn;
2380             for (int k = 0; k < cn; k++)
2381             {
2382                 dst[k] = dst[k] + (FT(src[idxp1 + k])>>4);
2383                 dst[k + cn] = dst[k + cn] + (FT(src[idxp1 + k])>>2) + (FT(src[idxp2 + k])>>4);
2384             }
2385         }
2386     }
2387 }
2388 template <>
2389 void hlineSmooth5N14641<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16*, int, ufixedpoint16* dst, int len, int borderType)
2390 {
2391     if (len == 1)
2392     {
2393         if (borderType == BORDER_CONSTANT)
2394             for (int k = 0; k < cn; k++)
2395                 dst[k] = (ufixedpoint16(src[k])>>3) * 3;
2396         else
2397         {
2398             for (int k = 0; k < cn; k++)
2399                 dst[k] = src[k];
2400         }
2401     }
2402     else if (len == 2)
2403     {
2404         if (borderType == BORDER_CONSTANT)
2405             for (int k = 0; k < cn; k++)
2406             {
2407                 dst[k] = (ufixedpoint16(src[k]) >> 4) * 6 + (ufixedpoint16(src[k + cn]) >> 2);
2408                 dst[k + cn] = (ufixedpoint16(src[k]) >> 2) + (ufixedpoint16(src[k + cn]) >> 4) * 6;
2409             }
2410         else
2411         {
2412             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2413             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2414             int idxp1 = borderInterpolate(2, len, borderType)*cn;
2415             int idxp2 = borderInterpolate(3, len, borderType)*cn;
2416             for (int k = 0; k < cn; k++)
2417             {
2418                 dst[k] = (ufixedpoint16(src[k]) >> 4) * 6 + (ufixedpoint16(src[k + idxm1]) >> 2) + (ufixedpoint16(src[k + cn]) >> 2) + (ufixedpoint16(src[k + idxp1]) >> 4) + (ufixedpoint16(src[k + idxm2]) >> 4);
2419                 dst[k + cn] = (ufixedpoint16(src[k + cn]) >> 4) * 6 + (ufixedpoint16(src[k]) >> 2) + (ufixedpoint16(src[k + idxp1]) >> 2) + (ufixedpoint16(src[k + idxm1]) >> 4) + (ufixedpoint16(src[k + idxp2]) >> 4);
2420             }
2421         }
2422     }
2423     else if (len == 3)
2424     {
2425         if (borderType == BORDER_CONSTANT)
2426             for (int k = 0; k < cn; k++)
2427             {
2428                 dst[k] = (ufixedpoint16(src[k]) >> 4) * 6 + (ufixedpoint16(src[k + cn]) >> 2) + (ufixedpoint16(src[k + 2 * cn]) >> 4);
2429                 dst[k + cn] = (ufixedpoint16(src[k + cn]) >> 4) * 6 + (ufixedpoint16(src[k]) >> 2) + (ufixedpoint16(src[k + 2 * cn]) >> 2);
2430                 dst[k + 2 * cn] = (ufixedpoint16(src[k + 2 * cn]) >> 4) * 6 + (ufixedpoint16(src[k + cn]) >> 2) + (ufixedpoint16(src[k]) >> 4);
2431             }
2432         else
2433         {
2434             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2435             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2436             int idxp1 = borderInterpolate(3, len, borderType)*cn;
2437             int idxp2 = borderInterpolate(4, len, borderType)*cn;
2438             for (int k = 0; k < cn; k++)
2439             {
2440                 dst[k] = (ufixedpoint16(src[k]) >> 4) * 6 + (ufixedpoint16(src[k + cn]) >> 2) + (ufixedpoint16(src[k + idxm1]) >> 2) + (ufixedpoint16(src[k + 2 * cn]) >> 4) + (ufixedpoint16(src[k + idxm2]) >> 4);
2441                 dst[k + cn] = (ufixedpoint16(src[k + cn]) >> 4) * 6 + (ufixedpoint16(src[k]) >> 2) + (ufixedpoint16(src[k + 2 * cn]) >> 2) + (ufixedpoint16(src[k + idxm1]) >> 4) + (ufixedpoint16(src[k + idxp1]) >> 4);
2442                 dst[k + 2 * cn] = (ufixedpoint16(src[k + 2 * cn]) >> 4) * 6 + (ufixedpoint16(src[k + cn]) >> 2) + (ufixedpoint16(src[k + idxp1]) >> 2) + (ufixedpoint16(src[k]) >> 4) + (ufixedpoint16(src[k + idxp2]) >> 4);
2443             }
2444         }
2445     }
2446     else
2447     {
2448         // Points that fall left from border
2449         for (int k = 0; k < cn; k++)
2450         {
2451             dst[k] = (ufixedpoint16(src[k]) >> 4) * 6 + (ufixedpoint16(src[cn + k]) >> 2) + (ufixedpoint16(src[2 * cn + k]) >> 4);
2452             dst[k + cn] = (ufixedpoint16(src[cn + k]) >> 4) * 6 + (ufixedpoint16(src[k]) >> 2) + (ufixedpoint16(src[2 * cn + k]) >> 2) + (ufixedpoint16(src[3 * cn + k]) >> 4);
2453         }
2454         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2455         {
2456             int idxm2 = borderInterpolate(-2, len, borderType)*cn;
2457             int idxm1 = borderInterpolate(-1, len, borderType)*cn;
2458             for (int k = 0; k < cn; k++)
2459             {
2460                 dst[k] = dst[k] + (ufixedpoint16(src[idxm2 + k]) >> 4) + (ufixedpoint16(src[idxm1 + k]) >> 2);
2461                 dst[k + cn] = dst[k + cn] + (ufixedpoint16(src[idxm1 + k]) >> 4);
2462             }
2463         }
2464
2465         src += 2 * cn; dst += 2 * cn;
2466         int i = 2 * cn, lencn = (len - 2)*cn;
2467         v_uint16x8 v_6 = v_setall_u16(6);
2468         for (; i < lencn - 15; i += 16, src += 16, dst += 16)
2469         {
2470             v_uint16x8 v_src00, v_src01, v_src10, v_src11, v_src20, v_src21, v_src30, v_src31, v_src40, v_src41;
2471             v_expand(v_load(src - 2*cn), v_src00, v_src01);
2472             v_expand(v_load(src - cn), v_src10, v_src11);
2473             v_expand(v_load(src), v_src20, v_src21);
2474             v_expand(v_load(src + cn), v_src30, v_src31);
2475             v_expand(v_load(src + 2*cn), v_src40, v_src41);
2476             v_store((uint16_t*)dst, (v_src20 * v_6 + ((v_src10 + v_src30) << 2) + v_src00 + v_src40) << 4);
2477             v_store((uint16_t*)dst + 8, (v_src21 * v_6 + ((v_src11 + v_src31) << 2) + v_src01 + v_src41) << 4);
2478         }
2479         for (; i < lencn; i++, src++, dst++)
2480             *((uint16_t*)dst) = (uint16_t(src[0]) * 6 + ((uint16_t(src[-cn]) + uint16_t(src[cn])) << 2) + uint16_t(src[-2 * cn]) + uint16_t(src[2 * cn])) << 4;
2481
2482         // Points that fall right from border
2483         for (int k = 0; k < cn; k++)
2484         {
2485             dst[k] = (ufixedpoint16(src[k]) >> 4) * 6 + (ufixedpoint16(src[k - cn]) >> 2) + (ufixedpoint16(src[k + cn]) >> 2) + (ufixedpoint16(src[k - 2 * cn]) >> 4);
2486             dst[k + cn] = (ufixedpoint16(src[k + cn]) >> 4) * 6 + (ufixedpoint16(src[k]) >> 2) + (ufixedpoint16(src[k - cn]) >> 4);
2487         }
2488         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2489         {
2490             int idxp1 = (borderInterpolate(len, len, borderType) - (len - 2))*cn;
2491             int idxp2 = (borderInterpolate(len + 1, len, borderType) - (len - 2))*cn;
2492             for (int k = 0; k < cn; k++)
2493             {
2494                 dst[k] = dst[k] + (ufixedpoint16(src[idxp1 + k]) >> 4);
2495                 dst[k + cn] = dst[k + cn] + (ufixedpoint16(src[idxp1 + k]) >> 2) + (ufixedpoint16(src[idxp2 + k]) >> 4);
2496             }
2497         }
2498     }
2499 }
2500 template <typename ET, typename FT>
2501 void hlineSmooth(const ET* src, int cn, const FT* m, int n, FT* dst, int len, int borderType)
2502 {
2503     int pre_shift = n / 2;
2504     int post_shift = n - pre_shift;
2505     int i = 0;
2506     for (; i < min(pre_shift, len); i++, dst += cn) // Points that fall left from border
2507     {
2508         for (int k = 0; k < cn; k++)
2509             dst[k] = m[pre_shift-i] * src[k];
2510         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2511             for (int j = i - pre_shift, mid = 0; j < 0; j++, mid++)
2512             {
2513                 int src_idx = borderInterpolate(j, len, borderType);
2514                 for (int k = 0; k < cn; k++)
2515                     dst[k] = dst[k] + m[mid] * src[src_idx*cn + k];
2516             }
2517         int j, mid;
2518         for (j = 1, mid = pre_shift - i + 1; j < min(i + post_shift, len); j++, mid++)
2519             for (int k = 0; k < cn; k++)
2520                 dst[k] = dst[k] + m[mid] * src[j*cn + k];
2521         if (borderType != BORDER_CONSTANT)
2522             for (; j < i + post_shift; j++, mid++)
2523             {
2524                 int src_idx = borderInterpolate(j, len, borderType);
2525                 for (int k = 0; k < cn; k++)
2526                     dst[k] = dst[k] + m[mid] * src[src_idx*cn + k];
2527             }
2528     }
2529     i *= cn;
2530     for (; i < (len - post_shift + 1)*cn; i++, src++, dst++)
2531     {
2532         *dst = m[0] * src[0];
2533         for (int j = 1; j < n; j++)
2534             *dst = *dst + m[j] * src[j*cn];
2535     }
2536     i /= cn;
2537     for (i -= pre_shift; i < len - pre_shift; i++, src += cn, dst += cn) // Points that fall right from border
2538     {
2539         for (int k = 0; k < cn; k++)
2540             dst[k] = m[0] * src[k];
2541         int j = 1;
2542         for (; j < len - i; j++)
2543             for (int k = 0; k < cn; k++)
2544                 dst[k] = dst[k] + m[j] * src[j*cn + k];
2545         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2546             for (; j < n; j++)
2547             {
2548                 int src_idx = borderInterpolate(i + j, len, borderType) - i;
2549                 for (int k = 0; k < cn; k++)
2550                     dst[k] = dst[k] + m[j] * src[src_idx*cn + k];
2551             }
2552     }
2553 }
2554 template <>
2555 void hlineSmooth<uint8_t, ufixedpoint16>(const uint8_t* src, int cn, const ufixedpoint16* m, int n, ufixedpoint16* dst, int len, int borderType)
2556 {
2557     int pre_shift = n / 2;
2558     int post_shift = n - pre_shift;
2559     int i = 0;
2560     for (; i < min(pre_shift, len); i++, dst += cn) // Points that fall left from border
2561     {
2562         for (int k = 0; k < cn; k++)
2563             dst[k] = m[pre_shift - i] * src[k];
2564         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2565             for (int j = i - pre_shift, mid = 0; j < 0; j++, mid++)
2566             {
2567                 int src_idx = borderInterpolate(j, len, borderType);
2568                 for (int k = 0; k < cn; k++)
2569                     dst[k] = dst[k] + m[mid] * src[src_idx*cn + k];
2570             }
2571         int j, mid;
2572         for (j = 1, mid = pre_shift - i + 1; j < min(i + post_shift, len); j++, mid++)
2573             for (int k = 0; k < cn; k++)
2574                 dst[k] = dst[k] + m[mid] * src[j*cn + k];
2575         if (borderType != BORDER_CONSTANT)
2576             for (; j < i + post_shift; j++, mid++)
2577             {
2578                 int src_idx = borderInterpolate(j, len, borderType);
2579                 for (int k = 0; k < cn; k++)
2580                     dst[k] = dst[k] + m[mid] * src[src_idx*cn + k];
2581             }
2582     }
2583     i *= cn;
2584     int lencn = (len - post_shift + 1)*cn;
2585     for (; i < lencn - 15; i+=16, src+=16, dst+=16)
2586     {
2587         v_uint16x8 v_src00, v_src01, v_src10, v_src11;
2588         v_int16x8 v_tmp0, v_tmp1;
2589
2590         v_int16x8 v_mul = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)m)));
2591
2592         v_expand(v_load(src), v_src00, v_src01);
2593         v_expand(v_load(src+cn), v_src10, v_src11);
2594         v_zip(v_reinterpret_as_s16(v_src00), v_reinterpret_as_s16(v_src10), v_tmp0, v_tmp1);
2595         v_int32x4 v_res0 = v_dotprod(v_tmp0, v_mul);
2596         v_int32x4 v_res1 = v_dotprod(v_tmp1, v_mul);
2597         v_zip(v_reinterpret_as_s16(v_src01), v_reinterpret_as_s16(v_src11), v_tmp0, v_tmp1);
2598         v_int32x4 v_res2 = v_dotprod(v_tmp0, v_mul);
2599         v_int32x4 v_res3 = v_dotprod(v_tmp1, v_mul);
2600
2601         int j = 2;
2602         for (; j < n - 1; j += 2)
2603         {
2604             v_mul = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)(m + j))));
2605
2606             v_expand(v_load(src + j * cn), v_src00, v_src01);
2607             v_expand(v_load(src + (j + 1) * cn), v_src10, v_src11);
2608             v_zip(v_reinterpret_as_s16(v_src00), v_reinterpret_as_s16(v_src10), v_tmp0, v_tmp1);
2609             v_res0 += v_dotprod(v_tmp0, v_mul);
2610             v_res1 += v_dotprod(v_tmp1, v_mul);
2611             v_zip(v_reinterpret_as_s16(v_src01), v_reinterpret_as_s16(v_src11), v_tmp0, v_tmp1);
2612             v_res2 += v_dotprod(v_tmp0, v_mul);
2613             v_res3 += v_dotprod(v_tmp1, v_mul);
2614         }
2615         if (j < n)
2616         {
2617             v_int32x4 v_resj0, v_resj1, v_resj2, v_resj3;
2618             v_mul = v_reinterpret_as_s16(v_setall_u16(*((uint16_t*)(m + j))));
2619             v_expand(v_load(src + j * cn), v_src00, v_src01);
2620             v_mul_expand(v_reinterpret_as_s16(v_src00), v_mul, v_resj0, v_resj1);
2621             v_mul_expand(v_reinterpret_as_s16(v_src01), v_mul, v_resj2, v_resj3);
2622             v_res0 += v_resj0;
2623             v_res1 += v_resj1;
2624             v_res2 += v_resj2;
2625             v_res3 += v_resj3;
2626         }
2627
2628         v_store((uint16_t*)dst, v_pack(v_reinterpret_as_u32(v_res0), v_reinterpret_as_u32(v_res1)));
2629         v_store((uint16_t*)dst+8, v_pack(v_reinterpret_as_u32(v_res2), v_reinterpret_as_u32(v_res3)));
2630     }
2631     for (; i < lencn; i++, src++, dst++)
2632     {
2633             *dst = m[0] * src[0];
2634             for (int j = 1; j < n; j++)
2635                 *dst = *dst + m[j] * src[j*cn];
2636     }
2637     i /= cn;
2638     for (i -= pre_shift; i < len - pre_shift; i++, src += cn, dst += cn) // Points that fall right from border
2639     {
2640         for (int k = 0; k < cn; k++)
2641             dst[k] = m[0] * src[k];
2642         int j = 1;
2643         for (; j < len - i; j++)
2644             for (int k = 0; k < cn; k++)
2645                 dst[k] = dst[k] + m[j] * src[j*cn + k];
2646         if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2647             for (; j < n; j++)
2648             {
2649                 int src_idx = borderInterpolate(i + j, len, borderType) - i;
2650                 for (int k = 0; k < cn; k++)
2651                     dst[k] = dst[k] + m[j] * src[src_idx*cn + k];
2652             }
2653     }
2654 }
2655 template <typename ET, typename FT>
2656 void vlineSmooth1N(const FT* const * src, const FT* m, int, ET* dst, int len)
2657 {
2658     const FT* src0 = src[0];
2659     for (int i = 0; i < len; i++)
2660         dst[i] = m * src0[i];
2661 }
2662 template <>
2663 void vlineSmooth1N<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16* m, int, uint8_t* dst, int len)
2664 {
2665     const ufixedpoint16* src0 = src[0];
2666     v_uint16x8 v_mul = v_setall_u16(*((uint16_t*)m));
2667     int i = 0;
2668     for (; i < len - 7; i += 8)
2669     {
2670         v_uint16x8 v_src0 = v_load((uint16_t*)src0 + i);
2671         v_uint32x4 v_res0, v_res1;
2672         v_mul_expand(v_src0, v_mul, v_res0, v_res1);
2673         v_pack_store(dst + i, v_rshr_pack<16>(v_res0, v_res1));
2674     }
2675     for (; i < len; i++)
2676         dst[i] = m[0] * src0[i];
2677 }
2678 template <typename ET, typename FT>
2679 void vlineSmooth1N1(const FT* const * src, const FT*, int, ET* dst, int len)
2680 {
2681     const FT* src0 = src[0];
2682     for (int i = 0; i < len; i++)
2683         dst[i] = src0[i];
2684 }
2685 template <>
2686 void vlineSmooth1N1<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16*, int, uint8_t* dst, int len)
2687 {
2688     const ufixedpoint16* src0 = src[0];
2689     int i = 0;
2690     for (; i < len - 7; i += 8)
2691         v_rshr_pack_store<8>(dst + i, v_load((uint16_t*)(src0 + i)));
2692     for (; i < len; i++)
2693         dst[i] = src0[i];
2694 }
2695 template <typename ET, typename FT>
2696 void vlineSmooth3N(const FT* const * src, const FT* m, int, ET* dst, int len)
2697 {
2698     for (int i = 0; i < len; i++)
2699         dst[i] = m[0] * src[0][i] + m[1] * src[1][i] + m[2] * src[2][i];
2700 }
2701 template <>
2702 void vlineSmooth3N<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16* m, int, uint8_t* dst, int len)
2703 {
2704     static const v_int16x8 v_128 = v_reinterpret_as_s16(v_setall_u16((uint16_t)1 << 15));
2705
2706     v_int32x4 v_128_4 = v_setall_s32(128 << 16);
2707     if (len > 7)
2708     {
2709         ufixedpoint32 val[] = { (m[0] + m[1] + m[2]) * ufixedpoint16((uint8_t)128) };
2710         v_128_4 = v_setall_s32(*((int32_t*)val));
2711     }
2712
2713     int i = 0;
2714     v_int16x8 v_mul01 = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)m)));
2715     v_int16x8 v_mul2 = v_reinterpret_as_s16(v_setall_u16(*((uint16_t*)(m + 2))));
2716     for (; i < len - 7; i += 8)
2717     {
2718         v_int16x8 v_src0, v_src1;
2719         v_int16x8 v_tmp0, v_tmp1;
2720
2721         v_src0 = v_load((int16_t*)(src[0]) + i);
2722         v_src1 = v_load((int16_t*)(src[1]) + i);
2723         v_zip(v_add_wrap(v_src0, v_128), v_add_wrap(v_src1, v_128), v_tmp0, v_tmp1);
2724         v_int32x4 v_res0 = v_dotprod(v_tmp0, v_mul01);
2725         v_int32x4 v_res1 = v_dotprod(v_tmp1, v_mul01);
2726
2727         v_int32x4 v_resj0, v_resj1;
2728         v_src0 = v_load((int16_t*)(src[2]) + i);
2729         v_mul_expand(v_add_wrap(v_src0, v_128), v_mul2, v_resj0, v_resj1);
2730         v_res0 += v_resj0;
2731         v_res1 += v_resj1;
2732
2733         v_res0 += v_128_4;
2734         v_res1 += v_128_4;
2735
2736         v_uint16x8 v_res = v_reinterpret_as_u16(v_rshr_pack<16>(v_res0, v_res1));
2737         v_pack_store(dst + i, v_res);
2738     }
2739     for (; i < len; i++)
2740         dst[i] = m[0] * src[0][i] + m[1] * src[1][i] + m[2] * src[2][i];
2741 }
2742 template <typename ET, typename FT>
2743 void vlineSmooth3N121(const FT* const * src, const FT*, int, ET* dst, int len)
2744 {
2745     for (int i = 0; i < len; i++)
2746         dst[i] = ((FT::WT(src[0][i]) + FT::WT(src[2][i])) >> 2) + (FT::WT(src[1][i]) >> 1);
2747 }
2748 template <>
2749 void vlineSmooth3N121<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16*, int, uint8_t* dst, int len)
2750 {
2751     int i = 0;
2752     for (; i < len - 7; i += 8)
2753     {
2754         v_uint32x4 v_src00, v_src01, v_src10, v_src11, v_src20, v_src21;
2755         v_expand(v_load((uint16_t*)(src[0]) + i), v_src00, v_src01);
2756         v_expand(v_load((uint16_t*)(src[1]) + i), v_src10, v_src11);
2757         v_expand(v_load((uint16_t*)(src[2]) + i), v_src20, v_src21);
2758         v_uint16x8 v_res = v_rshr_pack<10>(v_src00 + v_src20 + (v_src10 << 1), v_src01 + v_src21 + (v_src11 << 1));
2759         v_pack_store(dst + i, v_res);
2760     }
2761     for (; i < len; i++)
2762         dst[i] = (((uint32_t)(((uint16_t*)(src[0]))[i]) + (uint32_t)(((uint16_t*)(src[2]))[i]) + ((uint32_t)(((uint16_t*)(src[1]))[i]) << 1)) + (1 << 9)) >> 10;
2763 }
2764 template <typename ET, typename FT>
2765 void vlineSmooth5N(const FT* const * src, const FT* m, int, ET* dst, int len)
2766 {
2767     for (int i = 0; i < len; i++)
2768         dst[i] = m[0] * src[0][i] + m[1] * src[1][i] + m[2] * src[2][i] + m[3] * src[3][i] + m[4] * src[4][i];
2769 }
2770 template <>
2771 void vlineSmooth5N<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16* m, int, uint8_t* dst, int len)
2772 {
2773     static const v_int16x8 v_128 = v_reinterpret_as_s16(v_setall_u16((uint16_t)1 << 15));
2774
2775     v_int32x4 v_128_4 = v_setall_s32(128 << 16);
2776     if (len > 7)
2777     {
2778         ufixedpoint32 val[] = { (m[0] + m[1] + m[2] + m[3] + m[4]) * ufixedpoint16((uint8_t)128) };
2779         v_128_4 = v_setall_s32(*((int32_t*)val));
2780     }
2781
2782     int i = 0;
2783     v_int16x8 v_mul01 = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)m)));
2784     v_int16x8 v_mul23 = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)(m + 2))));
2785     v_int16x8 v_mul4 = v_reinterpret_as_s16(v_setall_u16(*((uint16_t*)(m + 4))));
2786     for (; i < len - 7; i += 8)
2787     {
2788         v_int16x8 v_src0, v_src1;
2789         v_int16x8 v_tmp0, v_tmp1;
2790
2791         v_src0 = v_load((int16_t*)(src[0]) + i);
2792         v_src1 = v_load((int16_t*)(src[1]) + i);
2793         v_zip(v_add_wrap(v_src0, v_128), v_add_wrap(v_src1, v_128), v_tmp0, v_tmp1);
2794         v_int32x4 v_res0 = v_dotprod(v_tmp0, v_mul01);
2795         v_int32x4 v_res1 = v_dotprod(v_tmp1, v_mul01);
2796
2797         v_src0 = v_load((int16_t*)(src[2]) + i);
2798         v_src1 = v_load((int16_t*)(src[3]) + i);
2799         v_zip(v_add_wrap(v_src0, v_128), v_add_wrap(v_src1, v_128), v_tmp0, v_tmp1);
2800         v_res0 += v_dotprod(v_tmp0, v_mul23);
2801         v_res1 += v_dotprod(v_tmp1, v_mul23);
2802
2803         v_int32x4 v_resj0, v_resj1;
2804         v_src0 = v_load((int16_t*)(src[4]) + i);
2805         v_mul_expand(v_add_wrap(v_src0, v_128), v_mul4, v_resj0, v_resj1);
2806         v_res0 += v_resj0;
2807         v_res1 += v_resj1;
2808
2809         v_res0 += v_128_4;
2810         v_res1 += v_128_4;
2811
2812         v_uint16x8 v_res = v_reinterpret_as_u16(v_rshr_pack<16>(v_res0, v_res1));
2813         v_pack_store(dst + i, v_res);
2814     }
2815     for (; i < len; i++)
2816         dst[i] = m[0] * src[0][i] + m[1] * src[1][i] + m[2] * src[2][i] + m[3] * src[3][i] + m[4] * src[4][i];
2817 }
2818 template <typename ET, typename FT>
2819 void vlineSmooth5N14641(const FT* const * src, const FT*, int, ET* dst, int len)
2820 {
2821     for (int i = 0; i < len; i++)
2822         dst[i] = (FT::WT(src[2][i])*6 + ((FT::WT(src[1][i]) + FT::WT(src[3][i]))<<2) + FT::WT(src[0][i]) + FT::WT(src[4][i])) >> 4;
2823 }
2824 template <>
2825 void vlineSmooth5N14641<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16*, int, uint8_t* dst, int len)
2826 {
2827     int i = 0;
2828     v_uint32x4 v_6 = v_setall_u32(6);
2829     for (; i < len - 7; i += 8)
2830     {
2831         v_uint32x4 v_src00, v_src01, v_src10, v_src11, v_src20, v_src21, v_src30, v_src31, v_src40, v_src41;
2832         v_expand(v_load((uint16_t*)(src[0]) + i), v_src00, v_src01);
2833         v_expand(v_load((uint16_t*)(src[1]) + i), v_src10, v_src11);
2834         v_expand(v_load((uint16_t*)(src[2]) + i), v_src20, v_src21);
2835         v_expand(v_load((uint16_t*)(src[3]) + i), v_src30, v_src31);
2836         v_expand(v_load((uint16_t*)(src[4]) + i), v_src40, v_src41);
2837         v_uint16x8 v_res = v_rshr_pack<12>(v_src20*v_6 + ((v_src10 + v_src30) << 2) + v_src00 + v_src40,
2838                                            v_src21*v_6 + ((v_src11 + v_src31) << 2) + v_src01 + v_src41);
2839         v_pack_store(dst + i, v_res);
2840     }
2841     for (; i < len; i++)
2842         dst[i] = ((uint32_t)(((uint16_t*)(src[2]))[i]) * 6 +
2843                   (((uint32_t)(((uint16_t*)(src[1]))[i]) + (uint32_t)(((uint16_t*)(src[3]))[i])) << 2) +
2844                   (uint32_t)(((uint16_t*)(src[0]))[i]) + (uint32_t)(((uint16_t*)(src[4]))[i]) + (1 << 11)) >> 12;
2845 }
2846 template <typename ET, typename FT>
2847 void vlineSmooth(const FT* const * src, const FT* m, int n, ET* dst, int len)
2848 {
2849     for (int i = 0; i < len; i++)
2850     {
2851         typename FT::WT val = m[0] * src[0][i];
2852         for (int j = 1; j < n; j++)
2853             val = val + m[j] * src[j][i];
2854         dst[i] = val;
2855     }
2856 }
2857 template <>
2858 void vlineSmooth<uint8_t, ufixedpoint16>(const ufixedpoint16* const * src, const ufixedpoint16* m, int n, uint8_t* dst, int len)
2859 {
2860     static const v_int16x8 v_128 = v_reinterpret_as_s16(v_setall_u16((uint16_t)1 << 15));
2861
2862     v_int32x4 v_128_4 = v_setall_s32(128 << 16);
2863     if (len > 7)
2864     {
2865         ufixedpoint16 msum = m[0] + m[1];
2866         for (int j = 2; j < n; j++)
2867             msum = msum + m[j];
2868         ufixedpoint32 val[] = { msum * ufixedpoint16((uint8_t)128) };
2869         v_128_4 = v_setall_s32(*((int32_t*)val));
2870     }
2871
2872     int i = 0;
2873     for (; i < len - 7; i += 8)
2874     {
2875         v_int16x8 v_src0, v_src1;
2876         v_int16x8 v_tmp0, v_tmp1;
2877
2878         v_int16x8 v_mul = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)m)));
2879
2880         v_src0 = v_load((int16_t*)(src[0]) + i);
2881         v_src1 = v_load((int16_t*)(src[1]) + i);
2882         v_zip(v_add_wrap(v_src0, v_128), v_add_wrap(v_src1, v_128), v_tmp0, v_tmp1);
2883         v_int32x4 v_res0 = v_dotprod(v_tmp0, v_mul);
2884         v_int32x4 v_res1 = v_dotprod(v_tmp1, v_mul);
2885
2886         int j = 2;
2887         for (; j < n - 1; j+=2)
2888         {
2889             v_mul = v_reinterpret_as_s16(v_setall_u32(*((uint32_t*)(m+j))));
2890
2891             v_src0 = v_load((int16_t*)(src[j]) + i);
2892             v_src1 = v_load((int16_t*)(src[j+1]) + i);
2893             v_zip(v_add_wrap(v_src0, v_128), v_add_wrap(v_src1, v_128), v_tmp0, v_tmp1);
2894             v_res0 += v_dotprod(v_tmp0, v_mul);
2895             v_res1 += v_dotprod(v_tmp1, v_mul);
2896         }
2897         if(j < n)
2898         {
2899             v_int32x4 v_resj0, v_resj1;
2900             v_mul = v_reinterpret_as_s16(v_setall_u16(*((uint16_t*)(m + j))));
2901             v_src0 = v_load((int16_t*)(src[j]) + i);
2902             v_mul_expand(v_add_wrap(v_src0, v_128), v_mul, v_resj0, v_resj1);
2903             v_res0 += v_resj0;
2904             v_res1 += v_resj1;
2905         }
2906         v_res0 += v_128_4;
2907         v_res1 += v_128_4;
2908
2909         v_uint16x8 v_res = v_reinterpret_as_u16(v_rshr_pack<16>(v_res0, v_res1));
2910         v_pack_store(dst + i, v_res);
2911     }
2912     for (; i < len; i++)
2913     {
2914         ufixedpoint32 val = m[0] * src[0][i];
2915         for (int j = 1; j < n; j++)
2916         {
2917             val = val + m[j] * src[j][i];
2918         }
2919         dst[i] = val;
2920     }
2921 }
2922 template <typename ET, typename FT>
2923 class fixedSmoothInvoker : public ParallelLoopBody
2924 {
2925 public:
2926     fixedSmoothInvoker(const ET* _src, size_t _src_stride, ET* _dst, size_t _dst_stride,
2927                        int _width, int _height, int _cn, const FT* _kx, int _kxlen, const FT* _ky, int _kylen, int _borderType) : ParallelLoopBody(),
2928                        src(_src), dst(_dst), src_stride(_src_stride), dst_stride(_dst_stride),
2929                        width(_width), height(_height), cn(_cn), kx(_kx), ky(_ky), kxlen(_kxlen), kylen(_kylen), borderType(_borderType)
2930     {
2931         if (kxlen == 1)
2932         {
2933             if ((kx[0] - FT::one()).isZero())
2934                 hlineSmoothFunc = hlineSmooth1N1;
2935             else
2936                 hlineSmoothFunc = hlineSmooth1N;
2937         }
2938         else if (kxlen == 3)
2939         {
2940             if ((kx[0] - (FT::one()>>2)).isZero()&&(kx[1] - (FT::one()>>1)).isZero()&&(kx[2] - (FT::one()>>2)).isZero())
2941                 hlineSmoothFunc = hlineSmooth3N121;
2942             else
2943                 hlineSmoothFunc = hlineSmooth3N;
2944         }
2945         else if (kxlen == 5)
2946         {
2947             if ((kx[2] - (FT::one()*3>>3)).isZero()&&
2948                 (kx[1] - (FT::one()>>2)).isZero()&&(kx[3] - (FT::one()>>2)).isZero()&&
2949                 (kx[0] - (FT::one()>>4)).isZero()&&(kx[4] - (FT::one()>>4)).isZero())
2950                 hlineSmoothFunc = hlineSmooth5N14641;
2951             else
2952                 hlineSmoothFunc = hlineSmooth5N;
2953         }
2954         else
2955             hlineSmoothFunc = hlineSmooth;
2956         if (kylen == 1)
2957         {
2958             if ((ky[0] - FT::one()).isZero())
2959                 vlineSmoothFunc = vlineSmooth1N1;
2960             else
2961                 vlineSmoothFunc = vlineSmooth1N;
2962         }
2963         else if (kylen == 3)
2964         {
2965             if ((ky[0] - (FT::one() >> 2)).isZero() && (ky[1] - (FT::one() >> 1)).isZero() && (ky[2] - (FT::one() >> 2)).isZero())
2966                 vlineSmoothFunc = vlineSmooth3N121;
2967             else
2968                 vlineSmoothFunc = vlineSmooth3N;
2969         }
2970         else if (kylen == 5)
2971         {
2972             if ((ky[2] - (FT::one() * 3 >> 3)).isZero() &&
2973                 (ky[1] - (FT::one() >> 2)).isZero() && (ky[3] - (FT::one() >> 2)).isZero() &&
2974                 (ky[0] - (FT::one() >> 4)).isZero() && (ky[4] - (FT::one() >> 4)).isZero())
2975                 vlineSmoothFunc = vlineSmooth5N14641;
2976             else
2977                 vlineSmoothFunc = vlineSmooth5N;
2978         }
2979         else
2980             vlineSmoothFunc = vlineSmooth;
2981     }
2982     virtual void operator() (const Range& range) const
2983     {
2984         AutoBuffer<FT> _buf(width*cn*kylen);
2985         FT* buf = _buf;
2986         AutoBuffer<FT*> _ptrs(kylen*2);
2987         FT** ptrs = _ptrs;
2988
2989         if (kylen == 1)
2990         {
2991             ptrs[0] = buf;
2992             for (int i = range.start; i < range.end; i++)
2993             {
2994                 hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[0], width, borderType);
2995                 vlineSmoothFunc(ptrs, ky, kylen, dst + i * dst_stride, width*cn);
2996             }
2997         }
2998         else if (borderType != BORDER_CONSTANT)// If BORDER_CONSTANT out of border values are equal to zero and could be skipped
2999         {
3000             int pre_shift = kylen / 2;
3001             int post_shift = kylen - pre_shift - 1;
3002             // First line evaluation
3003             int idst = range.start;
3004             int ifrom = max(0, idst - pre_shift);
3005             int ito = idst + post_shift + 1;
3006             int i = ifrom;
3007             int bufline = 0;
3008             for (; i < min(ito, height); i++, bufline++)
3009             {
3010                 ptrs[bufline+kylen] = ptrs[bufline] = buf + bufline * width*cn;
3011                 hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3012             }
3013             for (; i < ito; i++, bufline++)
3014             {
3015                 int src_idx = borderInterpolate(i, height, borderType);
3016                 if (src_idx < ifrom)
3017                 {
3018                     ptrs[bufline + kylen] = ptrs[bufline] = buf + bufline * width*cn;
3019                     hlineSmoothFunc(src + src_idx * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3020                 }
3021                 else
3022                 {
3023                     ptrs[bufline + kylen] = ptrs[bufline] = ptrs[src_idx - ifrom];
3024                 }
3025             }
3026             for (int j = idst - pre_shift; j < 0; j++)
3027             {
3028                 int src_idx = borderInterpolate(j, height, borderType);
3029                 if (src_idx >= ito)
3030                 {
3031                     ptrs[2*kylen + j] = ptrs[kylen + j] = buf + (kylen + j) * width*cn;
3032                     hlineSmoothFunc(src + src_idx * src_stride, cn, kx, kxlen, ptrs[kylen + j], width, borderType);
3033                 }
3034                 else
3035                 {
3036                     ptrs[2*kylen + j] = ptrs[kylen + j] = ptrs[src_idx];
3037                 }
3038             }
3039             vlineSmoothFunc(ptrs + bufline, ky, kylen, dst + idst*dst_stride, width*cn); idst++;
3040
3041             // border mode dependent part evaluation
3042             // i points to last src row to evaluate in convolution
3043             bufline %= kylen; ito = min(height, range.end + post_shift);
3044             for (; i < min(kylen, ito); i++, idst++)
3045             {
3046                 ptrs[bufline + kylen] = ptrs[bufline] = buf + bufline * width*cn;
3047                 hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3048                 bufline = (bufline + 1) % kylen;
3049                 vlineSmoothFunc(ptrs + bufline, ky, kylen, dst + idst*dst_stride, width*cn);
3050             }
3051             // Points inside the border
3052             for (; i < ito; i++, idst++)
3053             {
3054                 hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3055                 bufline = (bufline + 1) % kylen;
3056                 vlineSmoothFunc(ptrs + bufline, ky, kylen, dst + idst*dst_stride, width*cn);
3057             }
3058             // Points that could fall below border
3059             for (; i < range.end + post_shift; i++, idst++)
3060             {
3061                 int src_idx = borderInterpolate(i, height, borderType);
3062                 if ((i - src_idx) > kylen)
3063                     hlineSmoothFunc(src + src_idx * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3064                 else
3065                     ptrs[bufline + kylen] = ptrs[bufline] = ptrs[(bufline + kylen - (i - src_idx)) % kylen];
3066                 bufline = (bufline + 1) % kylen;
3067                 vlineSmoothFunc(ptrs + bufline, ky, kylen, dst + idst*dst_stride, width*cn);
3068             }
3069         }
3070         else
3071         {
3072             int pre_shift = kylen / 2;
3073             int post_shift = kylen - pre_shift - 1;
3074             // First line evaluation
3075             int idst = range.start;
3076             int ifrom = idst - pre_shift;
3077             int ito = min(idst + post_shift + 1, height);
3078             int i = max(0, ifrom);
3079             int bufline = 0;
3080             for (; i < ito; i++, bufline++)
3081             {
3082                 ptrs[bufline + kylen] = ptrs[bufline] = buf + bufline * width*cn;
3083                 hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3084             }
3085
3086             if (bufline == 1)
3087                 vlineSmooth1N(ptrs, ky - min(ifrom, 0), bufline, dst + idst*dst_stride, width*cn);
3088             else if (bufline == 3)
3089                 vlineSmooth3N(ptrs, ky - min(ifrom, 0), bufline, dst + idst*dst_stride, width*cn);
3090             else if (bufline == 5)
3091                 vlineSmooth5N(ptrs, ky - min(ifrom, 0), bufline, dst + idst*dst_stride, width*cn);
3092             else
3093                 vlineSmooth(ptrs, ky - min(ifrom, 0), bufline, dst + idst*dst_stride, width*cn);
3094             idst++;
3095
3096             // border mode dependent part evaluation
3097             // i points to last src row to evaluate in convolution
3098             bufline %= kylen; ito = min(height, range.end + post_shift);
3099             for (; i < min(kylen, ito); i++, idst++)
3100             {
3101                 ptrs[bufline + kylen] = ptrs[bufline] = buf + bufline * width*cn;
3102                 hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3103                 bufline++;
3104                 if (bufline == 3)
3105                     vlineSmooth3N(ptrs, ky + kylen - bufline, i + 1, dst + idst*dst_stride, width*cn);
3106                 else if (bufline == 5)
3107                     vlineSmooth5N(ptrs, ky + kylen - bufline, i + 1, dst + idst*dst_stride, width*cn);
3108                 else
3109                     vlineSmooth(ptrs, ky + kylen - bufline, i + 1, dst + idst*dst_stride, width*cn);
3110                 bufline %= kylen;
3111             }
3112             // Points inside the border
3113             if (i - max(0, ifrom) >= kylen)
3114             {
3115                 for (; i < ito; i++, idst++)
3116                 {
3117                     hlineSmoothFunc(src + i * src_stride, cn, kx, kxlen, ptrs[bufline], width, borderType);
3118                     bufline = (bufline + 1) % kylen;
3119                     vlineSmoothFunc(ptrs + bufline, ky, kylen, dst + idst*dst_stride, width*cn);
3120                 }
3121
3122                 // Points that could fall below border
3123                 // i points to first src row to evaluate in convolution
3124                 bufline = (bufline + 1) % kylen;
3125                 for (i = idst - pre_shift; i < range.end - pre_shift; i++, idst++, bufline++)
3126                     if (height - i == 3)
3127                         vlineSmooth3N(ptrs + bufline, ky, height - i, dst + idst*dst_stride, width*cn);
3128                     else if (height - i == 5)
3129                         vlineSmooth5N(ptrs + bufline, ky, height - i, dst + idst*dst_stride, width*cn);
3130                     else
3131                         vlineSmooth(ptrs + bufline, ky, height - i, dst + idst*dst_stride, width*cn);
3132             }
3133             else
3134             {
3135                 // i points to first src row to evaluate in convolution
3136                 for (i = idst - pre_shift; i < min(range.end - pre_shift, 0); i++, idst++)
3137                     if (height == 3)
3138                         vlineSmooth3N(ptrs, ky - i, height, dst + idst*dst_stride, width*cn);
3139                     else if (height == 5)
3140                         vlineSmooth5N(ptrs, ky - i, height, dst + idst*dst_stride, width*cn);
3141                     else
3142                         vlineSmooth(ptrs, ky - i, height, dst + idst*dst_stride, width*cn);
3143                 for (; i < range.end - pre_shift; i++, idst++)
3144                     if (height - i == 3)
3145                         vlineSmooth3N(ptrs + i - max(0, ifrom), ky, height - i, dst + idst*dst_stride, width*cn);
3146                     else if (height - i == 5)
3147                         vlineSmooth5N(ptrs + i - max(0, ifrom), ky, height - i, dst + idst*dst_stride, width*cn);
3148                     else
3149                         vlineSmooth(ptrs + i - max(0, ifrom), ky, height - i, dst + idst*dst_stride, width*cn);
3150             }
3151         }
3152     }
3153 private:
3154     const ET* src;
3155     ET* dst;
3156     size_t src_stride, dst_stride;
3157     int width, height, cn;
3158     const FT *kx, *ky;
3159     int kxlen, kylen;
3160     int borderType;
3161     void(*hlineSmoothFunc)(const ET* src, int cn, const FT* m, int n, FT* dst, int len, int borderType);
3162     void(*vlineSmoothFunc)(const FT* const * src, const FT* m, int n, ET* dst, int len);
3163
3164     fixedSmoothInvoker(const fixedSmoothInvoker&);
3165     fixedSmoothInvoker& operator=(const fixedSmoothInvoker&);
3166 };
3167
3168 static void getGaussianKernel(int n, double sigma, int ktype, Mat& res) { res = getGaussianKernel(n, sigma, ktype); }
3169 template <typename T> static void getGaussianKernel(int n, double sigma, int, std::vector<T>& res) { res = getFixedpointGaussianKernel<T>(n, sigma); }
3170
3171 template <typename T>
3172 static void createGaussianKernels( T & kx, T & ky, int type, Size &ksize,
3173                                    double sigma1, double sigma2 )
3174 {
3175     int depth = CV_MAT_DEPTH(type);
3176     if( sigma2 <= 0 )
3177         sigma2 = sigma1;
3178
3179     // automatic detection of kernel size from sigma
3180     if( ksize.width <= 0 && sigma1 > 0 )
3181         ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
3182     if( ksize.height <= 0 && sigma2 > 0 )
3183         ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
3184
3185     CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 &&
3186         ksize.height > 0 && ksize.height % 2 == 1 );
3187
3188     sigma1 = std::max( sigma1, 0. );
3189     sigma2 = std::max( sigma2, 0. );
3190
3191     getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F), kx );
3192     if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )
3193         ky = kx;
3194     else
3195         getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F), ky );
3196 }
3197
3198 }
3199
3200 cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize,
3201                                         double sigma1, double sigma2,
3202                                         int borderType )
3203 {
3204     Mat kx, ky;
3205     createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
3206
3207     return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );
3208 }
3209
3210 namespace cv
3211 {
3212 #ifdef HAVE_OPENCL
3213
3214 static bool ocl_GaussianBlur_8UC1(InputArray _src, OutputArray _dst, Size ksize, int ddepth,
3215                                   InputArray _kernelX, InputArray _kernelY, int borderType)
3216 {
3217     const ocl::Device & dev = ocl::Device::getDefault();
3218     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
3219
3220     if ( !(dev.isIntel() && (type == CV_8UC1) &&
3221          (_src.offset() == 0) && (_src.step() % 4 == 0) &&
3222          ((ksize.width == 5 && (_src.cols() % 4 == 0)) ||
3223          (ksize.width == 3 && (_src.cols() % 16 == 0) && (_src.rows() % 2 == 0)))) )
3224         return false;
3225
3226     Mat kernelX = _kernelX.getMat().reshape(1, 1);
3227     if (kernelX.cols % 2 != 1)
3228         return false;
3229     Mat kernelY = _kernelY.getMat().reshape(1, 1);
3230     if (kernelY.cols % 2 != 1)
3231         return false;
3232
3233     if (ddepth < 0)
3234         ddepth = sdepth;
3235
3236     Size size = _src.size();
3237     size_t globalsize[2] = { 0, 0 };
3238     size_t localsize[2] = { 0, 0 };
3239
3240     if (ksize.width == 3)
3241     {
3242         globalsize[0] = size.width / 16;
3243         globalsize[1] = size.height / 2;
3244     }
3245     else if (ksize.width == 5)
3246     {
3247         globalsize[0] = size.width / 4;
3248         globalsize[1] = size.height / 1;
3249     }
3250
3251     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
3252     char build_opts[1024];
3253     sprintf(build_opts, "-D %s %s%s", borderMap[borderType & ~BORDER_ISOLATED],
3254             ocl::kernelToStr(kernelX, CV_32F, "KERNEL_MATRIX_X").c_str(),
3255             ocl::kernelToStr(kernelY, CV_32F, "KERNEL_MATRIX_Y").c_str());
3256
3257     ocl::Kernel kernel;
3258
3259     if (ksize.width == 3)
3260         kernel.create("gaussianBlur3x3_8UC1_cols16_rows2", cv::ocl::imgproc::gaussianBlur3x3_oclsrc, build_opts);
3261     else if (ksize.width == 5)
3262         kernel.create("gaussianBlur5x5_8UC1_cols4", cv::ocl::imgproc::gaussianBlur5x5_oclsrc, build_opts);
3263
3264     if (kernel.empty())
3265         return false;
3266
3267     UMat src = _src.getUMat();
3268     _dst.create(size, CV_MAKETYPE(ddepth, cn));
3269     if (!(_dst.offset() == 0 && _dst.step() % 4 == 0))
3270         return false;
3271     UMat dst = _dst.getUMat();
3272
3273     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
3274     idxArg = kernel.set(idxArg, (int)src.step);
3275     idxArg = kernel.set(idxArg, ocl::KernelArg::PtrWriteOnly(dst));
3276     idxArg = kernel.set(idxArg, (int)dst.step);
3277     idxArg = kernel.set(idxArg, (int)dst.rows);
3278     idxArg = kernel.set(idxArg, (int)dst.cols);
3279
3280     return kernel.run(2, globalsize, (localsize[0] == 0) ? NULL : localsize, false);
3281 }
3282
3283 #endif
3284
3285 #ifdef HAVE_OPENVX
3286
3287 namespace ovx {
3288     template <> inline bool skipSmallImages<VX_KERNEL_GAUSSIAN_3x3>(int w, int h) { return w*h < 320 * 240; }
3289 }
3290 static bool openvx_gaussianBlur(InputArray _src, OutputArray _dst, Size ksize,
3291                                 double sigma1, double sigma2, int borderType)
3292 {
3293     if (sigma2 <= 0)
3294         sigma2 = sigma1;
3295     // automatic detection of kernel size from sigma
3296     if (ksize.width <= 0 && sigma1 > 0)
3297         ksize.width = cvRound(sigma1*6 + 1) | 1;
3298     if (ksize.height <= 0 && sigma2 > 0)
3299         ksize.height = cvRound(sigma2*6 + 1) | 1;
3300
3301     if (_src.type() != CV_8UC1 ||
3302         _src.cols() < 3 || _src.rows() < 3 ||
3303         ksize.width != 3 || ksize.height != 3)
3304         return false;
3305
3306     sigma1 = std::max(sigma1, 0.);
3307     sigma2 = std::max(sigma2, 0.);
3308
3309     if (!(sigma1 == 0.0 || (sigma1 - 0.8) < DBL_EPSILON) || !(sigma2 == 0.0 || (sigma2 - 0.8) < DBL_EPSILON) ||
3310         ovx::skipSmallImages<VX_KERNEL_GAUSSIAN_3x3>(_src.cols(), _src.rows()))
3311         return false;
3312
3313     Mat src = _src.getMat();
3314     Mat dst = _dst.getMat();
3315
3316     if ((borderType & BORDER_ISOLATED) == 0 && src.isSubmatrix())
3317         return false; //Process isolated borders only
3318     vx_enum border;
3319     switch (borderType & ~BORDER_ISOLATED)
3320     {
3321     case BORDER_CONSTANT:
3322         border = VX_BORDER_CONSTANT;
3323         break;
3324     case BORDER_REPLICATE:
3325         border = VX_BORDER_REPLICATE;
3326         break;
3327     default:
3328         return false;
3329     }
3330
3331     try
3332     {
3333         ivx::Context ctx = ovx::getOpenVXContext();
3334
3335         Mat a;
3336         if (dst.data != src.data)
3337             a = src;
3338         else
3339             src.copyTo(a);
3340
3341         ivx::Image
3342             ia = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
3343                 ivx::Image::createAddressing(a.cols, a.rows, 1, (vx_int32)(a.step)), a.data),
3344             ib = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
3345                 ivx::Image::createAddressing(dst.cols, dst.rows, 1, (vx_int32)(dst.step)), dst.data);
3346
3347         //ATTENTION: VX_CONTEXT_IMMEDIATE_BORDER attribute change could lead to strange issues in multi-threaded environments
3348         //since OpenVX standard says nothing about thread-safety for now
3349         ivx::border_t prevBorder = ctx.immediateBorder();
3350         ctx.setImmediateBorder(border, (vx_uint8)(0));
3351         ivx::IVX_CHECK_STATUS(vxuGaussian3x3(ctx, ia, ib));
3352         ctx.setImmediateBorder(prevBorder);
3353     }
3354     catch (ivx::RuntimeError & e)
3355     {
3356         VX_DbgThrow(e.what());
3357     }
3358     catch (ivx::WrapperError & e)
3359     {
3360         VX_DbgThrow(e.what());
3361     }
3362     return true;
3363 }
3364
3365 #endif
3366
3367 #ifdef HAVE_IPP
3368 #if IPP_VERSION_X100 == 201702  // IW 2017u2 has bug which doesn't allow use of partial inMem with tiling
3369 #define IPP_GAUSSIANBLUR_PARALLEL 0
3370 #else
3371 #define IPP_GAUSSIANBLUR_PARALLEL 1
3372 #endif
3373
3374 #ifdef HAVE_IPP_IW
3375
3376 class ipp_gaussianBlurParallel: public ParallelLoopBody
3377 {
3378 public:
3379     ipp_gaussianBlurParallel(::ipp::IwiImage &src, ::ipp::IwiImage &dst, int kernelSize, float sigma, ::ipp::IwiBorderType &border, bool *pOk):
3380         m_src(src), m_dst(dst), m_kernelSize(kernelSize), m_sigma(sigma), m_border(border), m_pOk(pOk) {
3381         *m_pOk = true;
3382     }
3383     ~ipp_gaussianBlurParallel()
3384     {
3385     }
3386
3387     virtual void operator() (const Range& range) const
3388     {
3389         CV_INSTRUMENT_REGION_IPP()
3390
3391         if(!*m_pOk)
3392             return;
3393
3394         try
3395         {
3396             ::ipp::IwiTile tile = ::ipp::IwiRoi(0, range.start, m_dst.m_size.width, range.end - range.start);
3397             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterGaussian, m_src, m_dst, m_kernelSize, m_sigma, ::ipp::IwDefault(), m_border, tile);
3398         }
3399         catch(::ipp::IwException e)
3400         {
3401             *m_pOk = false;
3402             return;
3403         }
3404     }
3405 private:
3406     ::ipp::IwiImage &m_src;
3407     ::ipp::IwiImage &m_dst;
3408
3409     int m_kernelSize;
3410     float m_sigma;
3411     ::ipp::IwiBorderType &m_border;
3412
3413     volatile bool *m_pOk;
3414     const ipp_gaussianBlurParallel& operator= (const ipp_gaussianBlurParallel&);
3415 };
3416
3417 #endif
3418
3419 static bool ipp_GaussianBlur(InputArray _src, OutputArray _dst, Size ksize,
3420                    double sigma1, double sigma2, int borderType )
3421 {
3422 #ifdef HAVE_IPP_IW
3423     CV_INSTRUMENT_REGION_IPP()
3424
3425 #if IPP_VERSION_X100 < 201800 && ((defined _MSC_VER && defined _M_IX86) || (defined __GNUC__ && defined __i386__))
3426     CV_UNUSED(_src); CV_UNUSED(_dst); CV_UNUSED(ksize); CV_UNUSED(sigma1); CV_UNUSED(sigma2); CV_UNUSED(borderType);
3427     return false; // bug on ia32
3428 #else
3429     if(sigma1 != sigma2)
3430         return false;
3431
3432     if(sigma1 < FLT_EPSILON)
3433         return false;
3434
3435     if(ksize.width != ksize.height)
3436         return false;
3437
3438     // Acquire data and begin processing
3439     try
3440     {
3441         Mat src = _src.getMat();
3442         Mat dst = _dst.getMat();
3443         ::ipp::IwiImage       iwSrc      = ippiGetImage(src);
3444         ::ipp::IwiImage       iwDst      = ippiGetImage(dst);
3445         ::ipp::IwiBorderSize  borderSize = ::ipp::iwiSizeToBorderSize(ippiGetSize(ksize));
3446         ::ipp::IwiBorderType  ippBorder(ippiGetBorder(iwSrc, borderType, borderSize));
3447         if(!ippBorder)
3448             return false;
3449
3450         const int threads = ippiSuggestThreadsNum(iwDst, 2);
3451         if(IPP_GAUSSIANBLUR_PARALLEL && threads > 1) {
3452             bool ok;
3453             ipp_gaussianBlurParallel invoker(iwSrc, iwDst, ksize.width, (float) sigma1, ippBorder, &ok);
3454
3455             if(!ok)
3456                 return false;
3457             const Range range(0, (int) iwDst.m_size.height);
3458             parallel_for_(range, invoker, threads*4);
3459
3460             if(!ok)
3461                 return false;
3462         } else {
3463             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterGaussian, iwSrc, iwDst, ksize.width, sigma1, ::ipp::IwDefault(), ippBorder);
3464         }
3465     }
3466     catch (::ipp::IwException ex)
3467     {
3468         return false;
3469     }
3470
3471     return true;
3472 #endif
3473 #else
3474     CV_UNUSED(_src); CV_UNUSED(_dst); CV_UNUSED(ksize); CV_UNUSED(sigma1); CV_UNUSED(sigma2); CV_UNUSED(borderType);
3475     return false;
3476 #endif
3477 }
3478 #endif
3479 }
3480
3481 void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
3482                    double sigma1, double sigma2,
3483                    int borderType )
3484 {
3485     CV_INSTRUMENT_REGION()
3486
3487     int type = _src.type();
3488     Size size = _src.size();
3489     _dst.create( size, type );
3490
3491     if( (borderType & ~BORDER_ISOLATED) != BORDER_CONSTANT &&
3492         ((borderType & BORDER_ISOLATED) != 0 || !_src.getMat().isSubmatrix()) )
3493     {
3494         if( size.height == 1 )
3495             ksize.height = 1;
3496         if( size.width == 1 )
3497             ksize.width = 1;
3498     }
3499
3500     if( ksize.width == 1 && ksize.height == 1 )
3501     {
3502         _src.copyTo(_dst);
3503         return;
3504     }
3505
3506     bool useOpenCL = (ocl::isOpenCLActivated() && _dst.isUMat() && _src.dims() <= 2 &&
3507                ((ksize.width == 3 && ksize.height == 3) ||
3508                (ksize.width == 5 && ksize.height == 5)) &&
3509                _src.rows() > ksize.height && _src.cols() > ksize.width);
3510     (void)useOpenCL;
3511
3512     int sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
3513
3514     if(sdepth == CV_8U && ((borderType & BORDER_ISOLATED) || !_src.getMat().isSubmatrix()))
3515     {
3516         std::vector<ufixedpoint16> fkx, fky;
3517         createGaussianKernels(fkx, fky, type, ksize, sigma1, sigma2);
3518         Mat src = _src.getMat();
3519         Mat dst = _dst.getMat();
3520         if (src.data == dst.data)
3521             src = src.clone();
3522         fixedSmoothInvoker<uint8_t, ufixedpoint16> invoker(src.ptr<uint8_t>(), src.step1(), dst.ptr<uint8_t>(), dst.step1(), dst.cols, dst.rows, dst.channels(), &fkx[0], (int)fkx.size(), &fky[0], (int)fky.size(), borderType & ~BORDER_ISOLATED);
3523         parallel_for_(Range(0, dst.rows), invoker, dst.total() * cn / (double)(1 << 13));
3524         return;
3525     }
3526
3527
3528     Mat kx, ky;
3529     createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
3530
3531     CV_OCL_RUN(useOpenCL, ocl_GaussianBlur_8UC1(_src, _dst, ksize, CV_MAT_DEPTH(type), kx, ky, borderType));
3532
3533     CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2 && (size_t)_src.rows() > kx.total() && (size_t)_src.cols() > kx.total(),
3534                ocl_sepFilter2D(_src, _dst, sdepth, kx, ky, Point(-1, -1), 0, borderType))
3535
3536     Mat src = _src.getMat();
3537     Mat dst = _dst.getMat();
3538
3539     Point ofs;
3540     Size wsz(src.cols, src.rows);
3541     if(!(borderType & BORDER_ISOLATED))
3542         src.locateROI( wsz, ofs );
3543
3544     CALL_HAL(gaussianBlur, cv_hal_gaussianBlur, src.ptr(), src.step, dst.ptr(), dst.step, src.cols, src.rows, sdepth, cn,
3545              ofs.x, ofs.y, wsz.width - src.cols - ofs.x, wsz.height - src.rows - ofs.y, ksize.width, ksize.height,
3546              sigma1, sigma2, borderType&~BORDER_ISOLATED);
3547
3548     CV_OVX_RUN(true,
3549                openvx_gaussianBlur(src, dst, ksize, sigma1, sigma2, borderType))
3550
3551     CV_IPP_RUN_FAST(ipp_GaussianBlur(src, dst, ksize, sigma1, sigma2, borderType));
3552
3553     sepFilter2D(src, dst, sdepth, kx, ky, Point(-1, -1), 0, borderType);
3554 }
3555
3556 /****************************************************************************************\
3557                                       Median Filter
3558 \****************************************************************************************/
3559
3560 namespace cv
3561 {
3562 typedef ushort HT;
3563
3564 /**
3565  * This structure represents a two-tier histogram. The first tier (known as the
3566  * "coarse" level) is 4 bit wide and the second tier (known as the "fine" level)
3567  * is 8 bit wide. Pixels inserted in the fine level also get inserted into the
3568  * coarse bucket designated by the 4 MSBs of the fine bucket value.
3569  *
3570  * The structure is aligned on 16 bits, which is a prerequisite for SIMD
3571  * instructions. Each bucket is 16 bit wide, which means that extra care must be
3572  * taken to prevent overflow.
3573  */
3574 typedef struct
3575 {
3576     HT coarse[16];
3577     HT fine[16][16];
3578 } Histogram;
3579
3580
3581 #if CV_SIMD128
3582
3583 static inline void histogram_add_simd( const HT x[16], HT y[16] )
3584 {
3585     v_store(y, v_load(x) + v_load(y));
3586     v_store(y + 8, v_load(x + 8) + v_load(y + 8));
3587 }
3588
3589 static inline void histogram_sub_simd( const HT x[16], HT y[16] )
3590 {
3591     v_store(y, v_load(y) - v_load(x));
3592     v_store(y + 8, v_load(y + 8) - v_load(x + 8));
3593 }
3594
3595 #endif
3596
3597
3598 static inline void histogram_add( const HT x[16], HT y[16] )
3599 {
3600     int i;
3601     for( i = 0; i < 16; ++i )
3602         y[i] = (HT)(y[i] + x[i]);
3603 }
3604
3605 static inline void histogram_sub( const HT x[16], HT y[16] )
3606 {
3607     int i;
3608     for( i = 0; i < 16; ++i )
3609         y[i] = (HT)(y[i] - x[i]);
3610 }
3611
3612 static inline void histogram_muladd( int a, const HT x[16],
3613         HT y[16] )
3614 {
3615     for( int i = 0; i < 16; ++i )
3616         y[i] = (HT)(y[i] + a * x[i]);
3617 }
3618
3619 static void
3620 medianBlur_8u_O1( const Mat& _src, Mat& _dst, int ksize )
3621 {
3622 /**
3623  * HOP is short for Histogram OPeration. This macro makes an operation \a op on
3624  * histogram \a h for pixel value \a x. It takes care of handling both levels.
3625  */
3626 #define HOP(h,x,op) \
3627     h.coarse[x>>4] op, \
3628     *((HT*)h.fine + x) op
3629
3630 #define COP(c,j,x,op) \
3631     h_coarse[ 16*(n*c+j) + (x>>4) ] op, \
3632     h_fine[ 16 * (n*(16*c+(x>>4)) + j) + (x & 0xF) ] op
3633
3634     int cn = _dst.channels(), m = _dst.rows, r = (ksize-1)/2;
3635     CV_Assert(cn > 0 && cn <= 4);
3636     size_t sstep = _src.step, dstep = _dst.step;
3637     Histogram CV_DECL_ALIGNED(16) H[4];
3638     HT CV_DECL_ALIGNED(16) luc[4][16];
3639
3640     int STRIPE_SIZE = std::min( _dst.cols, 512/cn );
3641
3642     std::vector<HT> _h_coarse(1 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
3643     std::vector<HT> _h_fine(16 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
3644     HT* h_coarse = alignPtr(&_h_coarse[0], 16);
3645     HT* h_fine = alignPtr(&_h_fine[0], 16);
3646 #if CV_SIMD128
3647     volatile bool useSIMD = hasSIMD128();
3648 #endif
3649
3650     for( int x = 0; x < _dst.cols; x += STRIPE_SIZE )
3651     {
3652         int i, j, k, c, n = std::min(_dst.cols - x, STRIPE_SIZE) + r*2;
3653         const uchar* src = _src.ptr() + x*cn;
3654         uchar* dst = _dst.ptr() + (x - r)*cn;
3655
3656         memset( h_coarse, 0, 16*n*cn*sizeof(h_coarse[0]) );
3657         memset( h_fine, 0, 16*16*n*cn*sizeof(h_fine[0]) );
3658
3659         // First row initialization
3660         for( c = 0; c < cn; c++ )
3661         {
3662             for( j = 0; j < n; j++ )
3663                 COP( c, j, src[cn*j+c], += (cv::HT)(r+2) );
3664
3665             for( i = 1; i < r; i++ )
3666             {
3667                 const uchar* p = src + sstep*std::min(i, m-1);
3668                 for ( j = 0; j < n; j++ )
3669                     COP( c, j, p[cn*j+c], ++ );
3670             }
3671         }
3672
3673         for( i = 0; i < m; i++ )
3674         {
3675             const uchar* p0 = src + sstep * std::max( 0, i-r-1 );
3676             const uchar* p1 = src + sstep * std::min( m-1, i+r );
3677
3678             memset( H, 0, cn*sizeof(H[0]) );
3679             memset( luc, 0, cn*sizeof(luc[0]) );
3680             for( c = 0; c < cn; c++ )
3681             {
3682                 // Update column histograms for the entire row.
3683                 for( j = 0; j < n; j++ )
3684                 {
3685                     COP( c, j, p0[j*cn + c], -- );
3686                     COP( c, j, p1[j*cn + c], ++ );
3687                 }
3688
3689                 // First column initialization
3690                 for( k = 0; k < 16; ++k )
3691                     histogram_muladd( 2*r+1, &h_fine[16*n*(16*c+k)], &H[c].fine[k][0] );
3692
3693 #if CV_SIMD128
3694                 if( useSIMD )
3695                 {
3696                     for( j = 0; j < 2*r; ++j )
3697                         histogram_add_simd( &h_coarse[16*(n*c+j)], H[c].coarse );
3698
3699                     for( j = r; j < n-r; j++ )
3700                     {
3701                         int t = 2*r*r + 2*r, b, sum = 0;
3702                         HT* segment;
3703
3704                         histogram_add_simd( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
3705
3706                         // Find median at coarse level
3707                         for ( k = 0; k < 16 ; ++k )
3708                         {
3709                             sum += H[c].coarse[k];
3710                             if ( sum > t )
3711                             {
3712                                 sum -= H[c].coarse[k];
3713                                 break;
3714                             }
3715                         }
3716                         CV_Assert( k < 16 );
3717
3718                         /* Update corresponding histogram segment */
3719                         if ( luc[c][k] <= j-r )
3720                         {
3721                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
3722                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
3723                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
3724
3725                             if ( luc[c][k] < j+r+1 )
3726                             {
3727                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
3728                                 luc[c][k] = (HT)(j+r+1);
3729                             }
3730                         }
3731                         else
3732                         {
3733                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
3734                             {
3735                                 histogram_sub_simd( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
3736                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
3737                             }
3738                         }
3739
3740                         histogram_sub_simd( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
3741
3742                         /* Find median in segment */
3743                         segment = H[c].fine[k];
3744                         for ( b = 0; b < 16 ; b++ )
3745                         {
3746                             sum += segment[b];
3747                             if ( sum > t )
3748                             {
3749                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
3750                                 break;
3751                             }
3752                         }
3753                         CV_Assert( b < 16 );
3754                     }
3755                 }
3756                 else
3757 #endif
3758                 {
3759                     for( j = 0; j < 2*r; ++j )
3760                         histogram_add( &h_coarse[16*(n*c+j)], H[c].coarse );
3761
3762                     for( j = r; j < n-r; j++ )
3763                     {
3764                         int t = 2*r*r + 2*r, b, sum = 0;
3765                         HT* segment;
3766
3767                         histogram_add( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
3768
3769                         // Find median at coarse level
3770                         for ( k = 0; k < 16 ; ++k )
3771                         {
3772                             sum += H[c].coarse[k];
3773                             if ( sum > t )
3774                             {
3775                                 sum -= H[c].coarse[k];
3776                                 break;
3777                             }
3778                         }
3779                         CV_Assert( k < 16 );
3780
3781                         /* Update corresponding histogram segment */
3782                         if ( luc[c][k] <= j-r )
3783                         {
3784                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
3785                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
3786                                 histogram_add( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
3787
3788                             if ( luc[c][k] < j+r+1 )
3789                             {
3790                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
3791                                 luc[c][k] = (HT)(j+r+1);
3792                             }
3793                         }
3794                         else
3795                         {
3796                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
3797                             {
3798                                 histogram_sub( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
3799                                 histogram_add( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
3800                             }
3801                         }
3802
3803                         histogram_sub( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
3804
3805                         /* Find median in segment */
3806                         segment = H[c].fine[k];
3807                         for ( b = 0; b < 16 ; b++ )
3808                         {
3809                             sum += segment[b];
3810                             if ( sum > t )
3811                             {
3812                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
3813                                 break;
3814                             }
3815                         }
3816                         CV_Assert( b < 16 );
3817                     }
3818                 }
3819             }
3820         }
3821     }
3822
3823 #undef HOP
3824 #undef COP
3825 }
3826
3827 static void
3828 medianBlur_8u_Om( const Mat& _src, Mat& _dst, int m )
3829 {
3830     #define N  16
3831     int     zone0[4][N];
3832     int     zone1[4][N*N];
3833     int     x, y;
3834     int     n2 = m*m/2;
3835     Size    size = _dst.size();
3836     const uchar* src = _src.ptr();
3837     uchar*  dst = _dst.ptr();
3838     int     src_step = (int)_src.step, dst_step = (int)_dst.step;
3839     int     cn = _src.channels();
3840     const uchar*  src_max = src + size.height*src_step;
3841     CV_Assert(cn > 0 && cn <= 4);
3842
3843     #define UPDATE_ACC01( pix, cn, op ) \
3844     {                                   \
3845         int p = (pix);                  \
3846         zone1[cn][p] op;                \
3847         zone0[cn][p >> 4] op;           \
3848     }
3849
3850     //CV_Assert( size.height >= nx && size.width >= nx );
3851     for( x = 0; x < size.width; x++, src += cn, dst += cn )
3852     {
3853         uchar* dst_cur = dst;
3854         const uchar* src_top = src;
3855         const uchar* src_bottom = src;
3856         int k, c;
3857         int src_step1 = src_step, dst_step1 = dst_step;
3858
3859         if( x % 2 != 0 )
3860         {
3861             src_bottom = src_top += src_step*(size.height-1);
3862             dst_cur += dst_step*(size.height-1);
3863             src_step1 = -src_step1;
3864             dst_step1 = -dst_step1;
3865         }
3866
3867         // init accumulator
3868         memset( zone0, 0, sizeof(zone0[0])*cn );
3869         memset( zone1, 0, sizeof(zone1[0])*cn );
3870
3871         for( y = 0; y <= m/2; y++ )
3872         {
3873             for( c = 0; c < cn; c++ )
3874             {
3875                 if( y > 0 )
3876                 {
3877                     for( k = 0; k < m*cn; k += cn )
3878                         UPDATE_ACC01( src_bottom[k+c], c, ++ );
3879                 }
3880                 else
3881                 {
3882                     for( k = 0; k < m*cn; k += cn )
3883                         UPDATE_ACC01( src_bottom[k+c], c, += m/2+1 );
3884                 }
3885             }
3886
3887             if( (src_step1 > 0 && y < size.height-1) ||
3888                 (src_step1 < 0 && size.height-y-1 > 0) )
3889                 src_bottom += src_step1;
3890         }
3891
3892         for( y = 0; y < size.height; y++, dst_cur += dst_step1 )
3893         {
3894             // find median
3895             for( c = 0; c < cn; c++ )
3896             {
3897                 int s = 0;
3898                 for( k = 0; ; k++ )
3899                 {
3900                     int t = s + zone0[c][k];
3901                     if( t > n2 ) break;
3902                     s = t;
3903                 }
3904
3905                 for( k *= N; ;k++ )
3906                 {
3907                     s += zone1[c][k];
3908                     if( s > n2 ) break;
3909                 }
3910
3911                 dst_cur[c] = (uchar)k;
3912             }
3913
3914             if( y+1 == size.height )
3915                 break;
3916
3917             if( cn == 1 )
3918             {
3919                 for( k = 0; k < m; k++ )
3920                 {
3921                     int p = src_top[k];
3922                     int q = src_bottom[k];
3923                     zone1[0][p]--;
3924                     zone0[0][p>>4]--;
3925                     zone1[0][q]++;
3926                     zone0[0][q>>4]++;
3927                 }
3928             }
3929             else if( cn == 3 )
3930             {
3931                 for( k = 0; k < m*3; k += 3 )
3932                 {
3933                     UPDATE_ACC01( src_top[k], 0, -- );
3934                     UPDATE_ACC01( src_top[k+1], 1, -- );
3935                     UPDATE_ACC01( src_top[k+2], 2, -- );
3936
3937                     UPDATE_ACC01( src_bottom[k], 0, ++ );
3938                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
3939                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
3940                 }
3941             }
3942             else
3943             {
3944                 assert( cn == 4 );
3945                 for( k = 0; k < m*4; k += 4 )
3946                 {
3947                     UPDATE_ACC01( src_top[k], 0, -- );
3948                     UPDATE_ACC01( src_top[k+1], 1, -- );
3949                     UPDATE_ACC01( src_top[k+2], 2, -- );
3950                     UPDATE_ACC01( src_top[k+3], 3, -- );
3951
3952                     UPDATE_ACC01( src_bottom[k], 0, ++ );
3953                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
3954                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
3955                     UPDATE_ACC01( src_bottom[k+3], 3, ++ );
3956                 }
3957             }
3958
3959             if( (src_step1 > 0 && src_bottom + src_step1 < src_max) ||
3960                 (src_step1 < 0 && src_bottom + src_step1 >= src) )
3961                 src_bottom += src_step1;
3962
3963             if( y >= m/2 )
3964                 src_top += src_step1;
3965         }
3966     }
3967 #undef N
3968 #undef UPDATE_ACC
3969 }
3970
3971
3972 struct MinMax8u
3973 {
3974     typedef uchar value_type;
3975     typedef int arg_type;
3976     enum { SIZE = 1 };
3977     arg_type load(const uchar* ptr) { return *ptr; }
3978     void store(uchar* ptr, arg_type val) { *ptr = (uchar)val; }
3979     void operator()(arg_type& a, arg_type& b) const
3980     {
3981         int t = CV_FAST_CAST_8U(a - b);
3982         b += t; a -= t;
3983     }
3984 };
3985
3986 struct MinMax16u
3987 {
3988     typedef ushort value_type;
3989     typedef int arg_type;
3990     enum { SIZE = 1 };
3991     arg_type load(const ushort* ptr) { return *ptr; }
3992     void store(ushort* ptr, arg_type val) { *ptr = (ushort)val; }
3993     void operator()(arg_type& a, arg_type& b) const
3994     {
3995         arg_type t = a;
3996         a = std::min(a, b);
3997         b = std::max(b, t);
3998     }
3999 };
4000
4001 struct MinMax16s
4002 {
4003     typedef short value_type;
4004     typedef int arg_type;
4005     enum { SIZE = 1 };
4006     arg_type load(const short* ptr) { return *ptr; }
4007     void store(short* ptr, arg_type val) { *ptr = (short)val; }
4008     void operator()(arg_type& a, arg_type& b) const
4009     {
4010         arg_type t = a;
4011         a = std::min(a, b);
4012         b = std::max(b, t);
4013     }
4014 };
4015
4016 struct MinMax32f
4017 {
4018     typedef float value_type;
4019     typedef float arg_type;
4020     enum { SIZE = 1 };
4021     arg_type load(const float* ptr) { return *ptr; }
4022     void store(float* ptr, arg_type val) { *ptr = val; }
4023     void operator()(arg_type& a, arg_type& b) const
4024     {
4025         arg_type t = a;
4026         a = std::min(a, b);
4027         b = std::max(b, t);
4028     }
4029 };
4030
4031 #if CV_SIMD128
4032
4033 struct MinMaxVec8u
4034 {
4035     typedef uchar value_type;
4036     typedef v_uint8x16 arg_type;
4037     enum { SIZE = 16 };
4038     arg_type load(const uchar* ptr) { return v_load(ptr); }
4039     void store(uchar* ptr, const arg_type &val) { v_store(ptr, val); }
4040     void operator()(arg_type& a, arg_type& b) const
4041     {
4042         arg_type t = a;
4043         a = v_min(a, b);
4044         b = v_max(b, t);
4045     }
4046 };
4047
4048
4049 struct MinMaxVec16u
4050 {
4051     typedef ushort value_type;
4052     typedef v_uint16x8 arg_type;
4053     enum { SIZE = 8 };
4054     arg_type load(const ushort* ptr) { return v_load(ptr); }
4055     void store(ushort* ptr, const arg_type &val) { v_store(ptr, val); }
4056     void operator()(arg_type& a, arg_type& b) const
4057     {
4058         arg_type t = a;
4059         a = v_min(a, b);
4060         b = v_max(b, t);
4061     }
4062 };
4063
4064
4065 struct MinMaxVec16s
4066 {
4067     typedef short value_type;
4068     typedef v_int16x8 arg_type;
4069     enum { SIZE = 8 };
4070     arg_type load(const short* ptr) { return v_load(ptr); }
4071     void store(short* ptr, const arg_type &val) { v_store(ptr, val); }
4072     void operator()(arg_type& a, arg_type& b) const
4073     {
4074         arg_type t = a;
4075         a = v_min(a, b);
4076         b = v_max(b, t);
4077     }
4078 };
4079
4080
4081 struct MinMaxVec32f
4082 {
4083     typedef float value_type;
4084     typedef v_float32x4 arg_type;
4085     enum { SIZE = 4 };
4086     arg_type load(const float* ptr) { return v_load(ptr); }
4087     void store(float* ptr, const arg_type &val) { v_store(ptr, val); }
4088     void operator()(arg_type& a, arg_type& b) const
4089     {
4090         arg_type t = a;
4091         a = v_min(a, b);
4092         b = v_max(b, t);
4093     }
4094 };
4095
4096 #else
4097
4098 typedef MinMax8u MinMaxVec8u;
4099 typedef MinMax16u MinMaxVec16u;
4100 typedef MinMax16s MinMaxVec16s;
4101 typedef MinMax32f MinMaxVec32f;
4102
4103 #endif
4104
4105 template<class Op, class VecOp>
4106 static void
4107 medianBlur_SortNet( const Mat& _src, Mat& _dst, int m )
4108 {
4109     typedef typename Op::value_type T;
4110     typedef typename Op::arg_type WT;
4111     typedef typename VecOp::arg_type VT;
4112
4113     const T* src = _src.ptr<T>();
4114     T* dst = _dst.ptr<T>();
4115     int sstep = (int)(_src.step/sizeof(T));
4116     int dstep = (int)(_dst.step/sizeof(T));
4117     Size size = _dst.size();
4118     int i, j, k, cn = _src.channels();
4119     Op op;
4120     VecOp vop;
4121     volatile bool useSIMD = hasSIMD128();
4122
4123     if( m == 3 )
4124     {
4125         if( size.width == 1 || size.height == 1 )
4126         {
4127             int len = size.width + size.height - 1;
4128             int sdelta = size.height == 1 ? cn : sstep;
4129             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
4130             int ddelta = size.height == 1 ? cn : dstep;
4131
4132             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
4133                 for( j = 0; j < cn; j++, src++ )
4134                 {
4135                     WT p0 = src[i > 0 ? -sdelta : 0];
4136                     WT p1 = src[0];
4137                     WT p2 = src[i < len - 1 ? sdelta : 0];
4138
4139                     op(p0, p1); op(p1, p2); op(p0, p1);
4140                     dst[j] = (T)p1;
4141                 }
4142             return;
4143         }
4144
4145         size.width *= cn;
4146         for( i = 0; i < size.height; i++, dst += dstep )
4147         {
4148             const T* row0 = src + std::max(i - 1, 0)*sstep;
4149             const T* row1 = src + i*sstep;
4150             const T* row2 = src + std::min(i + 1, size.height-1)*sstep;
4151             int limit = useSIMD ? cn : size.width;
4152
4153             for(j = 0;; )
4154             {
4155                 for( ; j < limit; j++ )
4156                 {
4157                     int j0 = j >= cn ? j - cn : j;
4158                     int j2 = j < size.width - cn ? j + cn : j;
4159                     WT p0 = row0[j0], p1 = row0[j], p2 = row0[j2];
4160                     WT p3 = row1[j0], p4 = row1[j], p5 = row1[j2];
4161                     WT p6 = row2[j0], p7 = row2[j], p8 = row2[j2];
4162
4163                     op(p1, p2); op(p4, p5); op(p7, p8); op(p0, p1);
4164                     op(p3, p4); op(p6, p7); op(p1, p2); op(p4, p5);
4165                     op(p7, p8); op(p0, p3); op(p5, p8); op(p4, p7);
4166                     op(p3, p6); op(p1, p4); op(p2, p5); op(p4, p7);
4167                     op(p4, p2); op(p6, p4); op(p4, p2);
4168                     dst[j] = (T)p4;
4169                 }
4170
4171                 if( limit == size.width )
4172                     break;
4173
4174                 for( ; j <= size.width - VecOp::SIZE - cn; j += VecOp::SIZE )
4175                 {
4176                     VT p0 = vop.load(row0+j-cn), p1 = vop.load(row0+j), p2 = vop.load(row0+j+cn);
4177                     VT p3 = vop.load(row1+j-cn), p4 = vop.load(row1+j), p5 = vop.load(row1+j+cn);
4178                     VT p6 = vop.load(row2+j-cn), p7 = vop.load(row2+j), p8 = vop.load(row2+j+cn);
4179
4180                     vop(p1, p2); vop(p4, p5); vop(p7, p8); vop(p0, p1);
4181                     vop(p3, p4); vop(p6, p7); vop(p1, p2); vop(p4, p5);
4182                     vop(p7, p8); vop(p0, p3); vop(p5, p8); vop(p4, p7);
4183                     vop(p3, p6); vop(p1, p4); vop(p2, p5); vop(p4, p7);
4184                     vop(p4, p2); vop(p6, p4); vop(p4, p2);
4185                     vop.store(dst+j, p4);
4186                 }
4187
4188                 limit = size.width;
4189             }
4190         }
4191     }
4192     else if( m == 5 )
4193     {
4194         if( size.width == 1 || size.height == 1 )
4195         {
4196             int len = size.width + size.height - 1;
4197             int sdelta = size.height == 1 ? cn : sstep;
4198             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
4199             int ddelta = size.height == 1 ? cn : dstep;
4200
4201             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
4202                 for( j = 0; j < cn; j++, src++ )
4203                 {
4204                     int i1 = i > 0 ? -sdelta : 0;
4205                     int i0 = i > 1 ? -sdelta*2 : i1;
4206                     int i3 = i < len-1 ? sdelta : 0;
4207                     int i4 = i < len-2 ? sdelta*2 : i3;
4208                     WT p0 = src[i0], p1 = src[i1], p2 = src[0], p3 = src[i3], p4 = src[i4];
4209
4210                     op(p0, p1); op(p3, p4); op(p2, p3); op(p3, p4); op(p0, p2);
4211                     op(p2, p4); op(p1, p3); op(p1, p2);
4212                     dst[j] = (T)p2;
4213                 }
4214             return;
4215         }
4216
4217         size.width *= cn;
4218         for( i = 0; i < size.height; i++, dst += dstep )
4219         {
4220             const T* row[5];
4221             row[0] = src + std::max(i - 2, 0)*sstep;
4222             row[1] = src + std::max(i - 1, 0)*sstep;
4223             row[2] = src + i*sstep;
4224             row[3] = src + std::min(i + 1, size.height-1)*sstep;
4225             row[4] = src + std::min(i + 2, size.height-1)*sstep;
4226             int limit = useSIMD ? cn*2 : size.width;
4227
4228             for(j = 0;; )
4229             {
4230                 for( ; j < limit; j++ )
4231                 {
4232                     WT p[25];
4233                     int j1 = j >= cn ? j - cn : j;
4234                     int j0 = j >= cn*2 ? j - cn*2 : j1;
4235                     int j3 = j < size.width - cn ? j + cn : j;
4236                     int j4 = j < size.width - cn*2 ? j + cn*2 : j3;
4237                     for( k = 0; k < 5; k++ )
4238                     {
4239                         const T* rowk = row[k];
4240                         p[k*5] = rowk[j0]; p[k*5+1] = rowk[j1];
4241                         p[k*5+2] = rowk[j]; p[k*5+3] = rowk[j3];
4242                         p[k*5+4] = rowk[j4];
4243                     }
4244
4245                     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]);
4246                     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]);
4247                     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]);
4248                     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]);
4249                     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]);
4250                     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]);
4251                     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]);
4252                     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]);
4253                     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]);
4254                     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]);
4255                     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]);
4256                     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]);
4257                     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]);
4258                     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]);
4259                     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]);
4260                     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]);
4261                     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]);
4262                     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]);
4263                     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]);
4264                     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]);
4265                     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]);
4266                     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]);
4267                     op(p[7], p[11]); op(p[11], p[13]); op(p[11], p[12]);
4268                     dst[j] = (T)p[12];
4269                 }
4270
4271                 if( limit == size.width )
4272                     break;
4273
4274                 for( ; j <= size.width - VecOp::SIZE - cn*2; j += VecOp::SIZE )
4275                 {
4276                     VT p[25];
4277                     for( k = 0; k < 5; k++ )
4278                     {
4279                         const T* rowk = row[k];
4280                         p[k*5] = vop.load(rowk+j-cn*2); p[k*5+1] = vop.load(rowk+j-cn);
4281                         p[k*5+2] = vop.load(rowk+j); p[k*5+3] = vop.load(rowk+j+cn);
4282                         p[k*5+4] = vop.load(rowk+j+cn*2);
4283                     }
4284
4285                     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]);
4286                     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]);
4287                     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]);
4288                     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]);
4289                     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]);
4290                     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]);
4291                     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]);
4292                     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]);
4293                     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]);
4294                     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]);
4295                     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]);
4296                     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]);
4297                     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]);
4298                     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]);
4299                     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]);
4300                     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]);
4301                     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]);
4302                     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]);
4303                     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]);
4304                     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]);
4305                     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]);
4306                     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]);
4307                     vop(p[7], p[11]); vop(p[11], p[13]); vop(p[11], p[12]);
4308                     vop.store(dst+j, p[12]);
4309                 }
4310
4311                 limit = size.width;
4312             }
4313         }
4314     }
4315 }
4316
4317 #ifdef HAVE_OPENCL
4318
4319 static bool ocl_medianFilter(InputArray _src, OutputArray _dst, int m)
4320 {
4321     size_t localsize[2] = { 16, 16 };
4322     size_t globalsize[2];
4323     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
4324
4325     if ( !((depth == CV_8U || depth == CV_16U || depth == CV_16S || depth == CV_32F) && cn <= 4 && (m == 3 || m == 5)) )
4326         return false;
4327
4328     Size imgSize = _src.size();
4329     bool useOptimized = (1 == cn) &&
4330                         (size_t)imgSize.width >= localsize[0] * 8  &&
4331                         (size_t)imgSize.height >= localsize[1] * 8 &&
4332                         imgSize.width % 4 == 0 &&
4333                         imgSize.height % 4 == 0 &&
4334                         (ocl::Device::getDefault().isIntel());
4335
4336     cv::String kname = format( useOptimized ? "medianFilter%d_u" : "medianFilter%d", m) ;
4337     cv::String kdefs = useOptimized ?
4338                          format("-D T=%s -D T1=%s -D T4=%s%d -D cn=%d -D USE_4OPT", ocl::typeToStr(type),
4339                          ocl::typeToStr(depth), ocl::typeToStr(depth), cn*4, cn)
4340                          :
4341                          format("-D T=%s -D T1=%s -D cn=%d", ocl::typeToStr(type), ocl::typeToStr(depth), cn) ;
4342
4343     ocl::Kernel k(kname.c_str(), ocl::imgproc::medianFilter_oclsrc, kdefs.c_str() );
4344
4345     if (k.empty())
4346         return false;
4347
4348     UMat src = _src.getUMat();
4349     _dst.create(src.size(), type);
4350     UMat dst = _dst.getUMat();
4351
4352     k.args(ocl::KernelArg::ReadOnlyNoSize(src), ocl::KernelArg::WriteOnly(dst));
4353
4354     if( useOptimized )
4355     {
4356         globalsize[0] = DIVUP(src.cols / 4, localsize[0]) * localsize[0];
4357         globalsize[1] = DIVUP(src.rows / 4, localsize[1]) * localsize[1];
4358     }
4359     else
4360     {
4361         globalsize[0] = (src.cols + localsize[0] + 2) / localsize[0] * localsize[0];
4362         globalsize[1] = (src.rows + localsize[1] - 1) / localsize[1] * localsize[1];
4363     }
4364
4365     return k.run(2, globalsize, localsize, false);
4366 }
4367
4368 #endif
4369
4370 }
4371
4372 #ifdef HAVE_OPENVX
4373 namespace cv
4374 {
4375     namespace ovx {
4376         template <> inline bool skipSmallImages<VX_KERNEL_MEDIAN_3x3>(int w, int h) { return w*h < 1280 * 720; }
4377     }
4378     static bool openvx_medianFilter(InputArray _src, OutputArray _dst, int ksize)
4379     {
4380         if (_src.type() != CV_8UC1 || _dst.type() != CV_8U
4381 #ifndef VX_VERSION_1_1
4382             || ksize != 3
4383 #endif
4384             )
4385             return false;
4386
4387         Mat src = _src.getMat();
4388         Mat dst = _dst.getMat();
4389
4390         if (
4391 #ifdef VX_VERSION_1_1
4392              ksize != 3 ? ovx::skipSmallImages<VX_KERNEL_NON_LINEAR_FILTER>(src.cols, src.rows) :
4393 #endif
4394              ovx::skipSmallImages<VX_KERNEL_MEDIAN_3x3>(src.cols, src.rows)
4395            )
4396             return false;
4397
4398         try
4399         {
4400             ivx::Context ctx = ovx::getOpenVXContext();
4401 #ifdef VX_VERSION_1_1
4402             if ((vx_size)ksize > ctx.nonlinearMaxDimension())
4403                 return false;
4404 #endif
4405
4406             Mat a;
4407             if (dst.data != src.data)
4408                 a = src;
4409             else
4410                 src.copyTo(a);
4411
4412             ivx::Image
4413                 ia = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
4414                     ivx::Image::createAddressing(a.cols, a.rows, 1, (vx_int32)(a.step)), a.data),
4415                 ib = ivx::Image::createFromHandle(ctx, VX_DF_IMAGE_U8,
4416                     ivx::Image::createAddressing(dst.cols, dst.rows, 1, (vx_int32)(dst.step)), dst.data);
4417
4418             //ATTENTION: VX_CONTEXT_IMMEDIATE_BORDER attribute change could lead to strange issues in multi-threaded environments
4419             //since OpenVX standard says nothing about thread-safety for now
4420             ivx::border_t prevBorder = ctx.immediateBorder();
4421             ctx.setImmediateBorder(VX_BORDER_REPLICATE);
4422 #ifdef VX_VERSION_1_1
4423             if (ksize == 3)
4424 #endif
4425             {
4426                 ivx::IVX_CHECK_STATUS(vxuMedian3x3(ctx, ia, ib));
4427             }
4428 #ifdef VX_VERSION_1_1
4429             else
4430             {
4431                 ivx::Matrix mtx;
4432                 if(ksize == 5)
4433                     mtx = ivx::Matrix::createFromPattern(ctx, VX_PATTERN_BOX, ksize, ksize);
4434                 else
4435                 {
4436                     vx_size supportedSize;
4437                     ivx::IVX_CHECK_STATUS(vxQueryContext(ctx, VX_CONTEXT_NONLINEAR_MAX_DIMENSION, &supportedSize, sizeof(supportedSize)));
4438                     if ((vx_size)ksize > supportedSize)
4439                     {
4440                         ctx.setImmediateBorder(prevBorder);
4441                         return false;
4442                     }
4443                     Mat mask(ksize, ksize, CV_8UC1, Scalar(255));
4444                     mtx = ivx::Matrix::create(ctx, VX_TYPE_UINT8, ksize, ksize);
4445                     mtx.copyFrom(mask);
4446                 }
4447                 ivx::IVX_CHECK_STATUS(vxuNonLinearFilter(ctx, VX_NONLINEAR_FILTER_MEDIAN, ia, mtx, ib));
4448             }
4449 #endif
4450             ctx.setImmediateBorder(prevBorder);
4451         }
4452         catch (ivx::RuntimeError & e)
4453         {
4454             VX_DbgThrow(e.what());
4455         }
4456         catch (ivx::WrapperError & e)
4457         {
4458             VX_DbgThrow(e.what());
4459         }
4460
4461         return true;
4462     }
4463 }
4464 #endif
4465
4466 #ifdef HAVE_IPP
4467 namespace cv
4468 {
4469 static bool ipp_medianFilter(Mat &src0, Mat &dst, int ksize)
4470 {
4471     CV_INSTRUMENT_REGION_IPP()
4472
4473 #if IPP_VERSION_X100 < 201801
4474     // Degradations for big kernel
4475     if(ksize > 7)
4476         return false;
4477 #endif
4478
4479     {
4480         int         bufSize;
4481         IppiSize    dstRoiSize = ippiSize(dst.cols, dst.rows), maskSize = ippiSize(ksize, ksize);
4482         IppDataType ippType = ippiGetDataType(src0.type());
4483         int         channels = src0.channels();
4484         IppAutoBuffer<Ipp8u> buffer;
4485
4486         if(src0.isSubmatrix())
4487             return false;
4488
4489         Mat src;
4490         if(dst.data != src0.data)
4491             src = src0;
4492         else
4493             src0.copyTo(src);
4494
4495         if(ippiFilterMedianBorderGetBufferSize(dstRoiSize, maskSize, ippType, channels, &bufSize) < 0)
4496             return false;
4497
4498         buffer.allocate(bufSize);
4499
4500         switch(ippType)
4501         {
4502         case ipp8u:
4503             if(channels == 1)
4504                 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;
4505             else if(channels == 3)
4506                 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;
4507             else if(channels == 4)
4508                 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;
4509             else
4510                 return false;
4511         case ipp16u:
4512             if(channels == 1)
4513                 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;
4514             else if(channels == 3)
4515                 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;
4516             else if(channels == 4)
4517                 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;
4518             else
4519                 return false;
4520         case ipp16s:
4521             if(channels == 1)
4522                 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;
4523             else if(channels == 3)
4524                 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;
4525             else if(channels == 4)
4526                 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;
4527             else
4528                 return false;
4529         case ipp32f:
4530             if(channels == 1)
4531                 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;
4532             else
4533                 return false;
4534         default:
4535             return false;
4536         }
4537     }
4538 }
4539 }
4540 #endif
4541
4542 void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize )
4543 {
4544     CV_INSTRUMENT_REGION()
4545
4546     CV_Assert( (ksize % 2 == 1) && (_src0.dims() <= 2 ));
4547
4548     if( ksize <= 1 || _src0.empty() )
4549     {
4550         _src0.copyTo(_dst);
4551         return;
4552     }
4553
4554     CV_OCL_RUN(_dst.isUMat(),
4555                ocl_medianFilter(_src0,_dst, ksize))
4556
4557     Mat src0 = _src0.getMat();
4558     _dst.create( src0.size(), src0.type() );
4559     Mat dst = _dst.getMat();
4560
4561     CALL_HAL(medianBlur, cv_hal_medianBlur, src0.data, src0.step, dst.data, dst.step, src0.cols, src0.rows, src0.depth(),
4562              src0.channels(), ksize);
4563
4564     CV_OVX_RUN(true,
4565                openvx_medianFilter(_src0, _dst, ksize))
4566
4567     CV_IPP_RUN_FAST(ipp_medianFilter(src0, dst, ksize));
4568
4569 #ifdef HAVE_TEGRA_OPTIMIZATION
4570     if (tegra::useTegra() && tegra::medianBlur(src0, dst, ksize))
4571         return;
4572 #endif
4573
4574     bool useSortNet = ksize == 3 || (ksize == 5
4575 #if !(CV_SIMD128)
4576             && ( src0.depth() > CV_8U || src0.channels() == 2 || src0.channels() > 4 )
4577 #endif
4578         );
4579
4580     Mat src;
4581     if( useSortNet )
4582     {
4583         if( dst.data != src0.data )
4584             src = src0;
4585         else
4586             src0.copyTo(src);
4587
4588         if( src.depth() == CV_8U )
4589             medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize );
4590         else if( src.depth() == CV_16U )
4591             medianBlur_SortNet<MinMax16u, MinMaxVec16u>( src, dst, ksize );
4592         else if( src.depth() == CV_16S )
4593             medianBlur_SortNet<MinMax16s, MinMaxVec16s>( src, dst, ksize );
4594         else if( src.depth() == CV_32F )
4595             medianBlur_SortNet<MinMax32f, MinMaxVec32f>( src, dst, ksize );
4596         else
4597             CV_Error(CV_StsUnsupportedFormat, "");
4598
4599         return;
4600     }
4601     else
4602     {
4603         cv::copyMakeBorder( src0, src, 0, 0, ksize/2, ksize/2, BORDER_REPLICATE|BORDER_ISOLATED);
4604
4605         int cn = src0.channels();
4606         CV_Assert( src.depth() == CV_8U && (cn == 1 || cn == 3 || cn == 4) );
4607
4608         double img_size_mp = (double)(src0.total())/(1 << 20);
4609         if( ksize <= 3 + (img_size_mp < 1 ? 12 : img_size_mp < 4 ? 6 : 2)*
4610             (CV_SIMD128 && hasSIMD128() ? 1 : 3))
4611             medianBlur_8u_Om( src, dst, ksize );
4612         else
4613             medianBlur_8u_O1( src, dst, ksize );
4614     }
4615 }
4616
4617 /****************************************************************************************\
4618                                    Bilateral Filtering
4619 \****************************************************************************************/
4620
4621 namespace cv
4622 {
4623
4624 class BilateralFilter_8u_Invoker :
4625     public ParallelLoopBody
4626 {
4627 public:
4628     BilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, int _radius, int _maxk,
4629         int* _space_ofs, float *_space_weight, float *_color_weight) :
4630         temp(&_temp), dest(&_dest), radius(_radius),
4631         maxk(_maxk), space_ofs(_space_ofs), space_weight(_space_weight), color_weight(_color_weight)
4632     {
4633     }
4634
4635     virtual void operator() (const Range& range) const
4636     {
4637         int i, j, cn = dest->channels(), k;
4638         Size size = dest->size();
4639 #if CV_SIMD128
4640         int CV_DECL_ALIGNED(16) buf[4];
4641         bool haveSIMD128 = hasSIMD128();
4642 #endif
4643
4644         for( i = range.start; i < range.end; i++ )
4645         {
4646             const uchar* sptr = temp->ptr(i+radius) + radius*cn;
4647             uchar* dptr = dest->ptr(i);
4648
4649             if( cn == 1 )
4650             {
4651                 for( j = 0; j < size.width; j++ )
4652                 {
4653                     float sum = 0, wsum = 0;
4654                     int val0 = sptr[j];
4655                     k = 0;
4656 #if CV_SIMD128
4657                     if( haveSIMD128 )
4658                     {
4659                         v_float32x4 _val0 = v_setall_f32(static_cast<float>(val0));
4660                         v_float32x4 vsumw = v_setzero_f32();
4661                         v_float32x4 vsumc = v_setzero_f32();
4662
4663                         for( ; k <= maxk - 4; k += 4 )
4664                         {
4665                             v_float32x4 _valF = v_float32x4(sptr[j + space_ofs[k]],
4666                                 sptr[j + space_ofs[k + 1]],
4667                                 sptr[j + space_ofs[k + 2]],
4668                                 sptr[j + space_ofs[k + 3]]);
4669                             v_float32x4 _val = v_abs(_valF - _val0);
4670                             v_store(buf, v_round(_val));
4671
4672                             v_float32x4 _cw = v_float32x4(color_weight[buf[0]],
4673                                 color_weight[buf[1]],
4674                                 color_weight[buf[2]],
4675                                 color_weight[buf[3]]);
4676                             v_float32x4 _sw = v_load(space_weight+k);
4677 #if defined(_MSC_VER) && _MSC_VER == 1700/* MSVS 2012 */ && CV_AVX
4678                             // details: https://github.com/opencv/opencv/issues/11004
4679                             vsumw += _cw * _sw;
4680                             vsumc += _cw * _sw * _valF;
4681 #else
4682                             v_float32x4 _w = _cw * _sw;
4683                             _cw = _w * _valF;
4684
4685                             vsumw += _w;
4686                             vsumc += _cw;
4687 #endif
4688                         }
4689                         float *bufFloat = (float*)buf;
4690                         v_float32x4 sum4 = v_reduce_sum4(vsumw, vsumc, vsumw, vsumc);
4691                         v_store(bufFloat, sum4);
4692                         sum += bufFloat[1];
4693                         wsum += bufFloat[0];
4694                     }
4695 #endif
4696                     for( ; k < maxk; k++ )
4697                     {
4698                         int val = sptr[j + space_ofs[k]];
4699                         float w = space_weight[k]*color_weight[std::abs(val - val0)];
4700                         sum += val*w;
4701                         wsum += w;
4702                     }
4703                     // overflow is not possible here => there is no need to use cv::saturate_cast
4704                     dptr[j] = (uchar)cvRound(sum/wsum);
4705                 }
4706             }
4707             else
4708             {
4709                 assert( cn == 3 );
4710                 for( j = 0; j < size.width*3; j += 3 )
4711                 {
4712                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
4713                     int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
4714                     k = 0;
4715 #if CV_SIMD128
4716                     if( haveSIMD128 )
4717                     {
4718                         v_float32x4 vsumw = v_setzero_f32();
4719                         v_float32x4 vsumb = v_setzero_f32();
4720                         v_float32x4 vsumg = v_setzero_f32();
4721                         v_float32x4 vsumr = v_setzero_f32();
4722                         const v_float32x4 _b0 = v_setall_f32(static_cast<float>(b0));
4723                         const v_float32x4 _g0 = v_setall_f32(static_cast<float>(g0));
4724                         const v_float32x4 _r0 = v_setall_f32(static_cast<float>(r0));
4725
4726                         for( ; k <= maxk - 4; k += 4 )
4727                         {
4728                             const uchar* const sptr_k0  = sptr + j + space_ofs[k];
4729                             const uchar* const sptr_k1  = sptr + j + space_ofs[k+1];
4730                             const uchar* const sptr_k2  = sptr + j + space_ofs[k+2];
4731                             const uchar* const sptr_k3  = sptr + j + space_ofs[k+3];
4732
4733                             v_float32x4 __b = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k0)));
4734                             v_float32x4 __g = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k1)));
4735                             v_float32x4 __r = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k2)));
4736                             v_float32x4 __z = v_cvt_f32(v_reinterpret_as_s32(v_load_expand_q(sptr_k3)));
4737                             v_float32x4 _b, _g, _r, _z;
4738
4739                             v_transpose4x4(__b, __g, __r, __z, _b, _g, _r, _z);
4740
4741                             v_float32x4 bt = v_abs(_b -_b0);
4742                             v_float32x4 gt = v_abs(_g -_g0);
4743                             v_float32x4 rt = v_abs(_r -_r0);
4744
4745                             bt = rt + bt + gt;
4746                             v_store(buf, v_round(bt));
4747
4748                             v_float32x4 _w  = v_float32x4(color_weight[buf[0]],color_weight[buf[1]],
4749                                                     color_weight[buf[2]],color_weight[buf[3]]);
4750                             v_float32x4 _sw = v_load(space_weight+k);
4751
4752 #if defined(_MSC_VER) && _MSC_VER == 1700/* MSVS 2012 */ && CV_AVX
4753                             // details: https://github.com/opencv/opencv/issues/11004
4754                             vsumw += _w * _sw;
4755                             vsumb += _w * _sw * _b;
4756                             vsumg += _w * _sw * _g;
4757                             vsumr += _w * _sw * _r;
4758 #else
4759                             _w *= _sw;
4760                             _b *=  _w;
4761                             _g *=  _w;
4762                             _r *=  _w;
4763
4764                             vsumw += _w;
4765                             vsumb += _b;
4766                             vsumg += _g;
4767                             vsumr += _r;
4768 #endif
4769                         }
4770                         float *bufFloat = (float*)buf;
4771                         v_float32x4 sum4 = v_reduce_sum4(vsumw, vsumb, vsumg, vsumr);
4772                         v_store(bufFloat, sum4);
4773                         wsum += bufFloat[0];
4774                         sum_b += bufFloat[1];
4775                         sum_g += bufFloat[2];
4776                         sum_r += bufFloat[3];
4777                     }
4778 #endif
4779
4780                     for( ; k < maxk; k++ )
4781                     {
4782                         const uchar* sptr_k = sptr + j + space_ofs[k];
4783                         int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
4784                         float w = space_weight[k]*color_weight[std::abs(b - b0) +
4785                                                                std::abs(g - g0) + std::abs(r - r0)];
4786                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
4787                         wsum += w;
4788                     }
4789                     wsum = 1.f/wsum;
4790                     b0 = cvRound(sum_b*wsum);
4791                     g0 = cvRound(sum_g*wsum);
4792                     r0 = cvRound(sum_r*wsum);
4793                     dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
4794                 }
4795             }
4796         }
4797     }
4798
4799 private:
4800     const Mat *temp;
4801     Mat *dest;
4802     int radius, maxk, *space_ofs;
4803     float *space_weight, *color_weight;
4804 };
4805
4806 #ifdef HAVE_OPENCL
4807
4808 static bool ocl_bilateralFilter_8u(InputArray _src, OutputArray _dst, int d,
4809                                    double sigma_color, double sigma_space,
4810                                    int borderType)
4811 {
4812 #ifdef __ANDROID__
4813     if (ocl::Device::getDefault().isNVidia())
4814         return false;
4815 #endif
4816
4817     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
4818     int i, j, maxk, radius;
4819
4820     if (depth != CV_8U || cn > 4)
4821         return false;
4822
4823     if (sigma_color <= 0)
4824         sigma_color = 1;
4825     if (sigma_space <= 0)
4826         sigma_space = 1;
4827
4828     double gauss_color_coeff = -0.5 / (sigma_color * sigma_color);
4829     double gauss_space_coeff = -0.5 / (sigma_space * sigma_space);
4830
4831     if ( d <= 0 )
4832         radius = cvRound(sigma_space * 1.5);
4833     else
4834         radius = d / 2;
4835     radius = MAX(radius, 1);
4836     d = radius * 2 + 1;
4837
4838     UMat src = _src.getUMat(), dst = _dst.getUMat(), temp;
4839     if (src.u == dst.u)
4840         return false;
4841
4842     copyMakeBorder(src, temp, radius, radius, radius, radius, borderType);
4843     std::vector<float> _space_weight(d * d);
4844     std::vector<int> _space_ofs(d * d);
4845     float * const space_weight = &_space_weight[0];
4846     int * const space_ofs = &_space_ofs[0];
4847
4848     // initialize space-related bilateral filter coefficients
4849     for( i = -radius, maxk = 0; i <= radius; i++ )
4850         for( j = -radius; j <= radius; j++ )
4851         {
4852             double r = std::sqrt((double)i * i + (double)j * j);
4853             if ( r > radius )
4854                 continue;
4855             space_weight[maxk] = (float)std::exp(r * r * gauss_space_coeff);
4856             space_ofs[maxk++] = (int)(i * temp.step + j * cn);
4857         }
4858
4859     char cvt[3][40];
4860     String cnstr = cn > 1 ? format("%d", cn) : "";
4861     String kernelName("bilateral");
4862     size_t sizeDiv = 1;
4863     if ((ocl::Device::getDefault().isIntel()) &&
4864         (ocl::Device::getDefault().type() == ocl::Device::TYPE_GPU))
4865     {
4866             //Intel GPU
4867             if (dst.cols % 4 == 0 && cn == 1) // For single channel x4 sized images.
4868             {
4869                 kernelName = "bilateral_float4";
4870                 sizeDiv = 4;
4871             }
4872      }
4873      ocl::Kernel k(kernelName.c_str(), ocl::imgproc::bilateral_oclsrc,
4874             format("-D radius=%d -D maxk=%d -D cn=%d -D int_t=%s -D uint_t=uint%s -D convert_int_t=%s"
4875             " -D uchar_t=%s -D float_t=%s -D convert_float_t=%s -D convert_uchar_t=%s -D gauss_color_coeff=(float)%f",
4876             radius, maxk, cn, ocl::typeToStr(CV_32SC(cn)), cnstr.c_str(),
4877             ocl::convertTypeStr(CV_8U, CV_32S, cn, cvt[0]),
4878             ocl::typeToStr(type), ocl::typeToStr(CV_32FC(cn)),
4879             ocl::convertTypeStr(CV_32S, CV_32F, cn, cvt[1]),
4880             ocl::convertTypeStr(CV_32F, CV_8U, cn, cvt[2]), gauss_color_coeff));
4881     if (k.empty())
4882         return false;
4883
4884     Mat mspace_weight(1, d * d, CV_32FC1, space_weight);
4885     Mat mspace_ofs(1, d * d, CV_32SC1, space_ofs);
4886     UMat ucolor_weight, uspace_weight, uspace_ofs;
4887
4888     mspace_weight.copyTo(uspace_weight);
4889     mspace_ofs.copyTo(uspace_ofs);
4890
4891     k.args(ocl::KernelArg::ReadOnlyNoSize(temp), ocl::KernelArg::WriteOnly(dst),
4892            ocl::KernelArg::PtrReadOnly(uspace_weight),
4893            ocl::KernelArg::PtrReadOnly(uspace_ofs));
4894
4895     size_t globalsize[2] = { (size_t)dst.cols / sizeDiv, (size_t)dst.rows };
4896     return k.run(2, globalsize, NULL, false);
4897 }
4898
4899 #endif
4900 static void
4901 bilateralFilter_8u( const Mat& src, Mat& dst, int d,
4902     double sigma_color, double sigma_space,
4903     int borderType )
4904 {
4905     int cn = src.channels();
4906     int i, j, maxk, radius;
4907     Size size = src.size();
4908
4909     CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) && src.data != dst.data );
4910
4911     if( sigma_color <= 0 )
4912         sigma_color = 1;
4913     if( sigma_space <= 0 )
4914         sigma_space = 1;
4915
4916     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
4917     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
4918
4919     if( d <= 0 )
4920         radius = cvRound(sigma_space*1.5);
4921     else
4922         radius = d/2;
4923     radius = MAX(radius, 1);
4924     d = radius*2 + 1;
4925
4926     Mat temp;
4927     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
4928
4929     std::vector<float> _color_weight(cn*256);
4930     std::vector<float> _space_weight(d*d);
4931     std::vector<int> _space_ofs(d*d);
4932     float* color_weight = &_color_weight[0];
4933     float* space_weight = &_space_weight[0];
4934     int* space_ofs = &_space_ofs[0];
4935
4936     // initialize color-related bilateral filter coefficients
4937
4938     for( i = 0; i < 256*cn; i++ )
4939         color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);
4940
4941     // initialize space-related bilateral filter coefficients
4942     for( i = -radius, maxk = 0; i <= radius; i++ )
4943     {
4944         j = -radius;
4945
4946         for( ; j <= radius; j++ )
4947         {
4948             double r = std::sqrt((double)i*i + (double)j*j);
4949             if( r > radius )
4950                 continue;
4951             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
4952             space_ofs[maxk++] = (int)(i*temp.step + j*cn);
4953         }
4954     }
4955
4956     BilateralFilter_8u_Invoker body(dst, temp, radius, maxk, space_ofs, space_weight, color_weight);
4957     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
4958 }
4959
4960
4961 class BilateralFilter_32f_Invoker :
4962     public ParallelLoopBody
4963 {
4964 public:
4965
4966     BilateralFilter_32f_Invoker(int _cn, int _radius, int _maxk, int *_space_ofs,
4967         const Mat& _temp, Mat& _dest, float _scale_index, float *_space_weight, float *_expLUT) :
4968         cn(_cn), radius(_radius), maxk(_maxk), space_ofs(_space_ofs),
4969         temp(&_temp), dest(&_dest), scale_index(_scale_index), space_weight(_space_weight), expLUT(_expLUT)
4970     {
4971     }
4972
4973     virtual void operator() (const Range& range) const
4974     {
4975         int i, j, k;
4976         Size size = dest->size();
4977 #if CV_SIMD128
4978         int CV_DECL_ALIGNED(16) idxBuf[4];
4979         bool haveSIMD128 = hasSIMD128();
4980 #endif
4981
4982         for( i = range.start; i < range.end; i++ )
4983         {
4984             const float* sptr = temp->ptr<float>(i+radius) + radius*cn;
4985             float* dptr = dest->ptr<float>(i);
4986
4987             if( cn == 1 )
4988             {
4989                 for( j = 0; j < size.width; j++ )
4990                 {
4991                     float sum = 0, wsum = 0;
4992                     float val0 = sptr[j];
4993                     k = 0;
4994 #if CV_SIMD128
4995                     if( haveSIMD128 )
4996                     {
4997                         v_float32x4 vecwsum = v_setzero_f32();
4998                         v_float32x4 vecvsum = v_setzero_f32();
4999                         const v_float32x4 _val0 = v_setall_f32(sptr[j]);
5000                         const v_float32x4 _scale_index = v_setall_f32(scale_index);
5001
5002                         for (; k <= maxk - 4; k += 4)
5003                         {
5004                             v_float32x4 _sw = v_load(space_weight + k);
5005                             v_float32x4 _val = v_float32x4(sptr[j + space_ofs[k]],
5006                                 sptr[j + space_ofs[k + 1]],
5007                                 sptr[j + space_ofs[k + 2]],
5008                                 sptr[j + space_ofs[k + 3]]);
5009                             v_float32x4 _alpha = v_abs(_val - _val0) * _scale_index;
5010
5011                             v_int32x4 _idx = v_round(_alpha);
5012                             v_store(idxBuf, _idx);
5013                             _alpha -= v_cvt_f32(_idx);
5014
5015                             v_float32x4 _explut = v_float32x4(expLUT[idxBuf[0]],
5016                                 expLUT[idxBuf[1]],
5017                                 expLUT[idxBuf[2]],
5018                                 expLUT[idxBuf[3]]);
5019                             v_float32x4 _explut1 = v_float32x4(expLUT[idxBuf[0] + 1],
5020                                 expLUT[idxBuf[1] + 1],
5021                                 expLUT[idxBuf[2] + 1],
5022                                 expLUT[idxBuf[3] + 1]);
5023
5024                             v_float32x4 _w = _sw * (_explut + (_alpha * (_explut1 - _explut)));
5025                             _val *= _w;
5026
5027                             vecwsum += _w;
5028                             vecvsum += _val;
5029                         }
5030                         float *bufFloat = (float*)idxBuf;
5031                         v_float32x4 sum4 = v_reduce_sum4(vecwsum, vecvsum, vecwsum, vecvsum);
5032                         v_store(bufFloat, sum4);
5033                         sum += bufFloat[1];
5034                         wsum += bufFloat[0];
5035                     }
5036 #endif
5037
5038                     for( ; k < maxk; k++ )
5039                     {
5040                         float val = sptr[j + space_ofs[k]];
5041                         float alpha = (float)(std::abs(val - val0)*scale_index);
5042                         int idx = cvFloor(alpha);
5043                         alpha -= idx;
5044                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
5045                         sum += val*w;
5046                         wsum += w;
5047                     }
5048                     dptr[j] = (float)(sum/wsum);
5049                 }
5050             }
5051             else
5052             {
5053                 CV_Assert( cn == 3 );
5054                 for( j = 0; j < size.width*3; j += 3 )
5055                 {
5056                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
5057                     float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
5058                     k = 0;
5059 #if CV_SIMD128
5060                     if( haveSIMD128 )
5061                     {
5062                         v_float32x4 sumw = v_setzero_f32();
5063                         v_float32x4 sumb = v_setzero_f32();
5064                         v_float32x4 sumg = v_setzero_f32();
5065                         v_float32x4 sumr = v_setzero_f32();
5066                         const v_float32x4 _b0 = v_setall_f32(b0);
5067                         const v_float32x4 _g0 = v_setall_f32(g0);
5068                         const v_float32x4 _r0 = v_setall_f32(r0);
5069                         const v_float32x4 _scale_index = v_setall_f32(scale_index);
5070
5071                         for( ; k <= maxk-4; k += 4 )
5072                         {
5073                             v_float32x4 _sw = v_load(space_weight + k);
5074
5075                             const float* const sptr_k0 = sptr + j + space_ofs[k];
5076                             const float* const sptr_k1 = sptr + j + space_ofs[k+1];
5077                             const float* const sptr_k2 = sptr + j + space_ofs[k+2];
5078                             const float* const sptr_k3 = sptr + j + space_ofs[k+3];
5079
5080                             v_float32x4 _v0 = v_load(sptr_k0);
5081                             v_float32x4 _v1 = v_load(sptr_k1);
5082                             v_float32x4 _v2 = v_load(sptr_k2);
5083                             v_float32x4 _v3 = v_load(sptr_k3);
5084                             v_float32x4 _b, _g, _r, _dummy;
5085
5086                             v_transpose4x4(_v0, _v1, _v2, _v3, _b, _g, _r, _dummy);
5087
5088                             v_float32x4 _bt = v_abs(_b - _b0);
5089                             v_float32x4 _gt = v_abs(_g - _g0);
5090                             v_float32x4 _rt = v_abs(_r - _r0);
5091                             v_float32x4 _alpha = _scale_index * (_bt + _gt + _rt);
5092
5093                             v_int32x4 _idx = v_round(_alpha);
5094                             v_store((int*)idxBuf, _idx);
5095                             _alpha -= v_cvt_f32(_idx);
5096
5097                             v_float32x4 _explut = v_float32x4(expLUT[idxBuf[0]],
5098                                 expLUT[idxBuf[1]],
5099                                 expLUT[idxBuf[2]],
5100                                 expLUT[idxBuf[3]]);
5101                             v_float32x4 _explut1 = v_float32x4(expLUT[idxBuf[0] + 1],
5102                                 expLUT[idxBuf[1] + 1],
5103                                 expLUT[idxBuf[2] + 1],
5104                                 expLUT[idxBuf[3] + 1]);
5105
5106                             v_float32x4 _w = _sw * (_explut + (_alpha * (_explut1 - _explut)));
5107
5108                             _b *=  _w;
5109                             _g *=  _w;
5110                             _r *=  _w;
5111                             sumw += _w;
5112                             sumb += _b;
5113                             sumg += _g;
5114                             sumr += _r;
5115                         }
5116                         v_float32x4 sum4 = v_reduce_sum4(sumw, sumb, sumg, sumr);
5117                         float *bufFloat = (float*)idxBuf;
5118                         v_store(bufFloat, sum4);
5119                         wsum += bufFloat[0];
5120                         sum_b += bufFloat[1];
5121                         sum_g += bufFloat[2];
5122                         sum_r += bufFloat[3];
5123                     }
5124 #endif
5125
5126                     for(; k < maxk; k++ )
5127                     {
5128                         const float* sptr_k = sptr + j + space_ofs[k];
5129                         float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
5130                         float alpha = (float)((std::abs(b - b0) +
5131                             std::abs(g - g0) + std::abs(r - r0))*scale_index);
5132                         int idx = cvFloor(alpha);
5133                         alpha -= idx;
5134                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
5135                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
5136                         wsum += w;
5137                     }
5138                     wsum = 1.f/wsum;
5139                     b0 = sum_b*wsum;
5140                     g0 = sum_g*wsum;
5141                     r0 = sum_r*wsum;
5142                     dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
5143                 }
5144             }
5145         }
5146     }
5147
5148 private:
5149     int cn, radius, maxk, *space_ofs;
5150     const Mat* temp;
5151     Mat *dest;
5152     float scale_index, *space_weight, *expLUT;
5153 };
5154
5155
5156 static void
5157 bilateralFilter_32f( const Mat& src, Mat& dst, int d,
5158                      double sigma_color, double sigma_space,
5159                      int borderType )
5160 {
5161     int cn = src.channels();
5162     int i, j, maxk, radius;
5163     double minValSrc=-1, maxValSrc=1;
5164     const int kExpNumBinsPerChannel = 1 << 12;
5165     int kExpNumBins = 0;
5166     float lastExpVal = 1.f;
5167     float len, scale_index;
5168     Size size = src.size();
5169
5170     CV_Assert( (src.type() == CV_32FC1 || src.type() == CV_32FC3) && src.data != dst.data );
5171
5172     if( sigma_color <= 0 )
5173         sigma_color = 1;
5174     if( sigma_space <= 0 )
5175         sigma_space = 1;
5176
5177     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
5178     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
5179
5180     if( d <= 0 )
5181         radius = cvRound(sigma_space*1.5);
5182     else
5183         radius = d/2;
5184     radius = MAX(radius, 1);
5185     d = radius*2 + 1;
5186     // compute the min/max range for the input image (even if multichannel)
5187
5188     minMaxLoc( src.reshape(1), &minValSrc, &maxValSrc );
5189     if(std::abs(minValSrc - maxValSrc) < FLT_EPSILON)
5190     {
5191         src.copyTo(dst);
5192         return;
5193     }
5194
5195     // temporary copy of the image with borders for easy processing
5196     Mat temp;
5197     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
5198     const double insteadNaNValue = -5. * sigma_color;
5199     patchNaNs( temp, insteadNaNValue ); // this replacement of NaNs makes the assumption that depth values are nonnegative
5200                                         // TODO: make insteadNaNValue avalible in the outside function interface to control the cases breaking the assumption
5201     // allocate lookup tables
5202     std::vector<float> _space_weight(d*d);
5203     std::vector<int> _space_ofs(d*d);
5204     float* space_weight = &_space_weight[0];
5205     int* space_ofs = &_space_ofs[0];
5206
5207     // assign a length which is slightly more than needed
5208     len = (float)(maxValSrc - minValSrc) * cn;
5209     kExpNumBins = kExpNumBinsPerChannel * cn;
5210     std::vector<float> _expLUT(kExpNumBins+2);
5211     float* expLUT = &_expLUT[0];
5212
5213     scale_index = kExpNumBins/len;
5214
5215     // initialize the exp LUT
5216     for( i = 0; i < kExpNumBins+2; i++ )
5217     {
5218         if( lastExpVal > 0.f )
5219         {
5220             double val =  i / scale_index;
5221             expLUT[i] = (float)std::exp(val * val * gauss_color_coeff);
5222             lastExpVal = expLUT[i];
5223         }
5224         else
5225             expLUT[i] = 0.f;
5226     }
5227
5228     // initialize space-related bilateral filter coefficients
5229     for( i = -radius, maxk = 0; i <= radius; i++ )
5230         for( j = -radius; j <= radius; j++ )
5231         {
5232             double r = std::sqrt((double)i*i + (double)j*j);
5233             if( r > radius )
5234                 continue;
5235             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
5236             space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
5237         }
5238
5239     // parallel_for usage
5240
5241     BilateralFilter_32f_Invoker body(cn, radius, maxk, space_ofs, temp, dst, scale_index, space_weight, expLUT);
5242     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
5243 }
5244
5245 #ifdef HAVE_IPP
5246 #define IPP_BILATERAL_PARALLEL 1
5247
5248 #ifdef HAVE_IPP_IW
5249 class ipp_bilateralFilterParallel: public ParallelLoopBody
5250 {
5251 public:
5252     ipp_bilateralFilterParallel(::ipp::IwiImage &_src, ::ipp::IwiImage &_dst, int _radius, Ipp32f _valSquareSigma, Ipp32f _posSquareSigma, ::ipp::IwiBorderType _borderType, bool *_ok):
5253         src(_src), dst(_dst)
5254     {
5255         pOk = _ok;
5256
5257         radius          = _radius;
5258         valSquareSigma  = _valSquareSigma;
5259         posSquareSigma  = _posSquareSigma;
5260         borderType      = _borderType;
5261
5262         *pOk = true;
5263     }
5264     ~ipp_bilateralFilterParallel() {}
5265
5266     virtual void operator() (const Range& range) const
5267     {
5268         if(*pOk == false)
5269             return;
5270
5271         try
5272         {
5273             ::ipp::IwiTile tile = ::ipp::IwiRoi(0, range.start, dst.m_size.width, range.end - range.start);
5274             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBilateral, src, dst, radius, valSquareSigma, posSquareSigma, ::ipp::IwDefault(), borderType, tile);
5275         }
5276         catch(::ipp::IwException)
5277         {
5278             *pOk = false;
5279             return;
5280         }
5281     }
5282 private:
5283     ::ipp::IwiImage &src;
5284     ::ipp::IwiImage &dst;
5285
5286     int                  radius;
5287     Ipp32f               valSquareSigma;
5288     Ipp32f               posSquareSigma;
5289     ::ipp::IwiBorderType borderType;
5290
5291     bool  *pOk;
5292     const ipp_bilateralFilterParallel& operator= (const ipp_bilateralFilterParallel&);
5293 };
5294 #endif
5295
5296 static bool ipp_bilateralFilter(Mat &src, Mat &dst, int d, double sigmaColor, double sigmaSpace, int borderType)
5297 {
5298 #ifdef HAVE_IPP_IW
5299     CV_INSTRUMENT_REGION_IPP()
5300
5301     int         radius         = IPP_MAX(((d <= 0)?cvRound(sigmaSpace*1.5):d/2), 1);
5302     Ipp32f      valSquareSigma = (Ipp32f)((sigmaColor <= 0)?1:sigmaColor*sigmaColor);
5303     Ipp32f      posSquareSigma = (Ipp32f)((sigmaSpace <= 0)?1:sigmaSpace*sigmaSpace);
5304
5305     // Acquire data and begin processing
5306     try
5307     {
5308         ::ipp::IwiImage      iwSrc = ippiGetImage(src);
5309         ::ipp::IwiImage      iwDst = ippiGetImage(dst);
5310         ::ipp::IwiBorderSize borderSize(radius);
5311         ::ipp::IwiBorderType ippBorder(ippiGetBorder(iwSrc, borderType, borderSize));
5312         if(!ippBorder)
5313             return false;
5314
5315         const int threads = ippiSuggestThreadsNum(iwDst, 2);
5316         if(IPP_BILATERAL_PARALLEL && threads > 1) {
5317             bool  ok      = true;
5318             Range range(0, (int)iwDst.m_size.height);
5319             ipp_bilateralFilterParallel invoker(iwSrc, iwDst, radius, valSquareSigma, posSquareSigma, ippBorder, &ok);
5320             if(!ok)
5321                 return false;
5322
5323             parallel_for_(range, invoker, threads*4);
5324
5325             if(!ok)
5326                 return false;
5327         } else {
5328             CV_INSTRUMENT_FUN_IPP(::ipp::iwiFilterBilateral, iwSrc, iwDst, radius, valSquareSigma, posSquareSigma, ::ipp::IwDefault(), ippBorder);
5329         }
5330     }
5331     catch (::ipp::IwException)
5332     {
5333         return false;
5334     }
5335     return true;
5336 #else
5337     CV_UNUSED(src); CV_UNUSED(dst); CV_UNUSED(d); CV_UNUSED(sigmaColor); CV_UNUSED(sigmaSpace); CV_UNUSED(borderType);
5338     return false;
5339 #endif
5340 }
5341 #endif
5342
5343 }
5344
5345 void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,
5346                       double sigmaColor, double sigmaSpace,
5347                       int borderType )
5348 {
5349     CV_INSTRUMENT_REGION()
5350
5351     _dst.create( _src.size(), _src.type() );
5352
5353     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
5354                ocl_bilateralFilter_8u(_src, _dst, d, sigmaColor, sigmaSpace, borderType))
5355
5356     Mat src = _src.getMat(), dst = _dst.getMat();
5357
5358     CV_IPP_RUN_FAST(ipp_bilateralFilter(src, dst, d, sigmaColor, sigmaSpace, borderType));
5359
5360     if( src.depth() == CV_8U )
5361         bilateralFilter_8u( src, dst, d, sigmaColor, sigmaSpace, borderType );
5362     else if( src.depth() == CV_32F )
5363         bilateralFilter_32f( src, dst, d, sigmaColor, sigmaSpace, borderType );
5364     else
5365         CV_Error( CV_StsUnsupportedFormat,
5366         "Bilateral filtering is only implemented for 8u and 32f images" );
5367 }
5368
5369 //////////////////////////////////////////////////////////////////////////////////////////
5370
5371 CV_IMPL void
5372 cvSmooth( const void* srcarr, void* dstarr, int smooth_type,
5373           int param1, int param2, double param3, double param4 )
5374 {
5375     cv::Mat src = cv::cvarrToMat(srcarr), dst0 = cv::cvarrToMat(dstarr), dst = dst0;
5376
5377     CV_Assert( dst.size() == src.size() &&
5378         (smooth_type == CV_BLUR_NO_SCALE || dst.type() == src.type()) );
5379
5380     if( param2 <= 0 )
5381         param2 = param1;
5382
5383     if( smooth_type == CV_BLUR || smooth_type == CV_BLUR_NO_SCALE )
5384         cv::boxFilter( src, dst, dst.depth(), cv::Size(param1, param2), cv::Point(-1,-1),
5385             smooth_type == CV_BLUR, cv::BORDER_REPLICATE );
5386     else if( smooth_type == CV_GAUSSIAN )
5387         cv::GaussianBlur( src, dst, cv::Size(param1, param2), param3, param4, cv::BORDER_REPLICATE );
5388     else if( smooth_type == CV_MEDIAN )
5389         cv::medianBlur( src, dst, param1 );
5390     else
5391         cv::bilateralFilter( src, dst, param1, param3, param4, cv::BORDER_REPLICATE );
5392
5393     if( dst.data != dst0.data )
5394         CV_Error( CV_StsUnmatchedFormats, "The destination image does not have the proper type" );
5395 }
5396
5397 /* End of file. */