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 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "precomp.hpp"
  44 #include "opencl_kernels_imgproc.hpp"
  45
  46 /*
  47  * This file includes the code, contributed by Simon Perreault
  48  * (the function icvMedianBlur_8u_O1)
  49  *
  50  * Constant-time median filtering -- http://nomis80.org/ctmf.html
  51  * Copyright (C) 2006 Simon Perreault
  52  *
  53  * Contact:
  54  *  Laboratoire de vision et systemes numeriques
  55  *  Pavillon Adrien-Pouliot
  56  *  Universite Laval
  57  *  Sainte-Foy, Quebec, Canada
  58  *  G1K 7P4
  59  *
  60  *  perreaul@gel.ulaval.ca
  61  */
  62
  63 namespace cv
  64 {
  65
  66 /****************************************************************************************\
  67                                          Box Filter
  68 \****************************************************************************************/
  69
  70 template<typename T, typename ST>
  71 struct RowSum :
  72         public BaseRowFilter
  73 {
  74     RowSum( int _ksize, int _anchor ) :
  75         BaseRowFilter()
  76     {
  77         ksize = _ksize;
  78         anchor = _anchor;
  79     }
  80
  81     virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
  82     {
  83         const T* S = (const T*)src;
  84         ST* D = (ST*)dst;
  85         int i = 0, k, ksz_cn = ksize*cn;
  86
  87         width = (width - 1)*cn;
  88         for( k = 0; k < cn; k++, S++, D++ )
  89         {
  90             ST s = 0;
  91             for( i = 0; i < ksz_cn; i += cn )
  92                 s += S[i];
  93             D[0] = s;
  94             for( i = 0; i < width; i += cn )
  95             {
  96                 s += S[i + ksz_cn] - S[i];
  97                 D[i+cn] = s;
  98             }
  99         }
 100     }
 101 };
 102
 103
 104 template<typename ST, typename T>
 105 struct ColumnSum :
 106         public BaseColumnFilter
 107 {
 108     ColumnSum( int _ksize, int _anchor, double _scale ) :
 109         BaseColumnFilter()
 110     {
 111         ksize = _ksize;
 112         anchor = _anchor;
 113         scale = _scale;
 114         sumCount = 0;
 115     }
 116
 117     virtual void reset() { sumCount = 0; }
 118
 119     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 120     {
 121         int i;
 122         ST* SUM;
 123         bool haveScale = scale != 1;
 124         double _scale = scale;
 125
 126         if( width != (int)sum.size() )
 127         {
 128             sum.resize(width);
 129             sumCount = 0;
 130         }
 131
 132         SUM = &sum[0];
 133         if( sumCount == 0 )
 134         {
 135             memset((void*)SUM, 0, width*sizeof(ST));
 136
 137             for( ; sumCount < ksize - 1; sumCount++, src++ )
 138             {
 139                 const ST* Sp = (const ST*)src[0];
 140                 for( i = 0; i <= width - 2; i += 2 )
 141                 {
 142                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 143                     SUM[i] = s0; SUM[i+1] = s1;
 144                 }
 145
 146                 for( ; i < width; i++ )
 147                     SUM[i] += Sp[i];
 148             }
 149         }
 150         else
 151         {
 152             CV_Assert( sumCount == ksize-1 );
 153             src += ksize-1;
 154         }
 155
 156         for( ; count--; src++ )
 157         {
 158             const ST* Sp = (const ST*)src[0];
 159             const ST* Sm = (const ST*)src[1-ksize];
 160             T* D = (T*)dst;
 161             if( haveScale )
 162             {
 163                 for( i = 0; i <= width - 2; i += 2 )
 164                 {
 165                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 166                     D[i] = saturate_cast<T>(s0*_scale);
 167                     D[i+1] = saturate_cast<T>(s1*_scale);
 168                     s0 -= Sm[i]; s1 -= Sm[i+1];
 169                     SUM[i] = s0; SUM[i+1] = s1;
 170                 }
 171
 172                 for( ; i < width; i++ )
 173                 {
 174                     ST s0 = SUM[i] + Sp[i];
 175                     D[i] = saturate_cast<T>(s0*_scale);
 176                     SUM[i] = s0 - Sm[i];
 177                 }
 178             }
 179             else
 180             {
 181                 for( i = 0; i <= width - 2; i += 2 )
 182                 {
 183                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 184                     D[i] = saturate_cast<T>(s0);
 185                     D[i+1] = saturate_cast<T>(s1);
 186                     s0 -= Sm[i]; s1 -= Sm[i+1];
 187                     SUM[i] = s0; SUM[i+1] = s1;
 188                 }
 189
 190                 for( ; i < width; i++ )
 191                 {
 192                     ST s0 = SUM[i] + Sp[i];
 193                     D[i] = saturate_cast<T>(s0);
 194                     SUM[i] = s0 - Sm[i];
 195                 }
 196             }
 197             dst += dststep;
 198         }
 199     }
 200
 201     double scale;
 202     int sumCount;
 203     std::vector<ST> sum;
 204 };
 205
 206
 207 template<>
 208 struct ColumnSum<int, uchar> :
 209         public BaseColumnFilter
 210 {
 211     ColumnSum( int _ksize, int _anchor, double _scale ) :
 212         BaseColumnFilter()
 213     {
 214         ksize = _ksize;
 215         anchor = _anchor;
 216         scale = _scale;
 217         sumCount = 0;
 218     }
 219
 220     virtual void reset() { sumCount = 0; }
 221
 222     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 223     {
 224         int i;
 225         int* SUM;
 226         bool haveScale = scale != 1;
 227         double _scale = scale;
 228
 229         #if CV_SSE2
 230             bool haveSSE2 = checkHardwareSupport(CV_CPU_SSE2);
 231         #endif
 232
 233         if( width != (int)sum.size() )
 234         {
 235             sum.resize(width);
 236             sumCount = 0;
 237         }
 238
 239         SUM = &sum[0];
 240         if( sumCount == 0 )
 241         {
 242             memset((void*)SUM, 0, width*sizeof(int));
 243             for( ; sumCount < ksize - 1; sumCount++, src++ )
 244             {
 245                 const int* Sp = (const int*)src[0];
 246                 i = 0;
 247                 #if CV_SSE2
 248                 if(haveSSE2)
 249                 {
 250                     for( ; i <= width-4; i+=4 )
 251                     {
 252                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 253                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 254                         _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
 255                     }
 256                 }
 257                 #elif CV_NEON
 258                 for( ; i <= width - 4; i+=4 )
 259                     vst1q_s32(SUM + i, vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i)));
 260                 #endif
 261                 for( ; i < width; i++ )
 262                     SUM[i] += Sp[i];
 263             }
 264         }
 265         else
 266         {
 267             CV_Assert( sumCount == ksize-1 );
 268             src += ksize-1;
 269         }
 270
 271         for( ; count--; src++ )
 272         {
 273             const int* Sp = (const int*)src[0];
 274             const int* Sm = (const int*)src[1-ksize];
 275             uchar* D = (uchar*)dst;
 276             if( haveScale )
 277             {
 278                 i = 0;
 279                 #if CV_SSE2
 280                 if(haveSSE2)
 281                 {
 282                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 283                     for( ; i <= width-8; i+=8 )
 284                     {
 285                         __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
 286                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 287
 288                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 289                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 290                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 291                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 292
 293                         __m128i _s0T = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 294                         __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));
 295
 296                         _s0T = _mm_packs_epi32(_s0T, _s0T1);
 297
 298                         _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));
 299
 300                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 301                         _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
 302                     }
 303                 }
 304                 #elif CV_NEON
 305                 float32x4_t v_scale = vdupq_n_f32((float)_scale);
 306                 for( ; i <= width-8; i+=8 )
 307                 {
 308                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 309                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 310
 311                     uint32x4_t v_s0d = cv_vrndq_u32_f32(vmulq_f32(vcvtq_f32_s32(v_s0), v_scale));
 312                     uint32x4_t v_s01d = cv_vrndq_u32_f32(vmulq_f32(vcvtq_f32_s32(v_s01), v_scale));
 313
 314                     uint16x8_t v_dst = vcombine_u16(vqmovn_u32(v_s0d), vqmovn_u32(v_s01d));
 315                     vst1_u8(D + i, vqmovn_u16(v_dst));
 316
 317                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 318                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 319                 }
 320                 #endif
 321                 for( ; i < width; i++ )
 322                 {
 323                     int s0 = SUM[i] + Sp[i];
 324                     D[i] = saturate_cast<uchar>(s0*_scale);
 325                     SUM[i] = s0 - Sm[i];
 326                 }
 327             }
 328             else
 329             {
 330                 i = 0;
 331                 #if CV_SSE2
 332                 if(haveSSE2)
 333                 {
 334                     for( ; i <= width-8; i+=8 )
 335                     {
 336                         __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
 337                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 338
 339                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 340                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 341                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 342                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 343
 344                         __m128i _s0T = _mm_packs_epi32(_s0, _s01);
 345
 346                         _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));
 347
 348                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 349                         _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
 350                     }
 351                 }
 352                 #elif CV_NEON
 353                 for( ; i <= width-8; i+=8 )
 354                 {
 355                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 356                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 357
 358                     uint16x8_t v_dst = vcombine_u16(vqmovun_s32(v_s0), vqmovun_s32(v_s01));
 359                     vst1_u8(D + i, vqmovn_u16(v_dst));
 360
 361                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 362                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 363                 }
 364                 #endif
 365
 366                 for( ; i < width; i++ )
 367                 {
 368                     int s0 = SUM[i] + Sp[i];
 369                     D[i] = saturate_cast<uchar>(s0);
 370                     SUM[i] = s0 - Sm[i];
 371                 }
 372             }
 373             dst += dststep;
 374         }
 375     }
 376
 377     double scale;
 378     int sumCount;
 379     std::vector<int> sum;
 380 };
 381
 382 template<>
 383 struct ColumnSum<int, short> :
 384         public BaseColumnFilter
 385 {
 386     ColumnSum( int _ksize, int _anchor, double _scale ) :
 387         BaseColumnFilter()
 388     {
 389         ksize = _ksize;
 390         anchor = _anchor;
 391         scale = _scale;
 392         sumCount = 0;
 393     }
 394
 395     virtual void reset() { sumCount = 0; }
 396
 397     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 398     {
 399         int i;
 400         int* SUM;
 401         bool haveScale = scale != 1;
 402         double _scale = scale;
 403
 404         #if CV_SSE2
 405             bool haveSSE2 = checkHardwareSupport(CV_CPU_SSE2);
 406         #endif
 407
 408         if( width != (int)sum.size() )
 409         {
 410             sum.resize(width);
 411             sumCount = 0;
 412         }
 413         SUM = &sum[0];
 414         if( sumCount == 0 )
 415         {
 416             memset((void*)SUM, 0, width*sizeof(int));
 417             for( ; sumCount < ksize - 1; sumCount++, src++ )
 418             {
 419                 const int* Sp = (const int*)src[0];
 420                 i = 0;
 421                 #if CV_SSE2
 422                 if(haveSSE2)
 423                 {
 424                     for( ; i <= width-4; i+=4 )
 425                     {
 426                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 427                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 428                         _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
 429                     }
 430                 }
 431                 #elif CV_NEON
 432                 for( ; i <= width - 4; i+=4 )
 433                     vst1q_s32(SUM + i, vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i)));
 434                 #endif
 435                 for( ; i < width; i++ )
 436                     SUM[i] += Sp[i];
 437             }
 438         }
 439         else
 440         {
 441             CV_Assert( sumCount == ksize-1 );
 442             src += ksize-1;
 443         }
 444
 445         for( ; count--; src++ )
 446         {
 447             const int* Sp = (const int*)src[0];
 448             const int* Sm = (const int*)src[1-ksize];
 449             short* D = (short*)dst;
 450             if( haveScale )
 451             {
 452                 i = 0;
 453                 #if CV_SSE2
 454                 if(haveSSE2)
 455                 {
 456                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 457                     for( ; i <= width-8; i+=8 )
 458                     {
 459                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 460                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 461
 462                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 463                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 464                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 465                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 466
 467                         __m128i _s0T  = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 468                         __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));
 469
 470                         _mm_storeu_si128((__m128i*)(D+i), _mm_packs_epi32(_s0T, _s0T1));
 471
 472                         _mm_storeu_si128((__m128i*)(SUM+i),_mm_sub_epi32(_s0,_sm));
 473                         _mm_storeu_si128((__m128i*)(SUM+i+4), _mm_sub_epi32(_s01,_sm1));
 474                     }
 475                 }
 476                 #elif CV_NEON
 477                 float32x4_t v_scale = vdupq_n_f32((float)_scale);
 478                 for( ; i <= width-8; i+=8 )
 479                 {
 480                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 481                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 482
 483                     int32x4_t v_s0d = cv_vrndq_s32_f32(vmulq_f32(vcvtq_f32_s32(v_s0), v_scale));
 484                     int32x4_t v_s01d = cv_vrndq_s32_f32(vmulq_f32(vcvtq_f32_s32(v_s01), v_scale));
 485                     vst1q_s16(D + i, vcombine_s16(vqmovn_s32(v_s0d), vqmovn_s32(v_s01d)));
 486
 487                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 488                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 489                 }
 490                 #endif
 491                 for( ; i < width; i++ )
 492                 {
 493                     int s0 = SUM[i] + Sp[i];
 494                     D[i] = saturate_cast<short>(s0*_scale);
 495                     SUM[i] = s0 - Sm[i];
 496                 }
 497             }
 498             else
 499             {
 500                 i = 0;
 501                 #if CV_SSE2
 502                 if(haveSSE2)
 503                 {
 504                     for( ; i <= width-8; i+=8 )
 505                     {
 506
 507                         __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
 508                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 509
 510                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 511                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 512                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 513                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 514
 515                         _mm_storeu_si128((__m128i*)(D+i), _mm_packs_epi32(_s0, _s01));
 516
 517                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 518                         _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
 519                     }
 520                 }
 521                 #elif CV_NEON
 522                 for( ; i <= width-8; i+=8 )
 523                 {
 524                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 525                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 526
 527                     vst1q_s16(D + i, vcombine_s16(vqmovn_s32(v_s0), vqmovn_s32(v_s01)));
 528
 529                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 530                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 531                 }
 532                 #endif
 533
 534                 for( ; i < width; i++ )
 535                 {
 536                     int s0 = SUM[i] + Sp[i];
 537                     D[i] = saturate_cast<short>(s0);
 538                     SUM[i] = s0 - Sm[i];
 539                 }
 540             }
 541             dst += dststep;
 542         }
 543     }
 544
 545     double scale;
 546     int sumCount;
 547     std::vector<int> sum;
 548 };
 549
 550
 551 template<>
 552 struct ColumnSum<int, ushort> :
 553         public BaseColumnFilter
 554 {
 555     ColumnSum( int _ksize, int _anchor, double _scale ) :
 556         BaseColumnFilter()
 557     {
 558         ksize = _ksize;
 559         anchor = _anchor;
 560         scale = _scale;
 561         sumCount = 0;
 562     }
 563
 564     virtual void reset() { sumCount = 0; }
 565
 566     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 567     {
 568         int i;
 569         int* SUM;
 570         bool haveScale = scale != 1;
 571         double _scale = scale;
 572         #if CV_SSE2
 573                 bool haveSSE2 =  checkHardwareSupport(CV_CPU_SSE2);
 574         #endif
 575
 576         if( width != (int)sum.size() )
 577         {
 578             sum.resize(width);
 579             sumCount = 0;
 580         }
 581         SUM = &sum[0];
 582         if( sumCount == 0 )
 583         {
 584             memset((void*)SUM, 0, width*sizeof(int));
 585             for( ; sumCount < ksize - 1; sumCount++, src++ )
 586             {
 587                 const int* Sp = (const int*)src[0];
 588                 i = 0;
 589                 #if CV_SSE2
 590                 if(haveSSE2)
 591                 {
 592                     for( ; i < width-4; i+=4 )
 593                     {
 594                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 595                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 596                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_add_epi32(_sum, _sp));
 597                     }
 598                 }
 599                 #elif CV_NEON
 600                 for( ; i <= width - 4; i+=4 )
 601                     vst1q_s32(SUM + i, vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i)));
 602                 #endif
 603                 for( ; i < width; i++ )
 604                     SUM[i] += Sp[i];
 605             }
 606         }
 607         else
 608         {
 609             CV_Assert( sumCount == ksize-1 );
 610             src += ksize-1;
 611         }
 612
 613         for( ; count--; src++ )
 614         {
 615             const int* Sp = (const int*)src[0];
 616             const int* Sm = (const int*)src[1-ksize];
 617             ushort* D = (ushort*)dst;
 618             if( haveScale )
 619             {
 620                 i = 0;
 621                 #if CV_SSE2
 622                 if(haveSSE2)
 623                 {
 624                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 625                     const __m128i delta0 = _mm_set1_epi32(0x8000);
 626                     const __m128i delta1 = _mm_set1_epi32(0x80008000);
 627
 628                     for( ; i < width-4; i+=4)
 629                     {
 630                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 631                         __m128i _s0   = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 632                                                       _mm_loadu_si128((const __m128i*)(Sp+i)));
 633
 634                         __m128i _res = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 635
 636                         _res = _mm_sub_epi32(_res, delta0);
 637                         _res = _mm_add_epi16(_mm_packs_epi32(_res, _res), delta1);
 638
 639                         _mm_storel_epi64((__m128i*)(D+i), _res);
 640                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 641                     }
 642                 }
 643                 #elif CV_NEON
 644                 float32x4_t v_scale = vdupq_n_f32((float)_scale);
 645                 for( ; i <= width-8; i+=8 )
 646                 {
 647                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 648                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 649
 650                     uint32x4_t v_s0d = cv_vrndq_u32_f32(vmulq_f32(vcvtq_f32_s32(v_s0), v_scale));
 651                     uint32x4_t v_s01d = cv_vrndq_u32_f32(vmulq_f32(vcvtq_f32_s32(v_s01), v_scale));
 652                     vst1q_u16(D + i, vcombine_u16(vqmovn_u32(v_s0d), vqmovn_u32(v_s01d)));
 653
 654                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 655                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 656                 }
 657                 #endif
 658                 for( ; i < width; i++ )
 659                 {
 660                     int s0 = SUM[i] + Sp[i];
 661                     D[i] = saturate_cast<ushort>(s0*_scale);
 662                     SUM[i] = s0 - Sm[i];
 663                 }
 664             }
 665             else
 666             {
 667                 i = 0;
 668                 #if  CV_SSE2
 669                 if(haveSSE2)
 670                 {
 671                     const __m128i delta0 = _mm_set1_epi32(0x8000);
 672                     const __m128i delta1 = _mm_set1_epi32(0x80008000);
 673
 674                     for( ; i < width-4; i+=4 )
 675                     {
 676                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 677                         __m128i _s0   = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 678                                                       _mm_loadu_si128((const __m128i*)(Sp+i)));
 679
 680                         __m128i _res = _mm_sub_epi32(_s0, delta0);
 681                         _res = _mm_add_epi16(_mm_packs_epi32(_res, _res), delta1);
 682
 683                         _mm_storel_epi64((__m128i*)(D+i), _res);
 684                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 685                     }
 686                 }
 687                 #elif CV_NEON
 688                 for( ; i <= width-8; i+=8 )
 689                 {
 690                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 691                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 692
 693                     vst1q_u16(D + i, vcombine_u16(vqmovun_s32(v_s0), vqmovun_s32(v_s01)));
 694
 695                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 696                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 697                 }
 698                 #endif
 699
 700                 for( ; i < width; i++ )
 701                 {
 702                     int s0 = SUM[i] + Sp[i];
 703                     D[i] = saturate_cast<ushort>(s0);
 704                     SUM[i] = s0 - Sm[i];
 705                 }
 706             }
 707             dst += dststep;
 708         }
 709     }
 710
 711     double scale;
 712     int sumCount;
 713     std::vector<int> sum;
 714 };
 715
 716 template<>
 717 struct ColumnSum<int, float> :
 718         public BaseColumnFilter
 719 {
 720     ColumnSum( int _ksize, int _anchor, double _scale ) :
 721         BaseColumnFilter()
 722     {
 723         ksize = _ksize;
 724         anchor = _anchor;
 725         scale = _scale;
 726         sumCount = 0;
 727     }
 728
 729     virtual void reset() { sumCount = 0; }
 730
 731     virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 732     {
 733         int i;
 734         int* SUM;
 735         bool haveScale = scale != 1;
 736         double _scale = scale;
 737
 738         #if CV_SSE2
 739         bool haveSSE2 =  checkHardwareSupport(CV_CPU_SSE2);
 740         #endif
 741
 742         if( width != (int)sum.size() )
 743         {
 744             sum.resize(width);
 745             sumCount = 0;
 746         }
 747
 748         SUM = &sum[0];
 749         if( sumCount == 0 )
 750         {
 751             memset((void *)SUM, 0, sizeof(int) * width);
 752
 753             for( ; sumCount < ksize - 1; sumCount++, src++ )
 754             {
 755                 const int* Sp = (const int*)src[0];
 756                 i = 0;
 757
 758                 #if CV_SSE2
 759                 if(haveSSE2)
 760                 {
 761                     for( ; i < width-4; i+=4 )
 762                     {
 763                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 764                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 765                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_add_epi32(_sum, _sp));
 766                     }
 767                 }
 768                 #elif CV_NEON
 769                 for( ; i <= width - 4; i+=4 )
 770                     vst1q_s32(SUM + i, vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i)));
 771                 #endif
 772
 773                 for( ; i < width; i++ )
 774                     SUM[i] += Sp[i];
 775             }
 776         }
 777         else
 778         {
 779             CV_Assert( sumCount == ksize-1 );
 780             src += ksize-1;
 781         }
 782
 783         for( ; count--; src++ )
 784         {
 785             const int * Sp = (const int*)src[0];
 786             const int * Sm = (const int*)src[1-ksize];
 787             float* D = (float*)dst;
 788             if( haveScale )
 789             {
 790                 i = 0;
 791
 792                 #if CV_SSE2
 793                 if(haveSSE2)
 794                 {
 795                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 796
 797                     for( ; i < width-4; i+=4)
 798                     {
 799                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 800                         __m128i _s0   = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 801                                                       _mm_loadu_si128((const __m128i*)(Sp+i)));
 802
 803                         _mm_storeu_ps(D+i, _mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 804                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 805                     }
 806                 }
 807                 #elif CV_NEON
 808                 float32x4_t v_scale = vdupq_n_f32((float)_scale);
 809                 for( ; i <= width-8; i+=8 )
 810                 {
 811                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 812                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 813
 814                     vst1q_f32(D + i, vmulq_f32(vcvtq_f32_s32(v_s0), v_scale));
 815                     vst1q_f32(D + i + 4, vmulq_f32(vcvtq_f32_s32(v_s01), v_scale));
 816
 817                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 818                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 819                 }
 820                 #endif
 821
 822                 for( ; i < width; i++ )
 823                 {
 824                     int s0 = SUM[i] + Sp[i];
 825                     D[i] = (float)(s0*_scale);
 826                     SUM[i] = s0 - Sm[i];
 827                 }
 828             }
 829             else
 830             {
 831                 i = 0;
 832
 833                 #if CV_SSE2
 834                 if(haveSSE2)
 835                 {
 836                     for( ; i < width-4; i+=4)
 837                     {
 838                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 839                         __m128i _s0   = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 840                                                       _mm_loadu_si128((const __m128i*)(Sp+i)));
 841
 842                         _mm_storeu_ps(D+i, _mm_cvtepi32_ps(_s0));
 843                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 844                     }
 845                 }
 846                 #elif CV_NEON
 847                 for( ; i <= width-8; i+=8 )
 848                 {
 849                     int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
 850                     int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));
 851
 852                     vst1q_f32(D + i, vcvtq_f32_s32(v_s0));
 853                     vst1q_f32(D + i + 4, vcvtq_f32_s32(v_s01));
 854
 855                     vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
 856                     vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
 857                 }
 858                 #endif
 859
 860                 for( ; i < width; i++ )
 861                 {
 862                     int s0 = SUM[i] + Sp[i];
 863                     D[i] = (float)(s0);
 864                     SUM[i] = s0 - Sm[i];
 865                 }
 866             }
 867             dst += dststep;
 868         }
 869     }
 870
 871     double scale;
 872     int sumCount;
 873     std::vector<int> sum;
 874 };
 875
 876 #ifdef HAVE_OPENCL
 877
 878 #define DIVUP(total, grain) ((total + grain - 1) / (grain))
 879 #define ROUNDUP(sz, n)      ((sz) + (n) - 1 - (((sz) + (n) - 1) % (n)))
 880
 881 static bool ocl_boxFilter( InputArray _src, OutputArray _dst, int ddepth,
 882                            Size ksize, Point anchor, int borderType, bool normalize, bool sqr = false )
 883 {
 884     const ocl::Device & dev = ocl::Device::getDefault();
 885     int type = _src.type(), sdepth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type), esz = CV_ELEM_SIZE(type);
 886     bool doubleSupport = dev.doubleFPConfig() > 0;
 887
 888     if (ddepth < 0)
 889         ddepth = sdepth;
 890
 891     if (cn > 4 || (!doubleSupport && (sdepth == CV_64F || ddepth == CV_64F)) ||
 892         _src.offset() % esz != 0 || _src.step() % esz != 0)
 893         return false;
 894
 895     if (anchor.x < 0)
 896         anchor.x = ksize.width / 2;
 897     if (anchor.y < 0)
 898         anchor.y = ksize.height / 2;
 899
 900     int computeUnits = ocl::Device::getDefault().maxComputeUnits();
 901     float alpha = 1.0f / (ksize.height * ksize.width);
 902     Size size = _src.size(), wholeSize;
 903     bool isolated = (borderType & BORDER_ISOLATED) != 0;
 904     borderType &= ~BORDER_ISOLATED;
 905     int wdepth = std::max(CV_32F, std::max(ddepth, sdepth)),
 906         wtype = CV_MAKE_TYPE(wdepth, cn), dtype = CV_MAKE_TYPE(ddepth, cn);
 907
 908     const char * const borderMap[] = { "BORDER_CONSTANT", "BORDER_REPLICATE", "BORDER_REFLECT", 0, "BORDER_REFLECT_101" };
 909     size_t globalsize[2] = { size.width, size.height };
 910     size_t localsize_general[2] = { 0, 1 }, * localsize = NULL;
 911
 912     UMat src = _src.getUMat();
 913     if (!isolated)
 914     {
 915         Point ofs;
 916         src.locateROI(wholeSize, ofs);
 917     }
 918
 919     int h = isolated ? size.height : wholeSize.height;
 920     int w = isolated ? size.width : wholeSize.width;
 921
 922     size_t maxWorkItemSizes[32];
 923     ocl::Device::getDefault().maxWorkItemSizes(maxWorkItemSizes);
 924     int tryWorkItems = (int)maxWorkItemSizes[0];
 925
 926     ocl::Kernel kernel;
 927
 928     if (dev.isIntel() && !(dev.type() & ocl::Device::TYPE_CPU) &&
 929         ((ksize.width < 5 && ksize.height < 5 && esz <= 4) ||
 930          (ksize.width == 5 && ksize.height == 5 && cn == 1)))
 931     {
 932         if (w < ksize.width || h < ksize.height)
 933             return false;
 934
 935         // Figure out what vector size to use for loading the pixels.
 936         int pxLoadNumPixels = cn != 1 || size.width % 4 ? 1 : 4;
 937         int pxLoadVecSize = cn * pxLoadNumPixels;
 938
 939         // Figure out how many pixels per work item to compute in X and Y
 940         // directions.  Too many and we run out of registers.
 941         int pxPerWorkItemX = 1, pxPerWorkItemY = 1;
 942         if (cn <= 2 && ksize.width <= 4 && ksize.height <= 4)
 943         {
 944             pxPerWorkItemX = size.width % 8 ? size.width % 4 ? size.width % 2 ? 1 : 2 : 4 : 8;
 945             pxPerWorkItemY = size.height % 2 ? 1 : 2;
 946         }
 947         else if (cn < 4 || (ksize.width <= 4 && ksize.height <= 4))
 948         {
 949             pxPerWorkItemX = size.width % 2 ? 1 : 2;
 950             pxPerWorkItemY = size.height % 2 ? 1 : 2;
 951         }
 952         globalsize[0] = size.width / pxPerWorkItemX;
 953         globalsize[1] = size.height / pxPerWorkItemY;
 954
 955         // Need some padding in the private array for pixels
 956         int privDataWidth = ROUNDUP(pxPerWorkItemX + ksize.width - 1, pxLoadNumPixels);
 957
 958         // Make the global size a nice round number so the runtime can pick
 959         // from reasonable choices for the workgroup size
 960         const int wgRound = 256;
 961         globalsize[0] = ROUNDUP(globalsize[0], wgRound);
 962
 963         char build_options[1024], cvt[2][40];
 964         sprintf(build_options, "-D cn=%d "
 965                 "-D ANCHOR_X=%d -D ANCHOR_Y=%d -D KERNEL_SIZE_X=%d -D KERNEL_SIZE_Y=%d "
 966                 "-D PX_LOAD_VEC_SIZE=%d -D PX_LOAD_NUM_PX=%d "
 967                 "-D PX_PER_WI_X=%d -D PX_PER_WI_Y=%d -D PRIV_DATA_WIDTH=%d -D %s -D %s "
 968                 "-D PX_LOAD_X_ITERATIONS=%d -D PX_LOAD_Y_ITERATIONS=%d "
 969                 "-D srcT=%s -D srcT1=%s -D dstT=%s -D dstT1=%s -D WT=%s -D WT1=%s "
 970                 "-D convertToWT=%s -D convertToDstT=%s%s%s -D PX_LOAD_FLOAT_VEC_CONV=convert_%s -D OP_BOX_FILTER",
 971                 cn, anchor.x, anchor.y, ksize.width, ksize.height,
 972                 pxLoadVecSize, pxLoadNumPixels,
 973                 pxPerWorkItemX, pxPerWorkItemY, privDataWidth, borderMap[borderType],
 974                 isolated ? "BORDER_ISOLATED" : "NO_BORDER_ISOLATED",
 975                 privDataWidth / pxLoadNumPixels, pxPerWorkItemY + ksize.height - 1,
 976                 ocl::typeToStr(type), ocl::typeToStr(sdepth), ocl::typeToStr(dtype),
 977                 ocl::typeToStr(ddepth), ocl::typeToStr(wtype), ocl::typeToStr(wdepth),
 978                 ocl::convertTypeStr(sdepth, wdepth, cn, cvt[0]),
 979                 ocl::convertTypeStr(wdepth, ddepth, cn, cvt[1]),
 980                 normalize ? " -D NORMALIZE" : "", sqr ? " -D SQR" : "",
 981                 ocl::typeToStr(CV_MAKE_TYPE(wdepth, pxLoadVecSize)) //PX_LOAD_FLOAT_VEC_CONV
 982                 );
 983
 984
 985         if (!kernel.create("filterSmall", cv::ocl::imgproc::filterSmall_oclsrc, build_options))
 986             return false;
 987     }
 988     else
 989     {
 990         localsize = localsize_general;
 991         for ( ; ; )
 992         {
 993             int BLOCK_SIZE_X = tryWorkItems, BLOCK_SIZE_Y = std::min(ksize.height * 10, size.height);
 994
 995             while (BLOCK_SIZE_X > 32 && BLOCK_SIZE_X >= ksize.width * 2 && BLOCK_SIZE_X > size.width * 2)
 996                 BLOCK_SIZE_X /= 2;
 997             while (BLOCK_SIZE_Y < BLOCK_SIZE_X / 8 && BLOCK_SIZE_Y * computeUnits * 32 < size.height)
 998                 BLOCK_SIZE_Y *= 2;
 999
1000             if (ksize.width > BLOCK_SIZE_X || w < ksize.width || h < ksize.height)
1001                 return false;
1002
1003             char cvt[2][50];
1004             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"
1005                                  " -D ANCHOR_X=%d -D ANCHOR_Y=%d -D KERNEL_SIZE_X=%d -D KERNEL_SIZE_Y=%d -D %s%s%s%s%s"
1006                                  " -D ST1=%s -D DT1=%s -D cn=%d",
1007                                  BLOCK_SIZE_X, BLOCK_SIZE_Y, ocl::typeToStr(type), ocl::typeToStr(CV_MAKE_TYPE(ddepth, cn)),
1008                                  ocl::typeToStr(CV_MAKE_TYPE(wdepth, cn)),
1009                                  ocl::convertTypeStr(wdepth, ddepth, cn, cvt[0]),
1010                                  ocl::convertTypeStr(sdepth, wdepth, cn, cvt[1]),
1011                                  anchor.x, anchor.y, ksize.width, ksize.height, borderMap[borderType],
1012                                  isolated ? " -D BORDER_ISOLATED" : "", doubleSupport ? " -D DOUBLE_SUPPORT" : "",
1013                                  normalize ? " -D NORMALIZE" : "", sqr ? " -D SQR" : "",
1014                                  ocl::typeToStr(sdepth), ocl::typeToStr(ddepth), cn);
1015
1016             localsize[0] = BLOCK_SIZE_X;
1017             globalsize[0] = DIVUP(size.width, BLOCK_SIZE_X - (ksize.width - 1)) * BLOCK_SIZE_X;
1018             globalsize[1] = DIVUP(size.height, BLOCK_SIZE_Y);
1019
1020             kernel.create("boxFilter", cv::ocl::imgproc::boxFilter_oclsrc, opts);
1021             if (kernel.empty())
1022                 return false;
1023
1024             size_t kernelWorkGroupSize = kernel.workGroupSize();
1025             if (localsize[0] <= kernelWorkGroupSize)
1026                 break;
1027             if (BLOCK_SIZE_X < (int)kernelWorkGroupSize)
1028                 return false;
1029
1030             tryWorkItems = (int)kernelWorkGroupSize;
1031         }
1032     }
1033
1034     _dst.create(size, CV_MAKETYPE(ddepth, cn));
1035     UMat dst = _dst.getUMat();
1036
1037     int idxArg = kernel.set(0, ocl::KernelArg::PtrReadOnly(src));
1038     idxArg = kernel.set(idxArg, (int)src.step);
1039     int srcOffsetX = (int)((src.offset % src.step) / src.elemSize());
1040     int srcOffsetY = (int)(src.offset / src.step);
1041     int srcEndX = isolated ? srcOffsetX + size.width : wholeSize.width;
1042     int srcEndY = isolated ? srcOffsetY + size.height : wholeSize.height;
1043     idxArg = kernel.set(idxArg, srcOffsetX);
1044     idxArg = kernel.set(idxArg, srcOffsetY);
1045     idxArg = kernel.set(idxArg, srcEndX);
1046     idxArg = kernel.set(idxArg, srcEndY);
1047     idxArg = kernel.set(idxArg, ocl::KernelArg::WriteOnly(dst));
1048     if (normalize)
1049         idxArg = kernel.set(idxArg, (float)alpha);
1050
1051     return kernel.run(2, globalsize, localsize, false);
1052 }
1053
1054 #undef ROUNDUP
1055
1056 #endif
1057
1058 }
1059
1060
1061 cv::Ptr<cv::BaseRowFilter> cv::getRowSumFilter(int srcType, int sumType, int ksize, int anchor)
1062 {
1063     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
1064     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
1065
1066     if( anchor < 0 )
1067         anchor = ksize/2;
1068
1069     if( sdepth == CV_8U && ddepth == CV_32S )
1070         return makePtr<RowSum<uchar, int> >(ksize, anchor);
1071     if( sdepth == CV_8U && ddepth == CV_64F )
1072         return makePtr<RowSum<uchar, double> >(ksize, anchor);
1073     if( sdepth == CV_16U && ddepth == CV_32S )
1074         return makePtr<RowSum<ushort, int> >(ksize, anchor);
1075     if( sdepth == CV_16U && ddepth == CV_64F )
1076         return makePtr<RowSum<ushort, double> >(ksize, anchor);
1077     if( sdepth == CV_16S && ddepth == CV_32S )
1078         return makePtr<RowSum<short, int> >(ksize, anchor);
1079     if( sdepth == CV_32S && ddepth == CV_32S )
1080         return makePtr<RowSum<int, int> >(ksize, anchor);
1081     if( sdepth == CV_16S && ddepth == CV_64F )
1082         return makePtr<RowSum<short, double> >(ksize, anchor);
1083     if( sdepth == CV_32F && ddepth == CV_64F )
1084         return makePtr<RowSum<float, double> >(ksize, anchor);
1085     if( sdepth == CV_64F && ddepth == CV_64F )
1086         return makePtr<RowSum<double, double> >(ksize, anchor);
1087
1088     CV_Error_( CV_StsNotImplemented,
1089         ("Unsupported combination of source format (=%d), and buffer format (=%d)",
1090         srcType, sumType));
1091
1092     return Ptr<BaseRowFilter>();
1093 }
1094
1095
1096 cv::Ptr<cv::BaseColumnFilter> cv::getColumnSumFilter(int sumType, int dstType, int ksize,
1097                                                      int anchor, double scale)
1098 {
1099     int sdepth = CV_MAT_DEPTH(sumType), ddepth = CV_MAT_DEPTH(dstType);
1100     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(dstType) );
1101
1102     if( anchor < 0 )
1103         anchor = ksize/2;
1104
1105     if( ddepth == CV_8U && sdepth == CV_32S )
1106         return makePtr<ColumnSum<int, uchar> >(ksize, anchor, scale);
1107     if( ddepth == CV_8U && sdepth == CV_64F )
1108         return makePtr<ColumnSum<double, uchar> >(ksize, anchor, scale);
1109     if( ddepth == CV_16U && sdepth == CV_32S )
1110         return makePtr<ColumnSum<int, ushort> >(ksize, anchor, scale);
1111     if( ddepth == CV_16U && sdepth == CV_64F )
1112         return makePtr<ColumnSum<double, ushort> >(ksize, anchor, scale);
1113     if( ddepth == CV_16S && sdepth == CV_32S )
1114         return makePtr<ColumnSum<int, short> >(ksize, anchor, scale);
1115     if( ddepth == CV_16S && sdepth == CV_64F )
1116         return makePtr<ColumnSum<double, short> >(ksize, anchor, scale);
1117     if( ddepth == CV_32S && sdepth == CV_32S )
1118         return makePtr<ColumnSum<int, int> >(ksize, anchor, scale);
1119     if( ddepth == CV_32F && sdepth == CV_32S )
1120         return makePtr<ColumnSum<int, float> >(ksize, anchor, scale);
1121     if( ddepth == CV_32F && sdepth == CV_64F )
1122         return makePtr<ColumnSum<double, float> >(ksize, anchor, scale);
1123     if( ddepth == CV_64F && sdepth == CV_32S )
1124         return makePtr<ColumnSum<int, double> >(ksize, anchor, scale);
1125     if( ddepth == CV_64F && sdepth == CV_64F )
1126         return makePtr<ColumnSum<double, double> >(ksize, anchor, scale);
1127
1128     CV_Error_( CV_StsNotImplemented,
1129         ("Unsupported combination of sum format (=%d), and destination format (=%d)",
1130         sumType, dstType));
1131
1132     return Ptr<BaseColumnFilter>();
1133 }
1134
1135
1136 cv::Ptr<cv::FilterEngine> cv::createBoxFilter( int srcType, int dstType, Size ksize,
1137                     Point anchor, bool normalize, int borderType )
1138 {
1139     int sdepth = CV_MAT_DEPTH(srcType);
1140     int cn = CV_MAT_CN(srcType), sumType = CV_64F;
1141     if( sdepth <= CV_32S && (!normalize ||
1142         ksize.width*ksize.height <= (sdepth == CV_8U ? (1<<23) :
1143             sdepth == CV_16U ? (1 << 15) : (1 << 16))) )
1144         sumType = CV_32S;
1145     sumType = CV_MAKETYPE( sumType, cn );
1146
1147     Ptr<BaseRowFilter> rowFilter = getRowSumFilter(srcType, sumType, ksize.width, anchor.x );
1148     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
1149         dstType, ksize.height, anchor.y, normalize ? 1./(ksize.width*ksize.height) : 1);
1150
1151     return makePtr<FilterEngine>(Ptr<BaseFilter>(), rowFilter, columnFilter,
1152            srcType, dstType, sumType, borderType );
1153 }
1154
1155
1156 void cv::boxFilter( InputArray _src, OutputArray _dst, int ddepth,
1157                 Size ksize, Point anchor,
1158                 bool normalize, int borderType )
1159 {
1160     CV_OCL_RUN(_dst.isUMat(), ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize))
1161
1162     Mat src = _src.getMat();
1163     int stype = src.type(), sdepth = CV_MAT_DEPTH(stype), cn = CV_MAT_CN(stype);
1164     if( ddepth < 0 )
1165         ddepth = sdepth;
1166     _dst.create( src.size(), CV_MAKETYPE(ddepth, cn) );
1167     Mat dst = _dst.getMat();
1168     if( borderType != BORDER_CONSTANT && normalize && (borderType & BORDER_ISOLATED) != 0 )
1169     {
1170         if( src.rows == 1 )
1171             ksize.height = 1;
1172         if( src.cols == 1 )
1173             ksize.width = 1;
1174     }
1175 #ifdef HAVE_TEGRA_OPTIMIZATION
1176     if ( tegra::box(src, dst, ksize, anchor, normalize, borderType) )
1177         return;
1178 #endif
1179
1180 #if defined(HAVE_IPP)
1181     CV_IPP_CHECK()
1182     {
1183         int ippBorderType = borderType & ~BORDER_ISOLATED;
1184         Point ocvAnchor, ippAnchor;
1185         ocvAnchor.x = anchor.x < 0 ? ksize.width / 2 : anchor.x;
1186         ocvAnchor.y = anchor.y < 0 ? ksize.height / 2 : anchor.y;
1187         ippAnchor.x = ksize.width / 2 - (ksize.width % 2 == 0 ? 1 : 0);
1188         ippAnchor.y = ksize.height / 2 - (ksize.height % 2 == 0 ? 1 : 0);
1189
1190         if (normalize && !src.isSubmatrix() && ddepth == sdepth &&
1191             (/*ippBorderType == BORDER_REPLICATE ||*/ /* returns ippStsStepErr: Step value is not valid */
1192              ippBorderType == BORDER_CONSTANT) && ocvAnchor == ippAnchor &&
1193              dst.cols != ksize.width && dst.rows != ksize.height) // returns ippStsMaskSizeErr: mask has an illegal value
1194         {
1195             Ipp32s bufSize = 0;
1196             IppiSize roiSize = { dst.cols, dst.rows }, maskSize = { ksize.width, ksize.height };
1197
1198 #define IPP_FILTER_BOX_BORDER(ippType, ippDataType, flavor) \
1199             do \
1200             { \
1201                 if (ippiFilterBoxBorderGetBufferSize(roiSize, maskSize, ippDataType, cn, &bufSize) >= 0) \
1202                 { \
1203                     Ipp8u * buffer = ippsMalloc_8u(bufSize); \
1204                     ippType borderValue[4] = { 0, 0, 0, 0 }; \
1205                     ippBorderType = ippBorderType == BORDER_CONSTANT ? ippBorderConst : ippBorderRepl; \
1206                     IppStatus status = ippiFilterBoxBorder_##flavor(src.ptr<ippType>(), (int)src.step, dst.ptr<ippType>(), \
1207                                                                     (int)dst.step, roiSize, maskSize, \
1208                                                                     (IppiBorderType)ippBorderType, borderValue, buffer); \
1209                     ippsFree(buffer); \
1210                     if (status >= 0) \
1211                     { \
1212                         CV_IMPL_ADD(CV_IMPL_IPP); \
1213                         return; \
1214                     } \
1215                 } \
1216                 setIppErrorStatus(); \
1217             } while ((void)0, 0)
1218
1219             if (stype == CV_8UC1)
1220                 IPP_FILTER_BOX_BORDER(Ipp8u, ipp8u, 8u_C1R);
1221             else if (stype == CV_8UC3)
1222                 IPP_FILTER_BOX_BORDER(Ipp8u, ipp8u, 8u_C3R);
1223             else if (stype == CV_8UC4)
1224                 IPP_FILTER_BOX_BORDER(Ipp8u, ipp8u, 8u_C4R);
1225
1226             else if (stype == CV_16UC1)
1227                 IPP_FILTER_BOX_BORDER(Ipp16u, ipp16u, 16u_C1R);
1228             else if (stype == CV_16UC3)
1229                 IPP_FILTER_BOX_BORDER(Ipp16u, ipp16u, 16u_C3R);
1230             else if (stype == CV_16UC4)
1231                 IPP_FILTER_BOX_BORDER(Ipp16u, ipp16u, 16u_C4R);
1232
1233             else if (stype == CV_16SC1)
1234                 IPP_FILTER_BOX_BORDER(Ipp16s, ipp16s, 16s_C1R);
1235             else if (stype == CV_16SC3)
1236                 IPP_FILTER_BOX_BORDER(Ipp16s, ipp16s, 16s_C3R);
1237             else if (stype == CV_16SC4)
1238                 IPP_FILTER_BOX_BORDER(Ipp16s, ipp16s, 16s_C4R);
1239
1240             else if (stype == CV_32FC1)
1241                 IPP_FILTER_BOX_BORDER(Ipp32f, ipp32f, 32f_C1R);
1242             else if (stype == CV_32FC3)
1243                 IPP_FILTER_BOX_BORDER(Ipp32f, ipp32f, 32f_C3R);
1244             else if (stype == CV_32FC4)
1245                 IPP_FILTER_BOX_BORDER(Ipp32f, ipp32f, 32f_C4R);
1246         }
1247 #undef IPP_FILTER_BOX_BORDER
1248     }
1249 #endif
1250
1251     Ptr<FilterEngine> f = createBoxFilter( src.type(), dst.type(),
1252                         ksize, anchor, normalize, borderType );
1253     f->apply( src, dst );
1254 }
1255
1256 void cv::blur( InputArray src, OutputArray dst,
1257            Size ksize, Point anchor, int borderType )
1258 {
1259     boxFilter( src, dst, -1, ksize, anchor, true, borderType );
1260 }
1261
1262
1263 /****************************************************************************************\
1264                                     Squared Box Filter
1265 \****************************************************************************************/
1266
1267 namespace cv
1268 {
1269
1270 template<typename T, typename ST>
1271 struct SqrRowSum :
1272         public BaseRowFilter
1273 {
1274     SqrRowSum( int _ksize, int _anchor ) :
1275         BaseRowFilter()
1276     {
1277         ksize = _ksize;
1278         anchor = _anchor;
1279     }
1280
1281     virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
1282     {
1283         const T* S = (const T*)src;
1284         ST* D = (ST*)dst;
1285         int i = 0, k, ksz_cn = ksize*cn;
1286
1287         width = (width - 1)*cn;
1288         for( k = 0; k < cn; k++, S++, D++ )
1289         {
1290             ST s = 0;
1291             for( i = 0; i < ksz_cn; i += cn )
1292             {
1293                 ST val = (ST)S[i];
1294                 s += val*val;
1295             }
1296             D[0] = s;
1297             for( i = 0; i < width; i += cn )
1298             {
1299                 ST val0 = (ST)S[i], val1 = (ST)S[i + ksz_cn];
1300                 s += val1*val1 - val0*val0;
1301                 D[i+cn] = s;
1302             }
1303         }
1304     }
1305 };
1306
1307 static Ptr<BaseRowFilter> getSqrRowSumFilter(int srcType, int sumType, int ksize, int anchor)
1308 {
1309     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
1310     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
1311
1312     if( anchor < 0 )
1313         anchor = ksize/2;
1314
1315     if( sdepth == CV_8U && ddepth == CV_32S )
1316         return makePtr<SqrRowSum<uchar, int> >(ksize, anchor);
1317     if( sdepth == CV_8U && ddepth == CV_64F )
1318         return makePtr<SqrRowSum<uchar, double> >(ksize, anchor);
1319     if( sdepth == CV_16U && ddepth == CV_64F )
1320         return makePtr<SqrRowSum<ushort, double> >(ksize, anchor);
1321     if( sdepth == CV_16S && ddepth == CV_64F )
1322         return makePtr<SqrRowSum<short, double> >(ksize, anchor);
1323     if( sdepth == CV_32F && ddepth == CV_64F )
1324         return makePtr<SqrRowSum<float, double> >(ksize, anchor);
1325     if( sdepth == CV_64F && ddepth == CV_64F )
1326         return makePtr<SqrRowSum<double, double> >(ksize, anchor);
1327
1328     CV_Error_( CV_StsNotImplemented,
1329               ("Unsupported combination of source format (=%d), and buffer format (=%d)",
1330                srcType, sumType));
1331
1332     return Ptr<BaseRowFilter>();
1333 }
1334
1335 }
1336
1337 void cv::sqrBoxFilter( InputArray _src, OutputArray _dst, int ddepth,
1338                        Size ksize, Point anchor,
1339                        bool normalize, int borderType )
1340 {
1341     int srcType = _src.type(), sdepth = CV_MAT_DEPTH(srcType), cn = CV_MAT_CN(srcType);
1342     Size size = _src.size();
1343
1344     if( ddepth < 0 )
1345         ddepth = sdepth < CV_32F ? CV_32F : CV_64F;
1346
1347     if( borderType != BORDER_CONSTANT && normalize )
1348     {
1349         if( size.height == 1 )
1350             ksize.height = 1;
1351         if( size.width == 1 )
1352             ksize.width = 1;
1353     }
1354
1355     CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2,
1356                ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize, true))
1357
1358     int sumDepth = CV_64F;
1359     if( sdepth == CV_8U )
1360         sumDepth = CV_32S;
1361     int sumType = CV_MAKETYPE( sumDepth, cn ), dstType = CV_MAKETYPE(ddepth, cn);
1362
1363     Mat src = _src.getMat();
1364     _dst.create( size, dstType );
1365     Mat dst = _dst.getMat();
1366
1367     Ptr<BaseRowFilter> rowFilter = getSqrRowSumFilter(srcType, sumType, ksize.width, anchor.x );
1368     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
1369                                                             dstType, ksize.height, anchor.y,
1370                                                             normalize ? 1./(ksize.width*ksize.height) : 1);
1371
1372     Ptr<FilterEngine> f = makePtr<FilterEngine>(Ptr<BaseFilter>(), rowFilter, columnFilter,
1373                                                 srcType, dstType, sumType, borderType );
1374     f->apply( src, dst );
1375 }
1376
1377
1378 /****************************************************************************************\
1379                                      Gaussian Blur
1380 \****************************************************************************************/
1381
1382 cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )
1383 {
1384     const int SMALL_GAUSSIAN_SIZE = 7;
1385     static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
1386     {
1387         {1.f},
1388         {0.25f, 0.5f, 0.25f},
1389         {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
1390         {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
1391     };
1392
1393     const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
1394         small_gaussian_tab[n>>1] : 0;
1395
1396     CV_Assert( ktype == CV_32F || ktype == CV_64F );
1397     Mat kernel(n, 1, ktype);
1398     float* cf = kernel.ptr<float>();
1399     double* cd = kernel.ptr<double>();
1400
1401     double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
1402     double scale2X = -0.5/(sigmaX*sigmaX);
1403     double sum = 0;
1404
1405     int i;
1406     for( i = 0; i < n; i++ )
1407     {
1408         double x = i - (n-1)*0.5;
1409         double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);
1410         if( ktype == CV_32F )
1411         {
1412             cf[i] = (float)t;
1413             sum += cf[i];
1414         }
1415         else
1416         {
1417             cd[i] = t;
1418             sum += cd[i];
1419         }
1420     }
1421
1422     sum = 1./sum;
1423     for( i = 0; i < n; i++ )
1424     {
1425         if( ktype == CV_32F )
1426             cf[i] = (float)(cf[i]*sum);
1427         else
1428             cd[i] *= sum;
1429     }
1430
1431     return kernel;
1432 }
1433
1434 namespace cv {
1435
1436 static void createGaussianKernels( Mat & kx, Mat & ky, int type, Size ksize,
1437                                    double sigma1, double sigma2 )
1438 {
1439     int depth = CV_MAT_DEPTH(type);
1440     if( sigma2 <= 0 )
1441         sigma2 = sigma1;
1442
1443     // automatic detection of kernel size from sigma
1444     if( ksize.width <= 0 && sigma1 > 0 )
1445         ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
1446     if( ksize.height <= 0 && sigma2 > 0 )
1447         ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
1448
1449     CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 &&
1450         ksize.height > 0 && ksize.height % 2 == 1 );
1451
1452     sigma1 = std::max( sigma1, 0. );
1453     sigma2 = std::max( sigma2, 0. );
1454
1455     kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) );
1456     if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )
1457         ky = kx;
1458     else
1459         ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) );
1460 }
1461
1462 }
1463
1464 cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize,
1465                                         double sigma1, double sigma2,
1466                                         int borderType )
1467 {
1468     Mat kx, ky;
1469     createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
1470
1471     return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );
1472 }
1473
1474
1475 void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
1476                    double sigma1, double sigma2,
1477                    int borderType )
1478 {
1479     int type = _src.type();
1480     Size size = _src.size();
1481     _dst.create( size, type );
1482
1483     if( borderType != BORDER_CONSTANT && (borderType & BORDER_ISOLATED) != 0 )
1484     {
1485         if( size.height == 1 )
1486             ksize.height = 1;
1487         if( size.width == 1 )
1488             ksize.width = 1;
1489     }
1490
1491     if( ksize.width == 1 && ksize.height == 1 )
1492     {
1493         _src.copyTo(_dst);
1494         return;
1495     }
1496
1497 #ifdef HAVE_TEGRA_OPTIMIZATION
1498     if(sigma1 == 0 && sigma2 == 0 && tegra::gaussian(_src.getMat(), _dst.getMat(), ksize, borderType))
1499         return;
1500 #endif
1501
1502 #if IPP_VERSION_X100 >= 801 && 0 // these functions are slower in IPP 8.1
1503     CV_IPP_CHECK()
1504     {
1505         int depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
1506
1507         if ((depth == CV_8U || depth == CV_16U || depth == CV_16S || depth == CV_32F) && (cn == 1 || cn == 3) &&
1508                 sigma1 == sigma2 && ksize.width == ksize.height && sigma1 != 0.0 )
1509         {
1510             IppiBorderType ippBorder = ippiGetBorderType(borderType);
1511             if (ippBorderConst == ippBorder || ippBorderRepl == ippBorder)
1512             {
1513                 Mat src = _src.getMat(), dst = _dst.getMat();
1514                 IppiSize roiSize = { src.cols, src.rows };
1515                 IppDataType dataType = ippiGetDataType(depth);
1516                 Ipp32s specSize = 0, bufferSize = 0;
1517
1518                 if (ippiFilterGaussianGetBufferSize(roiSize, (Ipp32u)ksize.width, dataType, cn, &specSize, &bufferSize) >= 0)
1519                 {
1520                     IppFilterGaussianSpec * pSpec = (IppFilterGaussianSpec *)ippMalloc(specSize);
1521                     Ipp8u * pBuffer = (Ipp8u*)ippMalloc(bufferSize);
1522
1523                     if (ippiFilterGaussianInit(roiSize, (Ipp32u)ksize.width, (Ipp32f)sigma1, ippBorder, dataType, 1, pSpec, pBuffer) >= 0)
1524                     {
1525 #define IPP_FILTER_GAUSS(ippfavor, ippcn) \
1526         do \
1527         { \
1528             typedef Ipp##ippfavor ippType; \
1529             ippType borderValues[] = { 0, 0, 0 }; \
1530             IppStatus status = ippcn == 1 ? \
1531                 ippiFilterGaussianBorder_##ippfavor##_C1R(src.ptr<ippType>(), (int)src.step, \
1532                     dst.ptr<ippType>(), (int)dst.step, roiSize, borderValues[0], pSpec, pBuffer) : \
1533                 ippiFilterGaussianBorder_##ippfavor##_C3R(src.ptr<ippType>(), (int)src.step, \
1534                     dst.ptr<ippType>(), (int)dst.step, roiSize, borderValues, pSpec, pBuffer); \
1535             ippFree(pBuffer); \
1536             ippFree(pSpec); \
1537             if (status >= 0) \
1538             { \
1539                 CV_IMPL_ADD(CV_IMPL_IPP); \
1540                 return; \
1541             } \
1542         } while ((void)0, 0)
1543
1544                         if (type == CV_8UC1)
1545                             IPP_FILTER_GAUSS(8u, 1);
1546                         else if (type == CV_8UC3)
1547                             IPP_FILTER_GAUSS(8u, 3);
1548                         else if (type == CV_16UC1)
1549                             IPP_FILTER_GAUSS(16u, 1);
1550                         else if (type == CV_16UC3)
1551                             IPP_FILTER_GAUSS(16u, 3);
1552                         else if (type == CV_16SC1)
1553                             IPP_FILTER_GAUSS(16s, 1);
1554                         else if (type == CV_16SC3)
1555                             IPP_FILTER_GAUSS(16s, 3);
1556                         else if (type == CV_32FC1)
1557                             IPP_FILTER_GAUSS(32f, 1);
1558                         else if (type == CV_32FC3)
1559                             IPP_FILTER_GAUSS(32f, 3);
1560 #undef IPP_FILTER_GAUSS
1561                     }
1562                 }
1563                 setIppErrorStatus();
1564             }
1565         }
1566     }
1567 #endif
1568
1569     Mat kx, ky;
1570     createGaussianKernels(kx, ky, type, ksize, sigma1, sigma2);
1571     sepFilter2D(_src, _dst, CV_MAT_DEPTH(type), kx, ky, Point(-1,-1), 0, borderType );
1572 }
1573
1574 /****************************************************************************************\
1575                                       Median Filter
1576 \****************************************************************************************/
1577
1578 namespace cv
1579 {
1580 typedef ushort HT;
1581
1582 /**
1583  * This structure represents a two-tier histogram. The first tier (known as the
1584  * "coarse" level) is 4 bit wide and the second tier (known as the "fine" level)
1585  * is 8 bit wide. Pixels inserted in the fine level also get inserted into the
1586  * coarse bucket designated by the 4 MSBs of the fine bucket value.
1587  *
1588  * The structure is aligned on 16 bits, which is a prerequisite for SIMD
1589  * instructions. Each bucket is 16 bit wide, which means that extra care must be
1590  * taken to prevent overflow.
1591  */
1592 typedef struct
1593 {
1594     HT coarse[16];
1595     HT fine[16][16];
1596 } Histogram;
1597
1598
1599 #if CV_SSE2
1600 #define MEDIAN_HAVE_SIMD 1
1601
1602 static inline void histogram_add_simd( const HT x[16], HT y[16] )
1603 {
1604     const __m128i* rx = (const __m128i*)x;
1605     __m128i* ry = (__m128i*)y;
1606     __m128i r0 = _mm_add_epi16(_mm_load_si128(ry+0),_mm_load_si128(rx+0));
1607     __m128i r1 = _mm_add_epi16(_mm_load_si128(ry+1),_mm_load_si128(rx+1));
1608     _mm_store_si128(ry+0, r0);
1609     _mm_store_si128(ry+1, r1);
1610 }
1611
1612 static inline void histogram_sub_simd( const HT x[16], HT y[16] )
1613 {
1614     const __m128i* rx = (const __m128i*)x;
1615     __m128i* ry = (__m128i*)y;
1616     __m128i r0 = _mm_sub_epi16(_mm_load_si128(ry+0),_mm_load_si128(rx+0));
1617     __m128i r1 = _mm_sub_epi16(_mm_load_si128(ry+1),_mm_load_si128(rx+1));
1618     _mm_store_si128(ry+0, r0);
1619     _mm_store_si128(ry+1, r1);
1620 }
1621
1622 #elif CV_NEON
1623 #define MEDIAN_HAVE_SIMD 1
1624
1625 static inline void histogram_add_simd( const HT x[16], HT y[16] )
1626 {
1627     vst1q_u16(y, vaddq_u16(vld1q_u16(x), vld1q_u16(y)));
1628     vst1q_u16(y + 8, vaddq_u16(vld1q_u16(x + 8), vld1q_u16(y + 8)));
1629 }
1630
1631 static inline void histogram_sub_simd( const HT x[16], HT y[16] )
1632 {
1633     vst1q_u16(y, vsubq_u16(vld1q_u16(x), vld1q_u16(y)));
1634     vst1q_u16(y + 8, vsubq_u16(vld1q_u16(x + 8), vld1q_u16(y + 8)));
1635 }
1636
1637 #else
1638 #define MEDIAN_HAVE_SIMD 0
1639 #endif
1640
1641
1642 static inline void histogram_add( const HT x[16], HT y[16] )
1643 {
1644     int i;
1645     for( i = 0; i < 16; ++i )
1646         y[i] = (HT)(y[i] + x[i]);
1647 }
1648
1649 static inline void histogram_sub( const HT x[16], HT y[16] )
1650 {
1651     int i;
1652     for( i = 0; i < 16; ++i )
1653         y[i] = (HT)(y[i] - x[i]);
1654 }
1655
1656 static inline void histogram_muladd( int a, const HT x[16],
1657         HT y[16] )
1658 {
1659     for( int i = 0; i < 16; ++i )
1660         y[i] = (HT)(y[i] + a * x[i]);
1661 }
1662
1663 static void
1664 medianBlur_8u_O1( const Mat& _src, Mat& _dst, int ksize )
1665 {
1666 /**
1667  * HOP is short for Histogram OPeration. This macro makes an operation \a op on
1668  * histogram \a h for pixel value \a x. It takes care of handling both levels.
1669  */
1670 #define HOP(h,x,op) \
1671     h.coarse[x>>4] op, \
1672     *((HT*)h.fine + x) op
1673
1674 #define COP(c,j,x,op) \
1675     h_coarse[ 16*(n*c+j) + (x>>4) ] op, \
1676     h_fine[ 16 * (n*(16*c+(x>>4)) + j) + (x & 0xF) ] op
1677
1678     int cn = _dst.channels(), m = _dst.rows, r = (ksize-1)/2;
1679     size_t sstep = _src.step, dstep = _dst.step;
1680     Histogram CV_DECL_ALIGNED(16) H[4];
1681     HT CV_DECL_ALIGNED(16) luc[4][16];
1682
1683     int STRIPE_SIZE = std::min( _dst.cols, 512/cn );
1684
1685     std::vector<HT> _h_coarse(1 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
1686     std::vector<HT> _h_fine(16 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
1687     HT* h_coarse = alignPtr(&_h_coarse[0], 16);
1688     HT* h_fine = alignPtr(&_h_fine[0], 16);
1689 #if MEDIAN_HAVE_SIMD
1690     volatile bool useSIMD = checkHardwareSupport(CV_CPU_SSE2) || checkHardwareSupport(CV_CPU_NEON);
1691 #endif
1692
1693     for( int x = 0; x < _dst.cols; x += STRIPE_SIZE )
1694     {
1695         int i, j, k, c, n = std::min(_dst.cols - x, STRIPE_SIZE) + r*2;
1696         const uchar* src = _src.ptr() + x*cn;
1697         uchar* dst = _dst.ptr() + (x - r)*cn;
1698
1699         memset( h_coarse, 0, 16*n*cn*sizeof(h_coarse[0]) );
1700         memset( h_fine, 0, 16*16*n*cn*sizeof(h_fine[0]) );
1701
1702         // First row initialization
1703         for( c = 0; c < cn; c++ )
1704         {
1705             for( j = 0; j < n; j++ )
1706                 COP( c, j, src[cn*j+c], += (cv::HT)(r+2) );
1707
1708             for( i = 1; i < r; i++ )
1709             {
1710                 const uchar* p = src + sstep*std::min(i, m-1);
1711                 for ( j = 0; j < n; j++ )
1712                     COP( c, j, p[cn*j+c], ++ );
1713             }
1714         }
1715
1716         for( i = 0; i < m; i++ )
1717         {
1718             const uchar* p0 = src + sstep * std::max( 0, i-r-1 );
1719             const uchar* p1 = src + sstep * std::min( m-1, i+r );
1720
1721             memset( H, 0, cn*sizeof(H[0]) );
1722             memset( luc, 0, cn*sizeof(luc[0]) );
1723             for( c = 0; c < cn; c++ )
1724             {
1725                 // Update column histograms for the entire row.
1726                 for( j = 0; j < n; j++ )
1727                 {
1728                     COP( c, j, p0[j*cn + c], -- );
1729                     COP( c, j, p1[j*cn + c], ++ );
1730                 }
1731
1732                 // First column initialization
1733                 for( k = 0; k < 16; ++k )
1734                     histogram_muladd( 2*r+1, &h_fine[16*n*(16*c+k)], &H[c].fine[k][0] );
1735
1736             #if MEDIAN_HAVE_SIMD
1737                 if( useSIMD )
1738                 {
1739                     for( j = 0; j < 2*r; ++j )
1740                         histogram_add_simd( &h_coarse[16*(n*c+j)], H[c].coarse );
1741
1742                     for( j = r; j < n-r; j++ )
1743                     {
1744                         int t = 2*r*r + 2*r, b, sum = 0;
1745                         HT* segment;
1746
1747                         histogram_add_simd( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
1748
1749                         // Find median at coarse level
1750                         for ( k = 0; k < 16 ; ++k )
1751                         {
1752                             sum += H[c].coarse[k];
1753                             if ( sum > t )
1754                             {
1755                                 sum -= H[c].coarse[k];
1756                                 break;
1757                             }
1758                         }
1759                         assert( k < 16 );
1760
1761                         /* Update corresponding histogram segment */
1762                         if ( luc[c][k] <= j-r )
1763                         {
1764                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
1765                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
1766                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
1767
1768                             if ( luc[c][k] < j+r+1 )
1769                             {
1770                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
1771                                 luc[c][k] = (HT)(j+r+1);
1772                             }
1773                         }
1774                         else
1775                         {
1776                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
1777                             {
1778                                 histogram_sub_simd( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
1779                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
1780                             }
1781                         }
1782
1783                         histogram_sub_simd( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
1784
1785                         /* Find median in segment */
1786                         segment = H[c].fine[k];
1787                         for ( b = 0; b < 16 ; b++ )
1788                         {
1789                             sum += segment[b];
1790                             if ( sum > t )
1791                             {
1792                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
1793                                 break;
1794                             }
1795                         }
1796                         assert( b < 16 );
1797                     }
1798                 }
1799                 else
1800             #endif
1801                 {
1802                     for( j = 0; j < 2*r; ++j )
1803                         histogram_add( &h_coarse[16*(n*c+j)], H[c].coarse );
1804
1805                     for( j = r; j < n-r; j++ )
1806                     {
1807                         int t = 2*r*r + 2*r, b, sum = 0;
1808                         HT* segment;
1809
1810                         histogram_add( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
1811
1812                         // Find median at coarse level
1813                         for ( k = 0; k < 16 ; ++k )
1814                         {
1815                             sum += H[c].coarse[k];
1816                             if ( sum > t )
1817                             {
1818                                 sum -= H[c].coarse[k];
1819                                 break;
1820                             }
1821                         }
1822                         assert( k < 16 );
1823
1824                         /* Update corresponding histogram segment */
1825                         if ( luc[c][k] <= j-r )
1826                         {
1827                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
1828                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
1829                                 histogram_add( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
1830
1831                             if ( luc[c][k] < j+r+1 )
1832                             {
1833                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
1834                                 luc[c][k] = (HT)(j+r+1);
1835                             }
1836                         }
1837                         else
1838                         {
1839                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
1840                             {
1841                                 histogram_sub( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
1842                                 histogram_add( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
1843                             }
1844                         }
1845
1846                         histogram_sub( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
1847
1848                         /* Find median in segment */
1849                         segment = H[c].fine[k];
1850                         for ( b = 0; b < 16 ; b++ )
1851                         {
1852                             sum += segment[b];
1853                             if ( sum > t )
1854                             {
1855                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
1856                                 break;
1857                             }
1858                         }
1859                         assert( b < 16 );
1860                     }
1861                 }
1862             }
1863         }
1864     }
1865
1866 #undef HOP
1867 #undef COP
1868 }
1869
1870 static void
1871 medianBlur_8u_Om( const Mat& _src, Mat& _dst, int m )
1872 {
1873     #define N  16
1874     int     zone0[4][N];
1875     int     zone1[4][N*N];
1876     int     x, y;
1877     int     n2 = m*m/2;
1878     Size    size = _dst.size();
1879     const uchar* src = _src.ptr();
1880     uchar*  dst = _dst.ptr();
1881     int     src_step = (int)_src.step, dst_step = (int)_dst.step;
1882     int     cn = _src.channels();
1883     const uchar*  src_max = src + size.height*src_step;
1884
1885     #define UPDATE_ACC01( pix, cn, op ) \
1886     {                                   \
1887         int p = (pix);                  \
1888         zone1[cn][p] op;                \
1889         zone0[cn][p >> 4] op;           \
1890     }
1891
1892     //CV_Assert( size.height >= nx && size.width >= nx );
1893     for( x = 0; x < size.width; x++, src += cn, dst += cn )
1894     {
1895         uchar* dst_cur = dst;
1896         const uchar* src_top = src;
1897         const uchar* src_bottom = src;
1898         int k, c;
1899         int src_step1 = src_step, dst_step1 = dst_step;
1900
1901         if( x % 2 != 0 )
1902         {
1903             src_bottom = src_top += src_step*(size.height-1);
1904             dst_cur += dst_step*(size.height-1);
1905             src_step1 = -src_step1;
1906             dst_step1 = -dst_step1;
1907         }
1908
1909         // init accumulator
1910         memset( zone0, 0, sizeof(zone0[0])*cn );
1911         memset( zone1, 0, sizeof(zone1[0])*cn );
1912
1913         for( y = 0; y <= m/2; y++ )
1914         {
1915             for( c = 0; c < cn; c++ )
1916             {
1917                 if( y > 0 )
1918                 {
1919                     for( k = 0; k < m*cn; k += cn )
1920                         UPDATE_ACC01( src_bottom[k+c], c, ++ );
1921                 }
1922                 else
1923                 {
1924                     for( k = 0; k < m*cn; k += cn )
1925                         UPDATE_ACC01( src_bottom[k+c], c, += m/2+1 );
1926                 }
1927             }
1928
1929             if( (src_step1 > 0 && y < size.height-1) ||
1930                 (src_step1 < 0 && size.height-y-1 > 0) )
1931                 src_bottom += src_step1;
1932         }
1933
1934         for( y = 0; y < size.height; y++, dst_cur += dst_step1 )
1935         {
1936             // find median
1937             for( c = 0; c < cn; c++ )
1938             {
1939                 int s = 0;
1940                 for( k = 0; ; k++ )
1941                 {
1942                     int t = s + zone0[c][k];
1943                     if( t > n2 ) break;
1944                     s = t;
1945                 }
1946
1947                 for( k *= N; ;k++ )
1948                 {
1949                     s += zone1[c][k];
1950                     if( s > n2 ) break;
1951                 }
1952
1953                 dst_cur[c] = (uchar)k;
1954             }
1955
1956             if( y+1 == size.height )
1957                 break;
1958
1959             if( cn == 1 )
1960             {
1961                 for( k = 0; k < m; k++ )
1962                 {
1963                     int p = src_top[k];
1964                     int q = src_bottom[k];
1965                     zone1[0][p]--;
1966                     zone0[0][p>>4]--;
1967                     zone1[0][q]++;
1968                     zone0[0][q>>4]++;
1969                 }
1970             }
1971             else if( cn == 3 )
1972             {
1973                 for( k = 0; k < m*3; k += 3 )
1974                 {
1975                     UPDATE_ACC01( src_top[k], 0, -- );
1976                     UPDATE_ACC01( src_top[k+1], 1, -- );
1977                     UPDATE_ACC01( src_top[k+2], 2, -- );
1978
1979                     UPDATE_ACC01( src_bottom[k], 0, ++ );
1980                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
1981                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
1982                 }
1983             }
1984             else
1985             {
1986                 assert( cn == 4 );
1987                 for( k = 0; k < m*4; k += 4 )
1988                 {
1989                     UPDATE_ACC01( src_top[k], 0, -- );
1990                     UPDATE_ACC01( src_top[k+1], 1, -- );
1991                     UPDATE_ACC01( src_top[k+2], 2, -- );
1992                     UPDATE_ACC01( src_top[k+3], 3, -- );
1993
1994                     UPDATE_ACC01( src_bottom[k], 0, ++ );
1995                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
1996                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
1997                     UPDATE_ACC01( src_bottom[k+3], 3, ++ );
1998                 }
1999             }
2000
2001             if( (src_step1 > 0 && src_bottom + src_step1 < src_max) ||
2002                 (src_step1 < 0 && src_bottom + src_step1 >= src) )
2003                 src_bottom += src_step1;
2004
2005             if( y >= m/2 )
2006                 src_top += src_step1;
2007         }
2008     }
2009 #undef N
2010 #undef UPDATE_ACC
2011 }
2012
2013
2014 struct MinMax8u
2015 {
2016     typedef uchar value_type;
2017     typedef int arg_type;
2018     enum { SIZE = 1 };
2019     arg_type load(const uchar* ptr) { return *ptr; }
2020     void store(uchar* ptr, arg_type val) { *ptr = (uchar)val; }
2021     void operator()(arg_type& a, arg_type& b) const
2022     {
2023         int t = CV_FAST_CAST_8U(a - b);
2024         b += t; a -= t;
2025     }
2026 };
2027
2028 struct MinMax16u
2029 {
2030     typedef ushort value_type;
2031     typedef int arg_type;
2032     enum { SIZE = 1 };
2033     arg_type load(const ushort* ptr) { return *ptr; }
2034     void store(ushort* ptr, arg_type val) { *ptr = (ushort)val; }
2035     void operator()(arg_type& a, arg_type& b) const
2036     {
2037         arg_type t = a;
2038         a = std::min(a, b);
2039         b = std::max(b, t);
2040     }
2041 };
2042
2043 struct MinMax16s
2044 {
2045     typedef short value_type;
2046     typedef int arg_type;
2047     enum { SIZE = 1 };
2048     arg_type load(const short* ptr) { return *ptr; }
2049     void store(short* ptr, arg_type val) { *ptr = (short)val; }
2050     void operator()(arg_type& a, arg_type& b) const
2051     {
2052         arg_type t = a;
2053         a = std::min(a, b);
2054         b = std::max(b, t);
2055     }
2056 };
2057
2058 struct MinMax32f
2059 {
2060     typedef float value_type;
2061     typedef float arg_type;
2062     enum { SIZE = 1 };
2063     arg_type load(const float* ptr) { return *ptr; }
2064     void store(float* ptr, arg_type val) { *ptr = val; }
2065     void operator()(arg_type& a, arg_type& b) const
2066     {
2067         arg_type t = a;
2068         a = std::min(a, b);
2069         b = std::max(b, t);
2070     }
2071 };
2072
2073 #if CV_SSE2
2074
2075 struct MinMaxVec8u
2076 {
2077     typedef uchar value_type;
2078     typedef __m128i arg_type;
2079     enum { SIZE = 16 };
2080     arg_type load(const uchar* ptr) { return _mm_loadu_si128((const __m128i*)ptr); }
2081     void store(uchar* ptr, arg_type val) { _mm_storeu_si128((__m128i*)ptr, val); }
2082     void operator()(arg_type& a, arg_type& b) const
2083     {
2084         arg_type t = a;
2085         a = _mm_min_epu8(a, b);
2086         b = _mm_max_epu8(b, t);
2087     }
2088 };
2089
2090
2091 struct MinMaxVec16u
2092 {
2093     typedef ushort value_type;
2094     typedef __m128i arg_type;
2095     enum { SIZE = 8 };
2096     arg_type load(const ushort* ptr) { return _mm_loadu_si128((const __m128i*)ptr); }
2097     void store(ushort* ptr, arg_type val) { _mm_storeu_si128((__m128i*)ptr, val); }
2098     void operator()(arg_type& a, arg_type& b) const
2099     {
2100         arg_type t = _mm_subs_epu16(a, b);
2101         a = _mm_subs_epu16(a, t);
2102         b = _mm_adds_epu16(b, t);
2103     }
2104 };
2105
2106
2107 struct MinMaxVec16s
2108 {
2109     typedef short value_type;
2110     typedef __m128i arg_type;
2111     enum { SIZE = 8 };
2112     arg_type load(const short* ptr) { return _mm_loadu_si128((const __m128i*)ptr); }
2113     void store(short* ptr, arg_type val) { _mm_storeu_si128((__m128i*)ptr, val); }
2114     void operator()(arg_type& a, arg_type& b) const
2115     {
2116         arg_type t = a;
2117         a = _mm_min_epi16(a, b);
2118         b = _mm_max_epi16(b, t);
2119     }
2120 };
2121
2122
2123 struct MinMaxVec32f
2124 {
2125     typedef float value_type;
2126     typedef __m128 arg_type;
2127     enum { SIZE = 4 };
2128     arg_type load(const float* ptr) { return _mm_loadu_ps(ptr); }
2129     void store(float* ptr, arg_type val) { _mm_storeu_ps(ptr, val); }
2130     void operator()(arg_type& a, arg_type& b) const
2131     {
2132         arg_type t = a;
2133         a = _mm_min_ps(a, b);
2134         b = _mm_max_ps(b, t);
2135     }
2136 };
2137
2138 #elif CV_NEON
2139
2140 struct MinMaxVec8u
2141 {
2142     typedef uchar value_type;
2143     typedef uint8x16_t arg_type;
2144     enum { SIZE = 16 };
2145     arg_type load(const uchar* ptr) { return vld1q_u8(ptr); }
2146     void store(uchar* ptr, arg_type val) { vst1q_u8(ptr, val); }
2147     void operator()(arg_type& a, arg_type& b) const
2148     {
2149         arg_type t = a;
2150         a = vminq_u8(a, b);
2151         b = vmaxq_u8(b, t);
2152     }
2153 };
2154
2155
2156 struct MinMaxVec16u
2157 {
2158     typedef ushort value_type;
2159     typedef uint16x8_t arg_type;
2160     enum { SIZE = 8 };
2161     arg_type load(const ushort* ptr) { return vld1q_u16(ptr); }
2162     void store(ushort* ptr, arg_type val) { vst1q_u16(ptr, val); }
2163     void operator()(arg_type& a, arg_type& b) const
2164     {
2165         arg_type t = a;
2166         a = vminq_u16(a, b);
2167         b = vmaxq_u16(b, t);
2168     }
2169 };
2170
2171
2172 struct MinMaxVec16s
2173 {
2174     typedef short value_type;
2175     typedef int16x8_t arg_type;
2176     enum { SIZE = 8 };
2177     arg_type load(const short* ptr) { return vld1q_s16(ptr); }
2178     void store(short* ptr, arg_type val) { vst1q_s16(ptr, val); }
2179     void operator()(arg_type& a, arg_type& b) const
2180     {
2181         arg_type t = a;
2182         a = vminq_s16(a, b);
2183         b = vmaxq_s16(b, t);
2184     }
2185 };
2186
2187
2188 struct MinMaxVec32f
2189 {
2190     typedef float value_type;
2191     typedef float32x4_t arg_type;
2192     enum { SIZE = 4 };
2193     arg_type load(const float* ptr) { return vld1q_f32(ptr); }
2194     void store(float* ptr, arg_type val) { vst1q_f32(ptr, val); }
2195     void operator()(arg_type& a, arg_type& b) const
2196     {
2197         arg_type t = a;
2198         a = vminq_f32(a, b);
2199         b = vmaxq_f32(b, t);
2200     }
2201 };
2202
2203
2204 #else
2205
2206 typedef MinMax8u MinMaxVec8u;
2207 typedef MinMax16u MinMaxVec16u;
2208 typedef MinMax16s MinMaxVec16s;
2209 typedef MinMax32f MinMaxVec32f;
2210
2211 #endif
2212
2213 template<class Op, class VecOp>
2214 static void
2215 medianBlur_SortNet( const Mat& _src, Mat& _dst, int m )
2216 {
2217     typedef typename Op::value_type T;
2218     typedef typename Op::arg_type WT;
2219     typedef typename VecOp::arg_type VT;
2220
2221     const T* src = _src.ptr<T>();
2222     T* dst = _dst.ptr<T>();
2223     int sstep = (int)(_src.step/sizeof(T));
2224     int dstep = (int)(_dst.step/sizeof(T));
2225     Size size = _dst.size();
2226     int i, j, k, cn = _src.channels();
2227     Op op;
2228     VecOp vop;
2229     volatile bool useSIMD = checkHardwareSupport(CV_CPU_SSE2) || checkHardwareSupport(CV_CPU_NEON);
2230
2231     if( m == 3 )
2232     {
2233         if( size.width == 1 || size.height == 1 )
2234         {
2235             int len = size.width + size.height - 1;
2236             int sdelta = size.height == 1 ? cn : sstep;
2237             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
2238             int ddelta = size.height == 1 ? cn : dstep;
2239
2240             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
2241                 for( j = 0; j < cn; j++, src++ )
2242                 {
2243                     WT p0 = src[i > 0 ? -sdelta : 0];
2244                     WT p1 = src[0];
2245                     WT p2 = src[i < len - 1 ? sdelta : 0];
2246
2247                     op(p0, p1); op(p1, p2); op(p0, p1);
2248                     dst[j] = (T)p1;
2249                 }
2250             return;
2251         }
2252
2253         size.width *= cn;
2254         for( i = 0; i < size.height; i++, dst += dstep )
2255         {
2256             const T* row0 = src + std::max(i - 1, 0)*sstep;
2257             const T* row1 = src + i*sstep;
2258             const T* row2 = src + std::min(i + 1, size.height-1)*sstep;
2259             int limit = useSIMD ? cn : size.width;
2260
2261             for(j = 0;; )
2262             {
2263                 for( ; j < limit; j++ )
2264                 {
2265                     int j0 = j >= cn ? j - cn : j;
2266                     int j2 = j < size.width - cn ? j + cn : j;
2267                     WT p0 = row0[j0], p1 = row0[j], p2 = row0[j2];
2268                     WT p3 = row1[j0], p4 = row1[j], p5 = row1[j2];
2269                     WT p6 = row2[j0], p7 = row2[j], p8 = row2[j2];
2270
2271                     op(p1, p2); op(p4, p5); op(p7, p8); op(p0, p1);
2272                     op(p3, p4); op(p6, p7); op(p1, p2); op(p4, p5);
2273                     op(p7, p8); op(p0, p3); op(p5, p8); op(p4, p7);
2274                     op(p3, p6); op(p1, p4); op(p2, p5); op(p4, p7);
2275                     op(p4, p2); op(p6, p4); op(p4, p2);
2276                     dst[j] = (T)p4;
2277                 }
2278
2279                 if( limit == size.width )
2280                     break;
2281
2282                 for( ; j <= size.width - VecOp::SIZE - cn; j += VecOp::SIZE )
2283                 {
2284                     VT p0 = vop.load(row0+j-cn), p1 = vop.load(row0+j), p2 = vop.load(row0+j+cn);
2285                     VT p3 = vop.load(row1+j-cn), p4 = vop.load(row1+j), p5 = vop.load(row1+j+cn);
2286                     VT p6 = vop.load(row2+j-cn), p7 = vop.load(row2+j), p8 = vop.load(row2+j+cn);
2287
2288                     vop(p1, p2); vop(p4, p5); vop(p7, p8); vop(p0, p1);
2289                     vop(p3, p4); vop(p6, p7); vop(p1, p2); vop(p4, p5);
2290                     vop(p7, p8); vop(p0, p3); vop(p5, p8); vop(p4, p7);
2291                     vop(p3, p6); vop(p1, p4); vop(p2, p5); vop(p4, p7);
2292                     vop(p4, p2); vop(p6, p4); vop(p4, p2);
2293                     vop.store(dst+j, p4);
2294                 }
2295
2296                 limit = size.width;
2297             }
2298         }
2299     }
2300     else if( m == 5 )
2301     {
2302         if( size.width == 1 || size.height == 1 )
2303         {
2304             int len = size.width + size.height - 1;
2305             int sdelta = size.height == 1 ? cn : sstep;
2306             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
2307             int ddelta = size.height == 1 ? cn : dstep;
2308
2309             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
2310                 for( j = 0; j < cn; j++, src++ )
2311                 {
2312                     int i1 = i > 0 ? -sdelta : 0;
2313                     int i0 = i > 1 ? -sdelta*2 : i1;
2314                     int i3 = i < len-1 ? sdelta : 0;
2315                     int i4 = i < len-2 ? sdelta*2 : i3;
2316                     WT p0 = src[i0], p1 = src[i1], p2 = src[0], p3 = src[i3], p4 = src[i4];
2317
2318                     op(p0, p1); op(p3, p4); op(p2, p3); op(p3, p4); op(p0, p2);
2319                     op(p2, p4); op(p1, p3); op(p1, p2);
2320                     dst[j] = (T)p2;
2321                 }
2322             return;
2323         }
2324
2325         size.width *= cn;
2326         for( i = 0; i < size.height; i++, dst += dstep )
2327         {
2328             const T* row[5];
2329             row[0] = src + std::max(i - 2, 0)*sstep;
2330             row[1] = src + std::max(i - 1, 0)*sstep;
2331             row[2] = src + i*sstep;
2332             row[3] = src + std::min(i + 1, size.height-1)*sstep;
2333             row[4] = src + std::min(i + 2, size.height-1)*sstep;
2334             int limit = useSIMD ? cn*2 : size.width;
2335
2336             for(j = 0;; )
2337             {
2338                 for( ; j < limit; j++ )
2339                 {
2340                     WT p[25];
2341                     int j1 = j >= cn ? j - cn : j;
2342                     int j0 = j >= cn*2 ? j - cn*2 : j1;
2343                     int j3 = j < size.width - cn ? j + cn : j;
2344                     int j4 = j < size.width - cn*2 ? j + cn*2 : j3;
2345                     for( k = 0; k < 5; k++ )
2346                     {
2347                         const T* rowk = row[k];
2348                         p[k*5] = rowk[j0]; p[k*5+1] = rowk[j1];
2349                         p[k*5+2] = rowk[j]; p[k*5+3] = rowk[j3];
2350                         p[k*5+4] = rowk[j4];
2351                     }
2352
2353                     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]);
2354                     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]);
2355                     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]);
2356                     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]);
2357                     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]);
2358                     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]);
2359                     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]);
2360                     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]);
2361                     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]);
2362                     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]);
2363                     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]);
2364                     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]);
2365                     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]);
2366                     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]);
2367                     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]);
2368                     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]);
2369                     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]);
2370                     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]);
2371                     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]);
2372                     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]);
2373                     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]);
2374                     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]);
2375                     op(p[7], p[11]); op(p[11], p[13]); op(p[11], p[12]);
2376                     dst[j] = (T)p[12];
2377                 }
2378
2379                 if( limit == size.width )
2380                     break;
2381
2382                 for( ; j <= size.width - VecOp::SIZE - cn*2; j += VecOp::SIZE )
2383                 {
2384                     VT p[25];
2385                     for( k = 0; k < 5; k++ )
2386                     {
2387                         const T* rowk = row[k];
2388                         p[k*5] = vop.load(rowk+j-cn*2); p[k*5+1] = vop.load(rowk+j-cn);
2389                         p[k*5+2] = vop.load(rowk+j); p[k*5+3] = vop.load(rowk+j+cn);
2390                         p[k*5+4] = vop.load(rowk+j+cn*2);
2391                     }
2392
2393                     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]);
2394                     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]);
2395                     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]);
2396                     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]);
2397                     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]);
2398                     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]);
2399                     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]);
2400                     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]);
2401                     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]);
2402                     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]);
2403                     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]);
2404                     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]);
2405                     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]);
2406                     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]);
2407                     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]);
2408                     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]);
2409                     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]);
2410                     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]);
2411                     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]);
2412                     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]);
2413                     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]);
2414                     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]);
2415                     vop(p[7], p[11]); vop(p[11], p[13]); vop(p[11], p[12]);
2416                     vop.store(dst+j, p[12]);
2417                 }
2418
2419                 limit = size.width;
2420             }
2421         }
2422     }
2423 }
2424
2425 #ifdef HAVE_OPENCL
2426
2427 static bool ocl_medianFilter(InputArray _src, OutputArray _dst, int m)
2428 {
2429     size_t localsize[2] = { 16, 16 };
2430     size_t globalsize[2];
2431     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2432
2433     if ( !((depth == CV_8U || depth == CV_16U || depth == CV_16S || depth == CV_32F) && cn <= 4 && (m == 3 || m == 5)) )
2434         return false;
2435
2436     Size imgSize = _src.size();
2437     bool useOptimized = (1 == cn) &&
2438                         (size_t)imgSize.width >= localsize[0] * 8  &&
2439                         (size_t)imgSize.height >= localsize[1] * 8 &&
2440                         imgSize.width % 4 == 0 &&
2441                         imgSize.height % 4 == 0 &&
2442                         (ocl::Device::getDefault().isIntel());
2443
2444     cv::String kname = format( useOptimized ? "medianFilter%d_u" : "medianFilter%d", m) ;
2445     cv::String kdefs = useOptimized ?
2446                          format("-D T=%s -D T1=%s -D T4=%s%d -D cn=%d -D USE_4OPT", ocl::typeToStr(type),
2447                          ocl::typeToStr(depth), ocl::typeToStr(depth), cn*4, cn)
2448                          :
2449                          format("-D T=%s -D T1=%s -D cn=%d", ocl::typeToStr(type), ocl::typeToStr(depth), cn) ;
2450
2451     ocl::Kernel k(kname.c_str(), ocl::imgproc::medianFilter_oclsrc, kdefs.c_str() );
2452
2453     if (k.empty())
2454         return false;
2455
2456     UMat src = _src.getUMat();
2457     _dst.create(src.size(), type);
2458     UMat dst = _dst.getUMat();
2459
2460     k.args(ocl::KernelArg::ReadOnlyNoSize(src), ocl::KernelArg::WriteOnly(dst));
2461
2462     if( useOptimized )
2463     {
2464         globalsize[0] = DIVUP(src.cols / 4, localsize[0]) * localsize[0];
2465         globalsize[1] = DIVUP(src.rows / 4, localsize[1]) * localsize[1];
2466     }
2467     else
2468     {
2469         globalsize[0] = (src.cols + localsize[0] + 2) / localsize[0] * localsize[0];
2470         globalsize[1] = (src.rows + localsize[1] - 1) / localsize[1] * localsize[1];
2471     }
2472
2473     return k.run(2, globalsize, localsize, false);
2474 }
2475
2476 #endif
2477
2478 }
2479
2480 void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize )
2481 {
2482     CV_Assert( (ksize % 2 == 1) && (_src0.dims() <= 2 ));
2483
2484     if( ksize <= 1 )
2485     {
2486         _src0.copyTo(_dst);
2487         return;
2488     }
2489
2490     CV_OCL_RUN(_dst.isUMat(),
2491                ocl_medianFilter(_src0,_dst, ksize))
2492
2493     Mat src0 = _src0.getMat();
2494     _dst.create( src0.size(), src0.type() );
2495     Mat dst = _dst.getMat();
2496
2497 #if IPP_VERSION_X100 >= 801
2498     CV_IPP_CHECK()
2499     {
2500 #define IPP_FILTER_MEDIAN_BORDER(ippType, ippDataType, flavor) \
2501         do \
2502         { \
2503             if (ippiFilterMedianBorderGetBufferSize(dstRoiSize, maskSize, \
2504                 ippDataType, CV_MAT_CN(type), &bufSize) >= 0) \
2505             { \
2506                 Ipp8u * buffer = ippsMalloc_8u(bufSize); \
2507                 IppStatus status = ippiFilterMedianBorder_##flavor(src.ptr<ippType>(), (int)src.step, \
2508                     dst.ptr<ippType>(), (int)dst.step, dstRoiSize, maskSize, \
2509                     ippBorderRepl, (ippType)0, buffer); \
2510                 ippsFree(buffer); \
2511                 if (status >= 0) \
2512                 { \
2513                     CV_IMPL_ADD(CV_IMPL_IPP); \
2514                     return; \
2515                 } \
2516             } \
2517             setIppErrorStatus(); \
2518         } \
2519         while ((void)0, 0)
2520
2521         if( ksize <= 5 )
2522         {
2523             Ipp32s bufSize;
2524             IppiSize dstRoiSize = ippiSize(dst.cols, dst.rows), maskSize = ippiSize(ksize, ksize);
2525             Mat src;
2526             if( dst.data != src0.data )
2527                 src = src0;
2528             else
2529                 src0.copyTo(src);
2530
2531             int type = src0.type();
2532             if (type == CV_8UC1)
2533                 IPP_FILTER_MEDIAN_BORDER(Ipp8u, ipp8u, 8u_C1R);
2534             else if (type == CV_16UC1)
2535                 IPP_FILTER_MEDIAN_BORDER(Ipp16u, ipp16u, 16u_C1R);
2536             else if (type == CV_16SC1)
2537                 IPP_FILTER_MEDIAN_BORDER(Ipp16s, ipp16s, 16s_C1R);
2538             else if (type == CV_32FC1)
2539                 IPP_FILTER_MEDIAN_BORDER(Ipp32f, ipp32f, 32f_C1R);
2540         }
2541 #undef IPP_FILTER_MEDIAN_BORDER
2542     }
2543 #endif
2544
2545 #ifdef HAVE_TEGRA_OPTIMIZATION
2546     if (tegra::medianBlur(src0, dst, ksize))
2547         return;
2548 #endif
2549
2550     bool useSortNet = ksize == 3 || (ksize == 5
2551 #if !(CV_SSE2 || CV_NEON)
2552             && src0.depth() > CV_8U
2553 #endif
2554         );
2555
2556     Mat src;
2557     if( useSortNet )
2558     {
2559         if( dst.data != src0.data )
2560             src = src0;
2561         else
2562             src0.copyTo(src);
2563
2564         if( src.depth() == CV_8U )
2565             medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize );
2566         else if( src.depth() == CV_16U )
2567             medianBlur_SortNet<MinMax16u, MinMaxVec16u>( src, dst, ksize );
2568         else if( src.depth() == CV_16S )
2569             medianBlur_SortNet<MinMax16s, MinMaxVec16s>( src, dst, ksize );
2570         else if( src.depth() == CV_32F )
2571             medianBlur_SortNet<MinMax32f, MinMaxVec32f>( src, dst, ksize );
2572         else
2573             CV_Error(CV_StsUnsupportedFormat, "");
2574
2575         return;
2576     }
2577     else
2578     {
2579         cv::copyMakeBorder( src0, src, 0, 0, ksize/2, ksize/2, BORDER_REPLICATE );
2580
2581         int cn = src0.channels();
2582         CV_Assert( src.depth() == CV_8U && (cn == 1 || cn == 3 || cn == 4) );
2583
2584         double img_size_mp = (double)(src0.total())/(1 << 20);
2585         if( ksize <= 3 + (img_size_mp < 1 ? 12 : img_size_mp < 4 ? 6 : 2)*
2586             (MEDIAN_HAVE_SIMD && (checkHardwareSupport(CV_CPU_SSE2) || checkHardwareSupport(CV_CPU_NEON)) ? 1 : 3))
2587             medianBlur_8u_Om( src, dst, ksize );
2588         else
2589             medianBlur_8u_O1( src, dst, ksize );
2590     }
2591 }
2592
2593 /****************************************************************************************\
2594                                    Bilateral Filtering
2595 \****************************************************************************************/
2596
2597 namespace cv
2598 {
2599
2600 class BilateralFilter_8u_Invoker :
2601     public ParallelLoopBody
2602 {
2603 public:
2604     BilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, int _radius, int _maxk,
2605         int* _space_ofs, float *_space_weight, float *_color_weight) :
2606         temp(&_temp), dest(&_dest), radius(_radius),
2607         maxk(_maxk), space_ofs(_space_ofs), space_weight(_space_weight), color_weight(_color_weight)
2608     {
2609     }
2610
2611     virtual void operator() (const Range& range) const
2612     {
2613         int i, j, cn = dest->channels(), k;
2614         Size size = dest->size();
2615         #if CV_SSE3
2616         int CV_DECL_ALIGNED(16) buf[4];
2617         float CV_DECL_ALIGNED(16) bufSum[4];
2618         static const int CV_DECL_ALIGNED(16) bufSignMask[] = { 0x80000000, 0x80000000, 0x80000000, 0x80000000 };
2619         bool haveSSE3 = checkHardwareSupport(CV_CPU_SSE3);
2620         #endif
2621
2622         for( i = range.start; i < range.end; i++ )
2623         {
2624             const uchar* sptr = temp->ptr(i+radius) + radius*cn;
2625             uchar* dptr = dest->ptr(i);
2626
2627             if( cn == 1 )
2628             {
2629                 for( j = 0; j < size.width; j++ )
2630                 {
2631                     float sum = 0, wsum = 0;
2632                     int val0 = sptr[j];
2633                     k = 0;
2634                     #if CV_SSE3
2635                     if( haveSSE3 )
2636                     {
2637                         __m128 _val0 = _mm_set1_ps(static_cast<float>(val0));
2638                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
2639
2640                         for( ; k <= maxk - 4; k += 4 )
2641                         {
2642                             __m128 _valF = _mm_set_ps(sptr[j + space_ofs[k+3]], sptr[j + space_ofs[k+2]],
2643                                                       sptr[j + space_ofs[k+1]], sptr[j + space_ofs[k]]);
2644
2645                             __m128 _val = _mm_andnot_ps(_signMask, _mm_sub_ps(_valF, _val0));
2646                             _mm_store_si128((__m128i*)buf, _mm_cvtps_epi32(_val));
2647
2648                             __m128 _cw = _mm_set_ps(color_weight[buf[3]],color_weight[buf[2]],
2649                                                     color_weight[buf[1]],color_weight[buf[0]]);
2650                             __m128 _sw = _mm_loadu_ps(space_weight+k);
2651                             __m128 _w = _mm_mul_ps(_cw, _sw);
2652                              _cw = _mm_mul_ps(_w, _valF);
2653
2654                              _sw = _mm_hadd_ps(_w, _cw);
2655                              _sw = _mm_hadd_ps(_sw, _sw);
2656                              _mm_storel_pi((__m64*)bufSum, _sw);
2657
2658                              sum += bufSum[1];
2659                              wsum += bufSum[0];
2660                         }
2661                     }
2662                     #endif
2663                     for( ; k < maxk; k++ )
2664                     {
2665                         int val = sptr[j + space_ofs[k]];
2666                         float w = space_weight[k]*color_weight[std::abs(val - val0)];
2667                         sum += val*w;
2668                         wsum += w;
2669                     }
2670                     // overflow is not possible here => there is no need to use cv::saturate_cast
2671                     dptr[j] = (uchar)cvRound(sum/wsum);
2672                 }
2673             }
2674             else
2675             {
2676                 assert( cn == 3 );
2677                 for( j = 0; j < size.width*3; j += 3 )
2678                 {
2679                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
2680                     int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
2681                     k = 0;
2682                     #if CV_SSE3
2683                     if( haveSSE3 )
2684                     {
2685                         const __m128i izero = _mm_setzero_si128();
2686                         const __m128 _b0 = _mm_set1_ps(static_cast<float>(b0));
2687                         const __m128 _g0 = _mm_set1_ps(static_cast<float>(g0));
2688                         const __m128 _r0 = _mm_set1_ps(static_cast<float>(r0));
2689                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
2690
2691                         for( ; k <= maxk - 4; k += 4 )
2692                         {
2693                             const int* const sptr_k0  = reinterpret_cast<const int*>(sptr + j + space_ofs[k]);
2694                             const int* const sptr_k1  = reinterpret_cast<const int*>(sptr + j + space_ofs[k+1]);
2695                             const int* const sptr_k2  = reinterpret_cast<const int*>(sptr + j + space_ofs[k+2]);
2696                             const int* const sptr_k3  = reinterpret_cast<const int*>(sptr + j + space_ofs[k+3]);
2697
2698                             __m128 _b = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_cvtsi32_si128(sptr_k0[0]), izero), izero));
2699                             __m128 _g = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_cvtsi32_si128(sptr_k1[0]), izero), izero));
2700                             __m128 _r = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_cvtsi32_si128(sptr_k2[0]), izero), izero));
2701                             __m128 _z = _mm_cvtepi32_ps(_mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_cvtsi32_si128(sptr_k3[0]), izero), izero));
2702
2703                             _MM_TRANSPOSE4_PS(_b, _g, _r, _z);
2704
2705                             __m128 bt = _mm_andnot_ps(_signMask, _mm_sub_ps(_b,_b0));
2706                             __m128 gt = _mm_andnot_ps(_signMask, _mm_sub_ps(_g,_g0));
2707                             __m128 rt = _mm_andnot_ps(_signMask, _mm_sub_ps(_r,_r0));
2708
2709                             bt =_mm_add_ps(rt, _mm_add_ps(bt, gt));
2710                             _mm_store_si128((__m128i*)buf, _mm_cvtps_epi32(bt));
2711
2712                             __m128 _w  = _mm_set_ps(color_weight[buf[3]],color_weight[buf[2]],
2713                                                     color_weight[buf[1]],color_weight[buf[0]]);
2714                             __m128 _sw = _mm_loadu_ps(space_weight+k);
2715
2716                             _w = _mm_mul_ps(_w,_sw);
2717                             _b = _mm_mul_ps(_b, _w);
2718                             _g = _mm_mul_ps(_g, _w);
2719                             _r = _mm_mul_ps(_r, _w);
2720
2721                              _w = _mm_hadd_ps(_w, _b);
2722                              _g = _mm_hadd_ps(_g, _r);
2723
2724                              _w = _mm_hadd_ps(_w, _g);
2725                              _mm_store_ps(bufSum, _w);
2726
2727                              wsum  += bufSum[0];
2728                              sum_b += bufSum[1];
2729                              sum_g += bufSum[2];
2730                              sum_r += bufSum[3];
2731                          }
2732                     }
2733                     #endif
2734
2735                     for( ; k < maxk; k++ )
2736                     {
2737                         const uchar* sptr_k = sptr + j + space_ofs[k];
2738                         int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
2739                         float w = space_weight[k]*color_weight[std::abs(b - b0) +
2740                                                                std::abs(g - g0) + std::abs(r - r0)];
2741                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
2742                         wsum += w;
2743                     }
2744                     wsum = 1.f/wsum;
2745                     b0 = cvRound(sum_b*wsum);
2746                     g0 = cvRound(sum_g*wsum);
2747                     r0 = cvRound(sum_r*wsum);
2748                     dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
2749                 }
2750             }
2751         }
2752     }
2753
2754 private:
2755     const Mat *temp;
2756     Mat *dest;
2757     int radius, maxk, *space_ofs;
2758     float *space_weight, *color_weight;
2759 };
2760
2761 #if defined (HAVE_IPP) && !defined(HAVE_IPP_ICV_ONLY) && 0
2762 class IPPBilateralFilter_8u_Invoker :
2763     public ParallelLoopBody
2764 {
2765 public:
2766     IPPBilateralFilter_8u_Invoker(Mat &_src, Mat &_dst, double _sigma_color, double _sigma_space, int _radius, bool *_ok) :
2767       ParallelLoopBody(), src(_src), dst(_dst), sigma_color(_sigma_color), sigma_space(_sigma_space), radius(_radius), ok(_ok)
2768       {
2769           *ok = true;
2770       }
2771
2772       virtual void operator() (const Range& range) const
2773       {
2774           int d = radius * 2 + 1;
2775           IppiSize kernel = {d, d};
2776           IppiSize roi={dst.cols, range.end - range.start};
2777           int bufsize=0;
2778           if (0 > ippiFilterBilateralGetBufSize_8u_C1R( ippiFilterBilateralGauss, roi, kernel, &bufsize))
2779           {
2780               *ok = false;
2781               return;
2782           }
2783           AutoBuffer<uchar> buf(bufsize);
2784           IppiFilterBilateralSpec *pSpec = (IppiFilterBilateralSpec *)alignPtr(&buf[0], 32);
2785           if (0 > ippiFilterBilateralInit_8u_C1R( ippiFilterBilateralGauss, kernel, (Ipp32f)sigma_color, (Ipp32f)sigma_space, 1, pSpec ))
2786           {
2787               *ok = false;
2788               return;
2789           }
2790           if (0 > ippiFilterBilateral_8u_C1R( src.ptr<uchar>(range.start) + radius * ((int)src.step[0] + 1), (int)src.step[0], dst.ptr<uchar>(range.start), (int)dst.step[0], roi, kernel, pSpec ))
2791               *ok = false;
2792           else
2793           {
2794             CV_IMPL_ADD(CV_IMPL_IPP|CV_IMPL_MT);
2795           }
2796       }
2797 private:
2798     Mat &src;
2799     Mat &dst;
2800     double sigma_color;
2801     double sigma_space;
2802     int radius;
2803     bool *ok;
2804     const IPPBilateralFilter_8u_Invoker& operator= (const IPPBilateralFilter_8u_Invoker&);
2805 };
2806 #endif
2807
2808 #ifdef HAVE_OPENCL
2809
2810 static bool ocl_bilateralFilter_8u(InputArray _src, OutputArray _dst, int d,
2811                                    double sigma_color, double sigma_space,
2812                                    int borderType)
2813 {
2814     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2815     int i, j, maxk, radius;
2816
2817     if (depth != CV_8U || cn > 4)
2818         return false;
2819
2820     if (sigma_color <= 0)
2821         sigma_color = 1;
2822     if (sigma_space <= 0)
2823         sigma_space = 1;
2824
2825     double gauss_color_coeff = -0.5 / (sigma_color * sigma_color);
2826     double gauss_space_coeff = -0.5 / (sigma_space * sigma_space);
2827
2828     if ( d <= 0 )
2829         radius = cvRound(sigma_space * 1.5);
2830     else
2831         radius = d / 2;
2832     radius = MAX(radius, 1);
2833     d = radius * 2 + 1;
2834
2835     UMat src = _src.getUMat(), dst = _dst.getUMat(), temp;
2836     if (src.u == dst.u)
2837         return false;
2838
2839     copyMakeBorder(src, temp, radius, radius, radius, radius, borderType);
2840     std::vector<float> _space_weight(d * d);
2841     std::vector<int> _space_ofs(d * d);
2842     float * const space_weight = &_space_weight[0];
2843     int * const space_ofs = &_space_ofs[0];
2844
2845     // initialize space-related bilateral filter coefficients
2846     for( i = -radius, maxk = 0; i <= radius; i++ )
2847         for( j = -radius; j <= radius; j++ )
2848         {
2849             double r = std::sqrt((double)i * i + (double)j * j);
2850             if ( r > radius )
2851                 continue;
2852             space_weight[maxk] = (float)std::exp(r * r * gauss_space_coeff);
2853             space_ofs[maxk++] = (int)(i * temp.step + j * cn);
2854         }
2855
2856     char cvt[3][40];
2857     String cnstr = cn > 1 ? format("%d", cn) : "";
2858     String kernelName("bilateral");
2859     size_t sizeDiv = 1;
2860     if ((ocl::Device::getDefault().isIntel()) &&
2861         (ocl::Device::getDefault().type() == ocl::Device::TYPE_GPU))
2862     {
2863             //Intel GPU
2864             if (dst.cols % 4 == 0 && cn == 1) // For single channel x4 sized images.
2865             {
2866                 kernelName = "bilateral_float4";
2867                 sizeDiv = 4;
2868             }
2869      }
2870      ocl::Kernel k(kernelName.c_str(), ocl::imgproc::bilateral_oclsrc,
2871             format("-D radius=%d -D maxk=%d -D cn=%d -D int_t=%s -D uint_t=uint%s -D convert_int_t=%s"
2872             " -D uchar_t=%s -D float_t=%s -D convert_float_t=%s -D convert_uchar_t=%s -D gauss_color_coeff=%f",
2873             radius, maxk, cn, ocl::typeToStr(CV_32SC(cn)), cnstr.c_str(),
2874             ocl::convertTypeStr(CV_8U, CV_32S, cn, cvt[0]),
2875             ocl::typeToStr(type), ocl::typeToStr(CV_32FC(cn)),
2876             ocl::convertTypeStr(CV_32S, CV_32F, cn, cvt[1]),
2877             ocl::convertTypeStr(CV_32F, CV_8U, cn, cvt[2]), gauss_color_coeff));
2878     if (k.empty())
2879         return false;
2880
2881     Mat mspace_weight(1, d * d, CV_32FC1, space_weight);
2882     Mat mspace_ofs(1, d * d, CV_32SC1, space_ofs);
2883     UMat ucolor_weight, uspace_weight, uspace_ofs;
2884
2885     mspace_weight.copyTo(uspace_weight);
2886     mspace_ofs.copyTo(uspace_ofs);
2887
2888     k.args(ocl::KernelArg::ReadOnlyNoSize(temp), ocl::KernelArg::WriteOnly(dst),
2889            ocl::KernelArg::PtrReadOnly(uspace_weight),
2890            ocl::KernelArg::PtrReadOnly(uspace_ofs));
2891
2892     size_t globalsize[2] = { dst.cols / sizeDiv, dst.rows };
2893     return k.run(2, globalsize, NULL, false);
2894 }
2895
2896 #endif
2897 static void
2898 bilateralFilter_8u( const Mat& src, Mat& dst, int d,
2899     double sigma_color, double sigma_space,
2900     int borderType )
2901 {
2902     int cn = src.channels();
2903     int i, j, maxk, radius;
2904     Size size = src.size();
2905
2906     CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) && src.data != dst.data );
2907
2908     if( sigma_color <= 0 )
2909         sigma_color = 1;
2910     if( sigma_space <= 0 )
2911         sigma_space = 1;
2912
2913     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
2914     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
2915
2916     if( d <= 0 )
2917         radius = cvRound(sigma_space*1.5);
2918     else
2919         radius = d/2;
2920     radius = MAX(radius, 1);
2921     d = radius*2 + 1;
2922
2923     Mat temp;
2924     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
2925
2926 #if defined HAVE_IPP && (IPP_VERSION_MAJOR >= 7) && 0
2927     CV_IPP_CHECK()
2928     {
2929         if( cn == 1 )
2930         {
2931             bool ok;
2932             IPPBilateralFilter_8u_Invoker body(temp, dst, sigma_color * sigma_color, sigma_space * sigma_space, radius, &ok );
2933             parallel_for_(Range(0, dst.rows), body, dst.total()/(double)(1<<16));
2934             if( ok )
2935             {
2936                 CV_IMPL_ADD(CV_IMPL_IPP|CV_IMPL_MT);
2937                 return;
2938             }
2939             setIppErrorStatus();
2940         }
2941     }
2942 #endif
2943
2944     std::vector<float> _color_weight(cn*256);
2945     std::vector<float> _space_weight(d*d);
2946     std::vector<int> _space_ofs(d*d);
2947     float* color_weight = &_color_weight[0];
2948     float* space_weight = &_space_weight[0];
2949     int* space_ofs = &_space_ofs[0];
2950
2951     // initialize color-related bilateral filter coefficients
2952
2953     for( i = 0; i < 256*cn; i++ )
2954         color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);
2955
2956     // initialize space-related bilateral filter coefficients
2957     for( i = -radius, maxk = 0; i <= radius; i++ )
2958     {
2959         j = -radius;
2960
2961         for( ; j <= radius; j++ )
2962         {
2963             double r = std::sqrt((double)i*i + (double)j*j);
2964             if( r > radius )
2965                 continue;
2966             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
2967             space_ofs[maxk++] = (int)(i*temp.step + j*cn);
2968         }
2969     }
2970
2971     BilateralFilter_8u_Invoker body(dst, temp, radius, maxk, space_ofs, space_weight, color_weight);
2972     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
2973 }
2974
2975
2976 class BilateralFilter_32f_Invoker :
2977     public ParallelLoopBody
2978 {
2979 public:
2980
2981     BilateralFilter_32f_Invoker(int _cn, int _radius, int _maxk, int *_space_ofs,
2982         const Mat& _temp, Mat& _dest, float _scale_index, float *_space_weight, float *_expLUT) :
2983         cn(_cn), radius(_radius), maxk(_maxk), space_ofs(_space_ofs),
2984         temp(&_temp), dest(&_dest), scale_index(_scale_index), space_weight(_space_weight), expLUT(_expLUT)
2985     {
2986     }
2987
2988     virtual void operator() (const Range& range) const
2989     {
2990         int i, j, k;
2991         Size size = dest->size();
2992         #if CV_SSE3
2993         int CV_DECL_ALIGNED(16) idxBuf[4];
2994         float CV_DECL_ALIGNED(16) bufSum32[4];
2995         static const int CV_DECL_ALIGNED(16) bufSignMask[] = { 0x80000000, 0x80000000, 0x80000000, 0x80000000 };
2996         bool haveSSE3 = checkHardwareSupport(CV_CPU_SSE3);
2997         #endif
2998
2999         for( i = range.start; i < range.end; i++ )
3000         {
3001             const float* sptr = temp->ptr<float>(i+radius) + radius*cn;
3002             float* dptr = dest->ptr<float>(i);
3003
3004             if( cn == 1 )
3005             {
3006                 for( j = 0; j < size.width; j++ )
3007                 {
3008                     float sum = 0, wsum = 0;
3009                     float val0 = sptr[j];
3010                     k = 0;
3011                     #if CV_SSE3
3012                     if( haveSSE3 )
3013                     {
3014                         __m128 psum = _mm_setzero_ps();
3015                         const __m128 _val0 = _mm_set1_ps(sptr[j]);
3016                         const __m128 _scale_index = _mm_set1_ps(scale_index);
3017                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
3018
3019                         for( ; k <= maxk - 4 ; k += 4 )
3020                         {
3021                             __m128 _sw    = _mm_loadu_ps(space_weight + k);
3022                             __m128 _val   = _mm_set_ps(sptr[j + space_ofs[k+3]], sptr[j + space_ofs[k+2]],
3023                                                        sptr[j + space_ofs[k+1]], sptr[j + space_ofs[k]]);
3024                             __m128 _alpha = _mm_mul_ps(_mm_andnot_ps( _signMask, _mm_sub_ps(_val,_val0)), _scale_index);
3025
3026                             __m128i _idx = _mm_cvtps_epi32(_alpha);
3027                             _mm_store_si128((__m128i*)idxBuf, _idx);
3028                             _alpha = _mm_sub_ps(_alpha, _mm_cvtepi32_ps(_idx));
3029
3030                             __m128 _explut  = _mm_set_ps(expLUT[idxBuf[3]], expLUT[idxBuf[2]],
3031                                                          expLUT[idxBuf[1]], expLUT[idxBuf[0]]);
3032                             __m128 _explut1 = _mm_set_ps(expLUT[idxBuf[3]+1], expLUT[idxBuf[2]+1],
3033                                                          expLUT[idxBuf[1]+1], expLUT[idxBuf[0]+1]);
3034
3035                             __m128 _w = _mm_mul_ps(_sw, _mm_add_ps(_explut, _mm_mul_ps(_alpha, _mm_sub_ps(_explut1, _explut))));
3036                             _val = _mm_mul_ps(_w, _val);
3037
3038                              _sw = _mm_hadd_ps(_w, _val);
3039                              _sw = _mm_hadd_ps(_sw, _sw);
3040                              psum = _mm_add_ps(_sw, psum);
3041                         }
3042                         _mm_storel_pi((__m64*)bufSum32, psum);
3043
3044                         sum = bufSum32[1];
3045                         wsum = bufSum32[0];
3046                     }
3047                     #endif
3048
3049                     for( ; k < maxk; k++ )
3050                     {
3051                         float val = sptr[j + space_ofs[k]];
3052                         float alpha = (float)(std::abs(val - val0)*scale_index);
3053                         int idx = cvFloor(alpha);
3054                         alpha -= idx;
3055                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
3056                         sum += val*w;
3057                         wsum += w;
3058                     }
3059                     dptr[j] = (float)(sum/wsum);
3060                 }
3061             }
3062             else
3063             {
3064                 CV_Assert( cn == 3 );
3065                 for( j = 0; j < size.width*3; j += 3 )
3066                 {
3067                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
3068                     float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
3069                     k = 0;
3070                     #if  CV_SSE3
3071                     if( haveSSE3 )
3072                     {
3073                         __m128 sum = _mm_setzero_ps();
3074                         const __m128 _b0 = _mm_set1_ps(b0);
3075                         const __m128 _g0 = _mm_set1_ps(g0);
3076                         const __m128 _r0 = _mm_set1_ps(r0);
3077                         const __m128 _scale_index = _mm_set1_ps(scale_index);
3078                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
3079
3080                         for( ; k <= maxk-4; k += 4 )
3081                         {
3082                             __m128 _sw = _mm_loadu_ps(space_weight + k);
3083
3084                             const float* const sptr_k0 = sptr + j + space_ofs[k];
3085                             const float* const sptr_k1 = sptr + j + space_ofs[k+1];
3086                             const float* const sptr_k2 = sptr + j + space_ofs[k+2];
3087                             const float* const sptr_k3 = sptr + j + space_ofs[k+3];
3088
3089                             __m128 _b = _mm_loadu_ps(sptr_k0);
3090                             __m128 _g = _mm_loadu_ps(sptr_k1);
3091                             __m128 _r = _mm_loadu_ps(sptr_k2);
3092                             __m128 _z = _mm_loadu_ps(sptr_k3);
3093                             _MM_TRANSPOSE4_PS(_b, _g, _r, _z);
3094
3095                             __m128 _bt = _mm_andnot_ps(_signMask,_mm_sub_ps(_b,_b0));
3096                             __m128 _gt = _mm_andnot_ps(_signMask,_mm_sub_ps(_g,_g0));
3097                             __m128 _rt = _mm_andnot_ps(_signMask,_mm_sub_ps(_r,_r0));
3098
3099                             __m128 _alpha = _mm_mul_ps(_scale_index, _mm_add_ps(_rt,_mm_add_ps(_bt, _gt)));
3100
3101                             __m128i _idx  = _mm_cvtps_epi32(_alpha);
3102                             _mm_store_si128((__m128i*)idxBuf, _idx);
3103                             _alpha = _mm_sub_ps(_alpha, _mm_cvtepi32_ps(_idx));
3104
3105                             __m128 _explut  = _mm_set_ps(expLUT[idxBuf[3]], expLUT[idxBuf[2]], expLUT[idxBuf[1]], expLUT[idxBuf[0]]);
3106                             __m128 _explut1 = _mm_set_ps(expLUT[idxBuf[3]+1], expLUT[idxBuf[2]+1], expLUT[idxBuf[1]+1], expLUT[idxBuf[0]+1]);
3107
3108                             __m128 _w = _mm_mul_ps(_sw, _mm_add_ps(_explut, _mm_mul_ps(_alpha, _mm_sub_ps(_explut1, _explut))));
3109
3110                             _b = _mm_mul_ps(_b, _w);
3111                             _g = _mm_mul_ps(_g, _w);
3112                             _r = _mm_mul_ps(_r, _w);
3113
3114                              _w = _mm_hadd_ps(_w, _b);
3115                              _g = _mm_hadd_ps(_g, _r);
3116
3117                              _w = _mm_hadd_ps(_w, _g);
3118                              sum = _mm_add_ps(sum, _w);
3119                         }
3120                         _mm_store_ps(bufSum32, sum);
3121                         wsum  = bufSum32[0];
3122                         sum_b = bufSum32[1];
3123                         sum_g = bufSum32[2];
3124                         sum_r = bufSum32[3];
3125                     }
3126                     #endif
3127
3128                     for(; k < maxk; k++ )
3129                     {
3130                         const float* sptr_k = sptr + j + space_ofs[k];
3131                         float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
3132                         float alpha = (float)((std::abs(b - b0) +
3133                             std::abs(g - g0) + std::abs(r - r0))*scale_index);
3134                         int idx = cvFloor(alpha);
3135                         alpha -= idx;
3136                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
3137                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
3138                         wsum += w;
3139                     }
3140                     wsum = 1.f/wsum;
3141                     b0 = sum_b*wsum;
3142                     g0 = sum_g*wsum;
3143                     r0 = sum_r*wsum;
3144                     dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
3145                 }
3146             }
3147         }
3148     }
3149
3150 private:
3151     int cn, radius, maxk, *space_ofs;
3152     const Mat* temp;
3153     Mat *dest;
3154     float scale_index, *space_weight, *expLUT;
3155 };
3156
3157
3158 static void
3159 bilateralFilter_32f( const Mat& src, Mat& dst, int d,
3160                      double sigma_color, double sigma_space,
3161                      int borderType )
3162 {
3163     int cn = src.channels();
3164     int i, j, maxk, radius;
3165     double minValSrc=-1, maxValSrc=1;
3166     const int kExpNumBinsPerChannel = 1 << 12;
3167     int kExpNumBins = 0;
3168     float lastExpVal = 1.f;
3169     float len, scale_index;
3170     Size size = src.size();
3171
3172     CV_Assert( (src.type() == CV_32FC1 || src.type() == CV_32FC3) && src.data != dst.data );
3173
3174     if( sigma_color <= 0 )
3175         sigma_color = 1;
3176     if( sigma_space <= 0 )
3177         sigma_space = 1;
3178
3179     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
3180     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
3181
3182     if( d <= 0 )
3183         radius = cvRound(sigma_space*1.5);
3184     else
3185         radius = d/2;
3186     radius = MAX(radius, 1);
3187     d = radius*2 + 1;
3188     // compute the min/max range for the input image (even if multichannel)
3189
3190     minMaxLoc( src.reshape(1), &minValSrc, &maxValSrc );
3191     if(std::abs(minValSrc - maxValSrc) < FLT_EPSILON)
3192     {
3193         src.copyTo(dst);
3194         return;
3195     }
3196
3197     // temporary copy of the image with borders for easy processing
3198     Mat temp;
3199     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
3200     const double insteadNaNValue = -5. * sigma_color;
3201     patchNaNs( temp, insteadNaNValue ); // this replacement of NaNs makes the assumption that depth values are nonnegative
3202                                         // TODO: make insteadNaNValue avalible in the outside function interface to control the cases breaking the assumption
3203     // allocate lookup tables
3204     std::vector<float> _space_weight(d*d);
3205     std::vector<int> _space_ofs(d*d);
3206     float* space_weight = &_space_weight[0];
3207     int* space_ofs = &_space_ofs[0];
3208
3209     // assign a length which is slightly more than needed
3210     len = (float)(maxValSrc - minValSrc) * cn;
3211     kExpNumBins = kExpNumBinsPerChannel * cn;
3212     std::vector<float> _expLUT(kExpNumBins+2);
3213     float* expLUT = &_expLUT[0];
3214
3215     scale_index = kExpNumBins/len;
3216
3217     // initialize the exp LUT
3218     for( i = 0; i < kExpNumBins+2; i++ )
3219     {
3220         if( lastExpVal > 0.f )
3221         {
3222             double val =  i / scale_index;
3223             expLUT[i] = (float)std::exp(val * val * gauss_color_coeff);
3224             lastExpVal = expLUT[i];
3225         }
3226         else
3227             expLUT[i] = 0.f;
3228     }
3229
3230     // initialize space-related bilateral filter coefficients
3231     for( i = -radius, maxk = 0; i <= radius; i++ )
3232         for( j = -radius; j <= radius; j++ )
3233         {
3234             double r = std::sqrt((double)i*i + (double)j*j);
3235             if( r > radius )
3236                 continue;
3237             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
3238             space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
3239         }
3240
3241     // parallel_for usage
3242
3243     BilateralFilter_32f_Invoker body(cn, radius, maxk, space_ofs, temp, dst, scale_index, space_weight, expLUT);
3244     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
3245 }
3246
3247 }
3248
3249 void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,
3250                       double sigmaColor, double sigmaSpace,
3251                       int borderType )
3252 {
3253     _dst.create( _src.size(), _src.type() );
3254
3255     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
3256                ocl_bilateralFilter_8u(_src, _dst, d, sigmaColor, sigmaSpace, borderType))
3257
3258     Mat src = _src.getMat(), dst = _dst.getMat();
3259
3260     if( src.depth() == CV_8U )
3261         bilateralFilter_8u( src, dst, d, sigmaColor, sigmaSpace, borderType );
3262     else if( src.depth() == CV_32F )
3263         bilateralFilter_32f( src, dst, d, sigmaColor, sigmaSpace, borderType );
3264     else
3265         CV_Error( CV_StsUnsupportedFormat,
3266         "Bilateral filtering is only implemented for 8u and 32f images" );
3267 }
3268
3269 //////////////////////////////////////////////////////////////////////////////////////////
3270
3271 CV_IMPL void
3272 cvSmooth( const void* srcarr, void* dstarr, int smooth_type,
3273           int param1, int param2, double param3, double param4 )
3274 {
3275     cv::Mat src = cv::cvarrToMat(srcarr), dst0 = cv::cvarrToMat(dstarr), dst = dst0;
3276
3277     CV_Assert( dst.size() == src.size() &&
3278         (smooth_type == CV_BLUR_NO_SCALE || dst.type() == src.type()) );
3279
3280     if( param2 <= 0 )
3281         param2 = param1;
3282
3283     if( smooth_type == CV_BLUR || smooth_type == CV_BLUR_NO_SCALE )
3284         cv::boxFilter( src, dst, dst.depth(), cv::Size(param1, param2), cv::Point(-1,-1),
3285             smooth_type == CV_BLUR, cv::BORDER_REPLICATE );
3286     else if( smooth_type == CV_GAUSSIAN )
3287         cv::GaussianBlur( src, dst, cv::Size(param1, param2), param3, param4, cv::BORDER_REPLICATE );
3288     else if( smooth_type == CV_MEDIAN )
3289         cv::medianBlur( src, dst, param1 );
3290     else
3291         cv::bilateralFilter( src, dst, param1, param3, param4, cv::BORDER_REPLICATE );
3292
3293     if( dst.data != dst0.data )
3294         CV_Error( CV_StsUnmatchedFormats, "The destination image does not have the proper type" );
3295 }
3296
3297 /* End of file. */